Cavity Antiresonance

In this example we investigate a system of $N$ closely spaced quantum emitters inside a coherently driven single mode cavity. The model is described in D. Plankensteiner, et. al., Phys. Rev. Lett. 119, 093601 (2017). The Hamiltonian of this system is composed of three parts $H = H_c + H_a + H_{\mathrm{int}}$, the driven cavity $H_c$, the dipole-dipole interacting atoms $H_a$ and the atom-cavity interaction $H_\mathrm{int}$:

\[\begin{align} H_\mathrm{c} &= \hbar \Delta_c a^\dagger a + \hbar \eta (a^\dagger + a) \\ &\\ H_a &= \hbar \Delta_a \sum\limits_{j} \sigma_j^{22} + \hbar \sum\limits_{i \neq j} \Omega_{ij} \sigma_i^{21} \sigma_j^{12} &\\ H_\mathrm{int} &= \hbar \sum\limits_{j} g_j (a^\dagger \sigma_j^{12} + a \sigma_j^{21}) \end{align}\]

Additionally the system features two decay channels, the lossy cavity with photon decay rate $\kappa$ and collective atomic emission described by the decay-rate matrix $\Gamma_{ij}$.

We start by loading the packages.

using QuantumCumulants
using OrdinaryDiffEqTsit5, ModelingToolkitBase
using Plots

The Hilbert space for this system is given by one cavity mode and $N$ two-level atoms. Here we use symbolic indices, sums and double sums to define the system. The parameters $g_j, \, \Gamma_{ij}$ and $\Omega_{ij}$ are defined as indexed variables of atom $i$ and $j$. We will describe the system in first order mean-field. Because $\Gamma_{ij}$ and $\Omega_{ij}$ depend on both atom labels, they enter the double sums as index-dependent coefficients; this per-pair dependence is carried through the cumulant expansion and evaluate without losing the summation scope.

hc = FockSpace(:cavity) # Hilbert space
ha = NLevelSpace(Symbol(:atom), 2)
h = hc ⊗ ha

@variables N Δc η Δₐ κ # Parameters
g(i) = IndexedVariable(:g, i)
Γ(i, j) = DoubleIndexedVariable(:Γ, i, j)
Ω(i, j) = DoubleIndexedVariable(:Ω, i, j; identical = false)


i = Index(h, :i, N, ha) # Indices
j = Index(h, :j, N, ha)

The kwarg ’identical=false’ for the double indexed variable specifies that $\Omega_{ij} = 0$ for $i = j$. Now we create the operators on the composite Hilbert space using the $\texttt{IndexedOperator}$ constructor, which assigns each $\texttt{Transition}$ operator an $\texttt{Index}$.

@qnumbers a::Destroy(h)
σ(x, y, k) = IndexedOperator(Transition(h, :σ, x, y), k)

We define the Hamiltonian and Liouvillian. For the collective atomic decay we write the corresponding dissipative processes with a double indexed variable $R_{ij}$ and an indexed jump operator $J_j$, such that an operator average $\langle \mathcal{O} \rangle$ follows the equation

\[\begin{equation} \langle \dot{\mathcal{O}} \rangle = \sum_{ij} R_{ij} \left( \langle J_i^\dagger \mathcal{O} J_j \rangle - \frac{1}{2} \langle J_i^\dagger J_j \mathcal{O} \rangle - \frac{1}{2} \langle \mathcal{O} J_i^\dagger J_j \rangle \right). \end{equation}\]

The inner dipole-dipole sum excludes the diagonal i == j by passing the non_equal vector [i] to the inner Σ (SQA v0.5 replaced the old non_equal=true keyword with this explicit form).

Hc = Δc * a'a + η * (a' + a) # Hamiltonian
Ha = Δₐ * Σ(σ(2, 2, i), i) + Σ(Σ(Ω(i, j) * σ(2, 1, i) * σ(1, 2, j), j, [i]), i)
Hi = Σ(g(i) * (a' * σ(1, 2, i) + a * σ(2, 1, i)), i)
H = Hc + Ha + Hi

J = [a, σ(1, 2, i)] # Jump operators
rates = [κ, Γ(i, j)]

We derive the system of equations in first order mean-field.

eqs = meanfield(a, H, J; rates = rates, order = 1)
complete!(eqs)

\[ \begin{aligned} \partial_{t} \langle a \rangle &= \underset{i}{\overset{N}{\sum}} - \langle {\sigma}_{i}^{{12}} \rangle ~ g\left( i \right) ~ i - i \eta + \langle a \rangle \left( - \frac{1}{2} \kappa - i \mathtt{{\Delta}c} \right) \\[-0.0em] \partial_{t} \langle {\sigma}_{i_{2}}^{{12}} \rangle &= \underset{j,i_2{\neq}j}{\overset{N}{\sum}}\langle {\sigma}_{j}^{{12}} \rangle ~ \left( - 0.5 ~ \Gamma\left( \mathtt{i_{2}}, j \right) -\mathit{i} ~ \Omega\left( \mathtt{i_{2}}, j \right) \right) + \langle {\sigma}_{i_{2}}^{{22}} \rangle \underset{j,i_2{\neq}j}{\overset{N}{\sum}}\langle {\sigma}_{j}^{{12}} \rangle ~ \left( \Gamma\left( \mathtt{i_{2}}, j \right) + 2.0 ~ i ~ \Omega\left( \mathtt{i_{2}}, j \right) \right) - \langle a \rangle g\left( \mathtt{i_{2}} \right) i + \langle {\sigma}_{i_{2}}^{{12}} \rangle \left( - 0.5 \Gamma\left( \mathtt{i_{2}}, \mathtt{i_{2}} \right) - i \mathtt{\Delta_a} \right) + 2 \langle {\sigma}_{i_{2}}^{{22}} \rangle \langle a \rangle g\left( \mathtt{i_{2}} \right) i \\[-0.0em] \partial_{t} \langle {\sigma}_{i_{2}}^{{22}} \rangle &= \langle {\sigma}_{i_{2}}^{{12}} \rangle \underset{i{\neq}i_2}{\overset{N}{\sum}}\langle {\sigma}_{i}^{{21}} \rangle ~ \left( - 0.5 ~ \Gamma\left( i, \mathtt{i_{2}} \right) + i ~ \Omega\left( i, \mathtt{i_{2}} \right) \right) + \langle {\sigma}_{i_{2}}^{{21}} \rangle \underset{j,i_2{\neq}j}{\overset{N}{\sum}}\langle {\sigma}_{j}^{{12}} \rangle ~ \left( - 0.5 ~ \Gamma\left( \mathtt{i_{2}}, j \right) - i ~ \Omega\left( \mathtt{i_{2}}, j \right) \right) + \langle {\sigma}_{i_{2}}^{{22}} \rangle \left( \mathrm{real}\left( - \Gamma\left( \mathtt{i_{2}}, \mathtt{i_{2}} \right) \right) + i \mathrm{imag}\left( - \Gamma\left( \mathtt{i_{2}}, \mathtt{i_{2}} \right) \right) \right) + \langle {\sigma}_{i_{2}}^{{12}} \rangle \langle a^{\dagger} \rangle g\left( \mathtt{i_{2}} \right) i - \langle a \rangle \langle {\sigma}_{i_{2}}^{{21}} \rangle g\left( \mathtt{i_{2}} \right) i \end{aligned} \]

To create the equations for a specific number of atoms we use the function evaluate.

N_ = 2
eqs_ = evaluate(eqs; limits = (N => N_))
sys = mtkcompile(System(eqs_; name = :sys))

Finally we need to define the initial state of the system and the numerical parameters. In the end we want to obtain the transmission rate $T$ of our system. For this purpose we calculate the steady state photon number in the cavity $|\langle a \rangle|^2$ for different laser frequencies.

u0 = zeros(ComplexF64, length(eqs_.states))
Γ_ = 1.0 # parameter
d = 2π * 0.08 #0.08λ
θ = π / 2

Ωij(i, j) =
    i == j ? 0 :
    Γ_ * (-3 / 4) * ((1 - (cos(θ))^2) * cos(d) / d - (1 - 3 * (cos(θ))^2) * (sin(d) / (d^2) + (cos(d) / (d^3))))
Γij(i, j) =
    i == j ? Γ_ :
    Γ_ * (3 / 2) * ((1 - (cos(θ))^2) * sin(d) / d + (1 - 3 * (cos(θ))^2) * ((cos(d) / (d^2)) - (sin(d) / (d^3))))

g_ = 2Γ_
κ_ = 20Γ_
Δₐ_ = 0Γ_
Δc_ = 0Γ_
η_ = κ_ / 100

gi_ = [g_ * (-1)^k for k in 1:N_]
Γij_ = [Γij(k, l) for k in 1:N_, l in 1:N_]
Ωij_ = [Ωij(k, l) for k in 1:N_, l in 1:N_]

Δ_ls = [-10:0.05:10;] * Γ_
n_ls = zeros(length(Δ_ls))
for k in eachindex(Δ_ls)
    Δc_i = Δ_ls[k]
    Δₐ_i = Δc_i + Ωij(1, 2)
    p = parameter_map(
        eqs_, Dict(
            Δc => Δc_i, η => η_, Δₐ => Δₐ_i, κ => κ_,
            g(i) => gi_, Γ(i, j) => Γij_, Ω(i, j) => Ωij_,
        )
    )
    prob_ = ODEProblem(sys, merge(initial_values(eqs_, u0), Dict(p)), (0.0, 20.0))
    sol_ = solve(prob_, Tsit5())
    n_ls[k] = abs2(get_solution(sol_, a, eqs_)(sol_.t[end]))
end

The transmission rate $T$ with respect to the pump laser detuning is given by the relative steady state intra-cavity photon number $n(\Delta)/n_\mathrm{max}$. We qualitatively reproduce the antiresonance from D. Plankensteiner, et. al., Phys. Rev. Lett. 119, 093601 (2017) for two atoms.

T = n_ls ./ maximum(n_ls)
plot(Δ_ls, T, xlabel = "Δ/Γ", ylabel = "T", legend = false)

Package versions

These results were obtained using the following versions:

using InteractiveUtils
versioninfo()

using Pkg
Pkg.status(
    ["QuantumCumulants", "OrdinaryDiffEqTsit5", "ModelingToolkitBase", "Plots"],
    mode = PKGMODE_MANIFEST,
)

Julia Version 1.13.1
Commit 96ca370cf0e (2026-09-25 19:34 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 9V45 96-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver5)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)
Environment:
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_DEBUG = Documenter
Status `~/work/QuantumCumulants.jl/QuantumCumulants.jl/docs/Manifest.toml`
  [47edcb42] ADTypes v1.24.0
  [1520ce14] AbstractTrees v0.4.5
  [4fba245c] ArrayInterface v7.30.2
  [aae01518] BandedMatrices v1.13.1
  [caf10ac8] BipartiteGraphs v0.1.14
  [8e7c35d0] BlockArrays v1.10.1
⌅ [861a8166] Combinatorics v1.0.2
  [38540f10] CommonSolve v0.2.14
  [34da2185] Compat v4.18.1
  [187b0558] ConstructionBase v1.6.0
  [d38c429a] Contour v0.6.3
  [864edb3b] DataStructures v0.19.6
  [2b5f629d] DiffEqBase v7.21.3
  [459566f4] DiffEqCallbacks v4.19.4
  [b552c78f] DiffRules v1.16.0
  [a0c0ee7d] DifferentiationInterface v0.7.21
  [ffbed154] DocStringExtensions v0.9.5
  [5b8099bc] DomainSets v0.8.3
  [4e289a0a] EnumX v1.0.7
  [f151be2c] EnzymeCore v0.8.21
  [e2ba6199] ExprTools v0.1.11
  [c87230d0] FFMPEG v0.4.6
  [7034ab61] FastBroadcast v1.4.0
  [1a297f60] FillArrays v1.17.1
⌅ [53c48c17] FixedPointNumbers v0.8.6
  [f6369f11] ForwardDiff v1.4.6
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v1.13.0
  [28b8d3ca] GR v0.73.27
  [86223c79] Graphs v1.15.0
  [3263718b] ImplicitDiscreteSolve v2.3.0
  [8197267c] IntervalSets v0.7.14
  [1019f520] JLFzf v0.1.11
  [682c06a0] JSON v1.10.0
  [ccbc3e58] JumpProcesses v9.33.1
  [b964fa9f] LaTeXStrings v1.4.1
  [23fbe1c1] Latexify v0.16.12
  [1914dd2f] MacroTools v0.5.16
  [442fdcdd] Measures v0.3.3
  [7771a370] ModelingToolkitBase v1.77.2
⌅ [2e0e35c7] Moshi v0.3.9
  [46d2c3a1] MuladdMacro v0.2.7
  [77ba4419] NaNMath v1.1.4
  [be0214bd] NonlinearSolveBase v2.54.1
  [6fe1bfb0] OffsetArrays v1.17.0
  [bac558e1] OrderedCollections v2.0.1
  [bbf590c4] OrdinaryDiffEqCore v4.18.1
  [b1df2697] OrdinaryDiffEqTsit5 v2.1.5
  [ccf2f8ad] PlotThemes v3.3.0
  [995b91a9] PlotUtils v1.5.0
  [91a5bcdd] Plots v1.41.7
  [f517fe37] Polyester v0.7.19
  [d236fae5] PreallocationTools v1.7.1
  [aea7be01] PrecompileTools v1.3.4
  [21216c6a] Preferences v1.6.0
  [35bcea6d] QuantumCumulants v0.7.2 `~/work/QuantumCumulants.jl/QuantumCumulants.jl`
  [795d4caa] ReadOnlyDicts v1.0.1
  [3cdcf5f2] RecipesBase v1.4.0
  [01d81517] RecipesPipeline v0.6.12
  [731186ca] RecursiveArrayTools v4.5.3
  [189a3867] Reexport v1.2.2
  [05181044] RelocatableFolders v1.0.1
  [ae029012] Requires v1.3.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.27
  [9dfe8606] SCCNonlinearSolve v1.15.4
  [0bca4576] SciMLBase v3.57.0
  [431bcebd] SciMLPublic v1.3.0
  [53ae85a6] SciMLStructures v1.10.5
  [6c6a2e73] Scratch v1.3.0
  [f7aa4685] SecondQuantizedAlgebra v0.12.0
  [efcf1570] Setfield v1.1.2
  [992d4aef] Showoff v1.1.1
  [727e6d20] SimpleNonlinearSolve v2.14.6
  [276daf66] SpecialFunctions v2.9.0
  [90137ffa] StaticArrays v1.9.22
  [1e83bf80] StaticArraysCore v1.4.4
  [10745b16] Statistics v1.11.5
  [2913bbd2] StatsBase v0.34.13
  [2efcf032] SymbolicIndexingInterface v0.3.55
  [d1185830] SymbolicUtils v4.48.0
  [0c5d862f] Symbolics v7.41.1
  [ed4db957] TaskLocalValues v0.1.3
  [8ea1fca8] TermInterface v2.0.0
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [1cfade01] UnicodeFun v0.4.1
  [41fe7b60] Unzip v0.2.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [f43a241f] Downloads v1.7.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.13.0
  [44cfe95a] Pkg v1.13.0
  [de0858da] Printf v1.11.0
  [3fa0cd96] REPL v1.11.0
  [9a3f8284] Random v1.11.0
  [9e88b42a] Serialization v1.11.0
  [2f01184e] SparseArrays v1.13.0
  [fa267f1f] TOML v1.0.3
  [cf7118a7] UUIDs v1.11.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`

This page was generated using Literate.jl.