The Dicke Model

Inspirations taken from this QuTiP tutorial by J. R. Johansson.

This tutorial mainly demonstrates the use of

• jmat
• plot_wigner
• plot_fock_distribution
• entropy_mutual

to explore the Dicke model.

Introduction

The Dicke model describes a system where $N$ two-level atoms interact with a single quantized electromagnetic mode. The original microscopic form of the Dicke Hamiltonian is given by:

$${\hat{H}}_D=\omega {\hat{a}}^{†}\hat{a}+\sum _{i=1}^N\frac{{\omega }_0}{2}{\hat{\sigma }}_z^{(i)}+\sum _{i=1}^N\frac{g}{\sqrt{N}}(\hat{a}+{\hat{a}}^{†})({\hat{\sigma }}_{+}^{(i)}+{\hat{\sigma }}_{-}^{(i)})$$

where:

• $\hat{a}$ and ${\hat{a}}^{†}$ are the cavity (with frequency $\omega$) annihilation and creation operators, respectively.
• ${\hat{\sigma }}_z^{(i)},{\hat{\sigma }}_{\pm }^{(i)}$ are the Pauli matrices for the $i$-th two-level atom, with transition frequency ${\omega }_0$.
• $g$ represents the light-matter coupling strength.

This formulation explicitly treats each spin individually. However, when the atoms interact identically with the cavity, we can rewrite the Hamiltonian regarding collective spin operators.

$${\hat{J}}_z=\sum _{i=1}^N{\hat{\sigma }}_z^{(i)},{\hat{J}}_{\pm }=\sum _{i=1}^N{\hat{\sigma }}_{\pm }^{(i)},$$ which describe the total spin of the system as a pseudospin of length $j=N/2$. Using these collective operators, the reformulated Dicke Hamiltonian takes the form:

$${\hat{H}}_D=\omega {\hat{a}}^{†}\hat{a}+{\omega }_0{\hat{J}}_z+\frac{g}{\sqrt{N}}(\hat{a}+{\hat{a}}^{†})({\hat{J}}_{+}+{\hat{J}}_{-})$$

This formulation reduces complexity, as it allows us to work on a collective basis instead of the whole individual spin Hilbert space. The superradiant phase transition occurs when the coupling strength $g$ exceeds a critical threshold $g_c=0.5\sqrt{\omega /{\omega }_0}$, leading to a macroscopic population of the cavity mode.

Code demonstration

using QuantumToolbox
using CairoMakie

CairoMakie.activate!(type = "svg")

ω = 1
ω0 = 1

gc = √(ω/ω0)/2

κ = 0.05;

Here, we define some functions for later usage.

# M: cavity Hilbert space truncation, N: number of atoms
Jz(M, N) = (qeye(M) ⊗ jmat(N/2, :z))

a(M, N) = (destroy(M) ⊗ qeye(N+1))

function H(M, N, g)
    j = N / 2
    n = 2 * j + 1

    a_ = a(M, N)
    Jz_ = Jz(M, N)
    Jp = qeye(M) ⊗ jmat(j, :+)
    Jm = qeye(M) ⊗ jmat(j, :-);

    H0 = ω * a_' * a_ + ω0 * Jz_
    H1 = 1/ √N * (a_ + a_') * (Jp + Jm)

    return (H0 + g * H1)
end;

Take the example of 4 atoms.

M0, N0 = 16, 4

H0(g) = H(M0, N0, g)

a0 = a(M0, N0)
Jz0 = Jz(M0, N0);

gs = 0.0:0.05:1.0
ψGs = map(gs) do g
    H = H0(g)
    vals, vecs = eigenstates(H)
    vecs[1]
end

nvec = expect(a0'*a0, ψGs)
Jzvec = expect(Jz0, ψGs);

fig = Figure(size=(900, 350))
axn = Axis(
 fig[1,1],
 xlabel = "interaction strength",
 ylabel = L"\langle \hat{n} \rangle"
)
axJz = Axis(
 fig[1,2],
 xlabel = "interaction strength",
 ylabel = L"\langle \hat{J}_{z} \rangle"
)
ylims!(-N0/2, N0/2)
lines!(axn, gs, real(nvec))
lines!(axJz, gs, real(Jzvec))
fig

The expectation value of photon number and ${\hat{J}}_z$ showed a sudden increment around $g_c$.

# the indices in coupling strength list (gs)
# to display wigner and fock distribution
cases = 1:5:21

fig = Figure(size = (900,650))
for (hpos, idx) in enumerate(cases)
    g = gs[idx] # coupling strength
    ρcav = ptrace(ψGs[idx], 1) # cavity reduced state
    
    # plot wigner
    _, ax, hm = plot_wigner(ρcav, location = fig[1,hpos])
    ax.title = "g = $g"
    ax.aspect = 1
    
    # plot fock distribution
    _, ax2 = plot_fock_distribution(ρcav, location = fig[2,hpos])

    ax2.xticks = (0:2:size(ρcav,1)-1, string.(0:2:size(ρcav,1)-1))
    
    if hpos != 1
        hideydecorations!(ax, ticks=false)
        hideydecorations!(ax2, ticks=false)
        if hpos == 5 # Add colorbar with the last returned heatmap (_hm) 
            Colorbar(fig[1,6], hm)
        end
    end    
end

# plot average photon number with respect to coupling strength
ax3 = Axis(fig[3,1:6], height=200, xlabel=L"g", ylabel=L"\langle \hat{n} \rangle")
xlims!(ax3, -0.02, 1.02)
lines!(ax3, gs, real(nvec), color=:teal)
ax3.xlabelsize, ax3.ylabelsize = 20, 20
vlines!(ax3, gs[cases], color=:orange, linestyle = :dash, linewidth = 4)

fig

As $g$ increases, the cavity ground state’s wigner function plot looks more coherent than a thermal state.

Ns = 8:8:24
slists = map(Ns) do N
    slist = map(gs) do g
        H_ = H(M0, N, g)
        _, vecs = eigenstates(H_)
        entropy_mutual(vecs[1], 1, 2)
    end
end;

fig = Figure(size=(900, 400))
ax = Axis(fig[1,1])
ax.xlabel = "coupling strength"
ax.ylabel = "mutual entropy"

for (idx, slist) in enumerate(slists)
    lines!(gs, slist, label = string(Ns[idx]))
end

Legend(fig[1,2], ax, label = "number of atoms")
fig

We further consult mutual entropy between the cavity and the spins as a measure of their correlation; the result showed that as the number of atoms $N$ increases, the peak of mutual entropy moves closer to $g_c$.

Version Information

QuantumToolbox.versioninfo()

 QuantumToolbox.jl: Quantum Toolbox in Julia
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
Copyright © QuTiP team 2022 and later.
Current admin team:
    Alberto Mercurio and Yi-Te Huang

Package information:
====================================
Julia              Ver. 1.11.5
QuantumToolbox     Ver. 0.31.1
SciMLOperators     Ver. 0.3.13
LinearSolve        Ver. 3.17.0
OrdinaryDiffEqCore Ver. 1.26.0

System information:
====================================
OS       : Linux (x86_64-linux-gnu)
CPU      : 4 × AMD EPYC 7763 64-Core Processor
Memory   : 15.621 GB
WORD_SIZE: 64
LIBM     : libopenlibm
LLVM     : libLLVM-16.0.6 (ORCJIT, znver3)
BLAS     : libopenblas64_.so (ilp64)
Threads  : 4 (on 4 virtual cores)