GitHub

Mean-field Two-sided Cavity

Here we show how to solve the dynamics of the example Two-sided Cavity with Atoms in the Heisenberg picture with a higher-order mean-field approach (cumulant expansion), which is done with the package QuantumCumulants.jl.

using QuantumInputOutput
using SecondQuantizedAlgebra
using QuantumCumulants
using ModelingToolkit
using OrdinaryDiffEq
using QuantumOpticsBase
using Plots
@rnumbers E κ_L κ_R Δ g γ
Natoms = 2

hc = FockSpace(:cavity)
ha(i) = NLevelSpace("a_$i", 2)
h = hc ⊗ tensor([ha(i) for i = 1:Natoms]...);

a = Destroy(h, :a, 1) # cavity
σ(α, i, j) = Transition(h, "σ_$(α)", i, j, 1+α) # two-level atom α
∑σ(i, j) = sum(σ(α, i, j) for α = 1:Natoms) # collective atomic operator

Empty two-sided cavity

We couple a classical drive into the cavity through the left mirror $(\kappa_L)$. The decay through the right mirror can be added in several ways: with the concatenate rule, including it already in the initial cavity SLH triple or by simply including the decay term to the Lindblad by hand. In this example, we use the first option.

G_d = SLH(1, E, 0) # classical drive
H_cavity = -Δ*a'a
G_c_L = SLH(1, [√(κ_L)*a], H_cavity)

G_cav_L_drive = G_d ▷ G_c_L
G_c_R = SLH(1, [√(κ_R)*a], 0)

G_cav_L_R_drive = G_cav_L_drive ⊞ G_c_R

Note that one needs to be careful to not double-count the Hamiltonian terms with the concatenation rule.

H1 = get_hamiltonian(G_cav_L_R_drive)

\[-0.5 i \sqrt{\kappa_{L}} E a^\dagger + 0.5 i \sqrt{\kappa_{L}} E a -1 \Delta a^\dagger a\]

L1_L = get_lindblad(G_cav_L_R_drive)[1]

\[\sqrt{\kappa_{L}} a + E\]

L1_R = get_lindblad(G_cav_L_R_drive)[2]

\[\sqrt{\kappa_{R}} a\]

The typical cavity drive-term $\sqrt{\kappa_L} E (a^\dagger + a)$ is a combination of Hamiltonian term and Lindblad. We use the function meanfield to obtain the equation for the intra-cavity field, which leads to a closed set of equations in this particular case.

eqs_a = meanfield([a], H1, [L1_L, L1_R])

\[\begin{align} \frac{d}{dt} \langle a\rangle &= 1 i \Delta \langle a\rangle -1.0 \sqrt{\kappa{L}} E -0.5 \left( \left( \sqrt{\kappa{L}} \right)^{2} + \left( \sqrt{\kappa_{R}} \right)^{2} \right) \langle a\rangle \end{align}\]

We defined the numerical parameters and the initial state of the system, create the ODE problem and solve the dynamics.

# numerical parameters
En = 0.5
κ_Rn = 1.0
κ_Ln = 1.0
Δn = 0.0
Δn_ls = [-5.0:0.1:5.0;];
lΔ=length(Δn_ls)

p_sym = [E, κ_R, κ_L, Δ]
p_num = [En, κ_Rn, κ_Ln, Δn]
dict_p1 = Dict(p_sym .=> p_num)

# initial state
u0 = [0.0im]
# numerical system
T = [0:0.01:1;]*20
@named sys_a = System(eqs_a)
dict = merge(Dict(p_sym .=> p_num), Dict(unknowns(sys_a) .=> u0))
prob_a = ODEProblem(sys_a, dict, (0.0, T[end]))
sol_a = solve(prob_a, Tsit5(); saveat = T)
n_cavity = abs2.(get_solution(sol_a, a))
n_ref = abs2.(get_solution(sol_a, √(κ_Ln)*a) .+ En)
n_trans = abs2.(get_solution(sol_a, √(κ_Rn)*a))
p1 = plot(T, n_cavity; label = "", xlabel = "t", ylabel = "cavity photons", grid = true)
p2 = plot(T, n_ref; label = "reflection")
plot!(p2, T, n_trans; label = "transmission", ls = :dash)
plot!(p2; xlabel = "t", ylabel = "intensity rate", grid = true, legend = :best)
plot(p1, p2; layout = (2, 1), size = (600, 500))
Example block output

Now we scan the laser-cavity detuning $\Delta$ to plot the transmission and reflection spectrum.

dict_p_Δ(Δn) = merge(Dict(p_sym .=> [En, κ_Rn, κ_Ln, Δn]), Dict(unknowns(sys_a) .=> u0))
n_ref_Δ = zeros(lΔ)
n_trans_Δ = zeros(lΔ)

for it = 1:lΔ
    Δn_ = Δn_ls[it]
    prob_ = ODEProblem(sys_a, dict_p_Δ(Δn_), (0.0, T[end]))
    sol_ = solve(prob_, Tsit5(); saveat = T)

    n_ref_Δ[it] = abs2.(get_solution(sol_, √(κ_Ln)*a) .+ En)[end]
    n_trans_Δ[it] = abs2.(get_solution(sol_, √(κ_Ln)*a))[end]
end
p = plot(Δn_ls, n_ref_Δ; label = "reflection")
plot!(p, Δn_ls, n_trans_Δ; label = "transmission", ls = :dash)
plot!(p; xlabel = "Δ", grid = true, legend = :best, size = (600, 350))
p
Example block output

Two-sided cavity with atoms

In the following, we include $N=2$ two-level atoms in the cavity and simulate the transmission and reflection of a coherent Gaussian pulse with a mean photon number of $|\alpha|^2 = 1/10$. We assume that the atoms are on resonance with the cavity, i.e. $\Delta = \Delta_c = \Delta_a$.

@syms t::Real
@register_symbolic Et(t)

G_d_t = SLH(1, Et(t), 0)
H_ac = -Δ*(a'a + ∑σ(2, 2)) + g*(a'∑σ(1, 2) + a*∑σ(2, 1))
G_ac = SLH(1, √κ_L*a, H_ac)
G_ac_drive = (G_d_t ▷ G_ac) ⊞ SLH(1, √κ_R*a, 0)
H2 = G_ac_drive.hamiltonian

\[g a^\dagger {\sigma_2}^{{12}} + g a {\sigma_2}^{{21}} + g a^\dagger {\sigma_1}^{{12}} + g a {\sigma_1}^{{21}} + 0.5 i \sqrt{\kappa_{L}} \mathrm{Et}\left( t \right) a -0.5 i \sqrt{\kappa_{L}} \mathrm{Et}\left( t \right) a^\dagger -1 \Delta a^\dagger a -1 \Delta {\sigma_2}^{{22}} -1 \Delta {\sigma_1}^{{22}}\]

L2_L = G_ac_drive.lindblad[1]

\[\sqrt{\kappa_{L}} a + \mathrm{Et}\left( t \right)\]

L2_R = G_ac_drive.lindblad[2]

\[\sqrt{\kappa_{R}} a\]

We derive the equations of motion for system with a second-order mean-field approximation.

J_add = [√(γ)*σ(α, 1, 2) for α = 1:Natoms]
eqs2 = meanfield([a'a, σ(1, 2, 2)], H2, [L2_L, L2_R, J_add...]; order = 2)

\[\begin{align} \frac{d}{dt} \langle a^\dagger a\rangle &= -1 i g \langle a^\dagger {\sigma1}^{{12}}\rangle + \langle a^\dagger a\rangle \left( -1.0 \left( \sqrt{\kappa{L}} \right)^{2} -1.0 \left( \sqrt{\kappa{R}} \right)^{2} \right) -1 i g \langle a^\dagger {\sigma2}^{{12}}\rangle + 1 i \langle a {\sigma2}^{{21}}\rangle g -1.0 \sqrt{\kappa{L}} \langle a\rangle \mathrm{Et}\left( t \right) -1.0 \sqrt{\kappa{L}} \langle a^\dagger\rangle \mathrm{Et}\left( t \right) + 1 i g \langle a {\sigma1}^{{21}}\rangle \\ \frac{d}{dt} \langle {\sigma1}^{{22}}\rangle &= 1 i g \langle a^\dagger {\sigma1}^{{12}}\rangle -1.0 \left( \sqrt{\gamma} \right)^{2} \langle {\sigma1}^{{22}}\rangle -1 i g \langle a {\sigma1}^{{21}}\rangle \end{align}\]

eqs2_c = complete(eqs2)
length(eqs2_c)
18

Again, we defined the numerical parameters and the initial state of the system, create the ODE problem and solve the dynamics.

# numerical parameter
κ_Rn2 = κ_Ln2 = 2π*1
gn2 = 2π*0.4
Δn2 = 0.0
γn = 2π*0.1

p_sym2 = [κ_R, κ_L, Δ, g, γ]
p_num2 = [κ_Rn2, κ_Ln2, Δn2, gn2, γn]

σp = 10/κ_Ln2 # pulse width
Tp = 4σp # pulse peak
Tend = 3Tp
α0 = √(0.1) # √ of total photon number
Ω0 = α0*2*√(κ_Ln2)/(π^(1/4)*√(σp))
Ω1(t) = Ω0/2*exp(-(t-Tp)^2 / (2*σp^2))
Et(t) = Ω1(t)/√(κ_Ln2)

T2 = [0:0.001:1;]*Tend
ΔT = T2[2] - T2[1]
n_pulse = round(sum(abs2.(Et.(T2)))*ΔT, digits = 7)
@show n_pulse

dict_p2 = Dict(p_sym2 .=> p_num2)
n_pulse = 0.1
# initial state
bc1 = FockBasis(4)
ba = NLevelBasis(2)
b = bc1 ⊗ tensor([ba for i = 1:Natoms]...)
ψ0 = LazyKet(b, (fockstate(bc1, 0), [nlevelstate(ba, 1) for i = 1:Natoms]...))
u0_2 = initial_values(eqs2_c, ψ0) # initial state
@named sys2 = System(eqs2_c) # initial state
dict2 = merge(Dict(p_sym2 .=> p_num2), Dict(unknowns(sys2) .=> u0_2))
prob2 = ODEProblem(sys2, dict2, (0.0, T2[end]))
sol2 = solve(prob2, Tsit5(); saveat = T2)
n_ref2 = abs2.(get_solution(sol2, √(κ_Ln2)*a) + Et.(T2))
n_trans2 = abs2.(get_solution(sol2, √(κ_Rn2)*a))

p = plot(
    T2,
    n_trans2;
    label = "transmission = $(round(sum(n_trans2) * ΔT / n_pulse * 100))%",
)
plot!(
    p,
    T2,
    n_ref2;
    ls = :dash,
    label = "reflection = $(round(sum(n_ref2) * ΔT / n_pulse * 100))%",
)
plot!(p; xlabel = "t", legend = :best, grid = true, size = (600, 350))
p
Example block output

We can see that all results agree with the full quantum dynamics of the example Two-sided Cavity with Atoms, which means the second order cumulant expansion is good approximation here.

Package versions

These results were obtained using the following versions:

using InteractiveUtils
versioninfo()

using Pkg
Pkg.status(
    [
        "QuantumInputOutput",
        "SecondQuantizedAlgebra",
        "QuantumCumulants",
        "ModelingToolkit",
        "OrdinaryDiffEq",
        "QuantumOpticsBase",
        "Plots",
    ],
    mode = PKGMODE_MANIFEST,
)
Julia Version 1.12.5
Commit 5fe89b8ddc1 (2026-02-09 16:05 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 9V74 80-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 4 default, 1 interactive, 4 GC (on 4 virtual cores)
Environment:
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_DEBUG = Documenter
  JULIA_NUM_THREADS = auto
Status `~/work/QuantumInputOutput.jl/QuantumInputOutput.jl/docs/Manifest.toml`
⌅ [961ee093] ModelingToolkit v10.32.1
⌃ [1dea7af3] OrdinaryDiffEq v6.103.0
  [91a5bcdd] Plots v1.41.6
  [35bcea6d] QuantumCumulants v0.4.3
  [18f9eda6] QuantumInputOutput v0.1.0 `~/work/QuantumInputOutput.jl/QuantumInputOutput.jl`
  [4f57444f] QuantumOpticsBase v0.5.9
  [f7aa4685] SecondQuantizedAlgebra v0.4.4
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`

This page was generated using Literate.jl.