Laser with Filter Cavities

An intuitive and straightforward approach to calculate the spectrum of a laser is to filter the emitted light. We can do this by coupling filter cavities with different detunings to the main cavity and observe the photon number in these 'filters', see for example K. Debnath et al., Phys Rev A 98, 063837 (2018).

The main goal of this example is to combine two indexed Hilbert spaces, where one will be scaled and the other evaluated. The model is basically the same as for the superradiant laser example, but with the additional terms due to the filter cavities. The Hamiltonian of this system is

\[\begin{equation} H = - \Delta a^\dagger a + g \sum\limits_{j=1}^{N} (a^\dagger \sigma^{12}_{j} + a \sigma^{21}_{j}) - \sum\limits_{i=1}^{M} \delta_i b_i^\dagger b_i + g_f \sum\limits_{i=1}^{M} (a^\dagger b_i + a b_i^\dagger), \end{equation}\]

where $\delta_i$ is the detuning of the $i$-th filter cavity and $g_f$ the coupling to the normal cavity. Furthermore, their decay rate is $\kappa_f$.

We start by loading the packages.

using QuantumCumulants
using OrdinaryDiffEqTsit5, ModelingToolkitBase
using Plots

We create the parameters of the system including the $\texttt{IndexedVariable}$ $\delta_i$. For the atoms and filter cavities we only need one Hilbert space each. We define the indices for each Hilbert space and use them to create $\texttt{IndexedOperators}$.

@variables κ g gf κf R Γ Δ ν N M # Parameters
δ(i) = IndexedVariable(:δ, i)

hc = FockSpace(:cavity) # Hilbert spaces
hf = FockSpace(:filter)
ha = NLevelSpace(:atom, 2)
h = hc ⊗ hf ⊗ ha

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

@qnumbers a::Destroy(h, 1)
b(k) = IndexedOperator(Destroy(h, :b, 2), k)

i is bound by the Hamiltonian sums and the dissipator. The canonical free-index slot on the filter Hilbert space mints i_2 (lex-first declared index, suffix 2) to keep state names disjoint from H's bound scope. Per-atom state lookups thus use i_2(k) = i_2_k, not i(k).

b(k::Integer) = IndexedOperator(Destroy(h, :b, 2), i(2)(k))
σ(α, β, k) = IndexedOperator(Transition(h, :σ, α, β, 3), k)
σ(α, β, k::Integer) = IndexedOperator(Transition(h, :σ, α, β, 3), j(k))

We define the Hamiltonian using symbolic sums and define the individual dissipative processes. For an indexed jump operator the (symbolic) sum is build in the Liouvillian, in this case corresponding to individual decay processes.

H =
    Δ * Σ(σ(2, 2, j), j) +
    Σ(δ(i) * b(i)'b(i), i) +
    gf * (Σ(a' * b(i) + a * b(i)', i)) +
    g * (Σ(a' * σ(1, 2, j) + a * σ(2, 1, j), j)) # Hamiltonian

J = [a, b(i), σ(1, 2, j), σ(2, 1, j), σ(2, 2, j)] # Jumps & rates
rates = [κ, κf, Γ, R, ν]

We derive the equation for $\langle a^\dagger a \rangle$ and complete the system automatically in second order.

eqs = meanfield(a'a, H, J; rates = rates, order = 2)

eqs_c = complete!(deepcopy(eqs));

Now we assume that all atoms behave identically, but we want to obtain the equations for 20 different filter cavities. To this end we $\texttt{scale}$ the Hilbert space of the atoms and $\texttt{evaluate}$ the filter cavities. Specifying the Hilbert space is done with the kwarg $\texttt{h}$, which can either be the specific Hilbert space or its acts-on number. Evaluating the filter cavities requires a numeric upper bound for the used $\texttt{Index}$, we provide this with a dictionary on the kwarg $\texttt{limits}$.

M_ = 20
eqs_sc = scale(eqs_c; h = [3])
eqs_eval = evaluate(eqs_sc; limits = Dict(M => M_))
println("Number of eqs.: $(length(eqs_eval))")

Number of eqs.: 552

To calculate the dynamics of the system we create a system of ordinary differential equations, which can be used by DifferentialEquations.jl. Finally we need to define the numerical parameters and the initial state of the system.

sys = System(eqs_eval; name = :sys)
sys_c = mtkcompile(sys)

u0 = zeros(ComplexF64, length(eqs_eval)) # Initial state


N_ = 200 # System parameters
Γ_ = 1.0
Δ_ = 0Γ_
g_ = 1Γ_
κ_ = 100Γ_
R_ = 10Γ_
ν_ = 1Γ_

gf_ = 0.1Γ_
κf_ = 0.1Γ_
δ_ls = [0:(1 / M_):(1 - 1 / M_);] * 10Γ_

pmap = parameter_map(
    eqs_eval, Dict(
        Γ => Γ_, κ => κ_, g => g_, κf => κf_, gf => gf_, R => R_,
        δ(i) => δ_ls,
        Δ => Δ_, ν => ν_, N => N_,
    )
)
dict = merge(initial_values(eqs_eval, u0), pmap)
prob = ODEProblem(sys_c, dict, (0.0, 10.0 / κf_))

sol = solve(prob, Tsit5(); abstol = 1.0e-10, reltol = 1.0e-10, maxiters = 1.0e7) # Solve the numeric problem

t = sol.t
n = abs.(get_solution(sol, a'a, eqs_eval).(t))
n_b(i) = abs.(get_solution(sol, b(i)'b(i), eqs_eval).(t))
n_f = [abs(get_solution(sol, b(i)'b(i), eqs_eval)(t[end])) for i in 1:M_] ./
    (abs(get_solution(sol, b(1)'b(1), eqs_eval)(t[end])))

p1 = plot(t, n_b(1), alpha = 0.5, ylabel = "⟨bᵢ⁺bᵢ⟩", legend = false) # Plot results
for i in 2:M_
    plot!(t, n_b(i), alpha = 0.5, legend = false)
end

p2 = plot(
    [-reverse(δ_ls); δ_ls],
    [reverse(n_f); n_f],
    xlabel = "δ/Γ",
    ylabel = "intensity",
    legend = false,
)
plot(p1, p2, layout = (1, 2), size = (700, 300))

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.0
Commit d1c37793dd2 (2026-09-09 19:00 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 7763 64-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver3)
  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.12.0
  [caf10ac8] BipartiteGraphs v0.1.14
  [8e7c35d0] BlockArrays v1.10.0
⌅ [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.1
  [459566f4] DiffEqCallbacks v4.19.4
  [b552c78f] DiffRules v1.16.0
  [a0c0ee7d] DifferentiationInterface v0.7.21
  [ffbed154] DocStringExtensions v0.9.5
  [5b8099bc] DomainSets v0.8.1
  [4e289a0a] EnumX v1.0.7
  [f151be2c] EnzymeCore v0.8.21
  [e2ba6199] ExprTools v0.1.11
  [c87230d0] FFMPEG v0.4.5
  [7034ab61] FastBroadcast v1.4.0
  [1a297f60] FillArrays v1.17.0
⌅ [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.8.0
  [ccbc3e58] JumpProcesses v9.32.3
  [b964fa9f] LaTeXStrings v1.4.1
  [23fbe1c1] Latexify v0.16.12
  [1914dd2f] MacroTools v0.5.16
  [442fdcdd] Measures v0.3.3
  [7771a370] ModelingToolkitBase v1.71.2
⌅ [2e0e35c7] Moshi v0.3.9
  [46d2c3a1] MuladdMacro v0.2.7
  [77ba4419] NaNMath v1.1.4
  [be0214bd] NonlinearSolveBase v2.49.5
  [6fe1bfb0] OffsetArrays v1.17.0
  [bac558e1] OrderedCollections v2.0.1
  [bbf590c4] OrdinaryDiffEqCore v4.17.2
  [b1df2697] OrdinaryDiffEqTsit5 v2.1.4
  [ccf2f8ad] PlotThemes v3.3.0
  [995b91a9] PlotUtils v1.4.4
  [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.3.4
  [01d81517] RecipesPipeline v0.6.12
  [731186ca] RecursiveArrayTools v4.5.1
  [189a3867] Reexport v1.2.2
  [05181044] RelocatableFolders v1.0.1
  [ae029012] Requires v1.3.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.26
  [9dfe8606] SCCNonlinearSolve v1.15.3
  [0bca4576] SciMLBase v3.54.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.5
  [276daf66] SpecialFunctions v2.9.0
  [90137ffa] StaticArrays v1.9.20
  [1e83bf80] StaticArraysCore v1.4.4
  [10745b16] Statistics v1.11.5
  [2913bbd2] StatsBase v0.34.13
  [2efcf032] SymbolicIndexingInterface v0.3.55
  [d1185830] SymbolicUtils v4.46.6
  [0c5d862f] Symbolics v7.39.2
  [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.