Many-atom laser

This example describes a second order laser system consisting of $N$ three-level atoms coupled to a single mode cavity. An auxiliary state $|3\rangle$, which quickly decays into the upper lasing state $|2\rangle$, is coherently pumped to achieve population inversion on the lasing transition $|1\rangle \leftrightarrow |2\rangle$. The Hamiltonian of this system is given by

\[H = -\Delta_{c} a^{\dagger} a - \sum_{i=1}^N \left[ \Delta_3^i \sigma_i^{33} + g_i (a^{\dagger} \sigma_i^{12} + a\sigma_i^{21}) + \Omega_i (\sigma_i^{31} + \sigma_i^{13}) \right].\]

Including dissipative processes as, e.g. the atomic decay or photon losses through the cavity mirrors, makes it an open quantum system. In the Schrödinger picture we would compute the dynamics of such open quantum systems with a density matrix $\rho$ according to a master equation (see e.g. https://docs.qojulia.org/),

\[\frac{d}{dt} \rho = - \frac{i}{\hbar} \left[ H, \rho \right] + \mathcal{L}[\rho],\]

with $\mathcal{L}[\rho] = \frac{\gamma}{2} (2 J \rho J^\dagger - J^\dagger J \rho - \rho J^\dagger J)$ the Liouvillian superoperator in standard Lindblad form for a dissipative process with jump operator $J$ and rate $R$.

With QuantumCumulants.jl we describe the system dynamics with averages, which are deduced from the operator equations of motion in the Heisenberg picture. In the Heisenberg picture open systems are described by the quantum Langevin equation. Assuming white noise, we can omit the stochastic terms of the quantum Langevin equation when computing averages. Thus we get the following equation for the time evolution of a system operator average $\langle O \rangle$ (if $O$ is not explicitly time dependent):

\[\frac{d}{dt} \langle O \rangle = \frac{i}{\hbar} \left[ H, O \right] + \bar{\mathcal{L}}[O].\]

The superoperator $\bar{\mathcal{L}}[O]$ is similar to the Lindblad term in the Schrödinger picture, except that $J$ and $J^\dagger$ are swapped in the first term, i.e. $\bar{\mathcal{L}}[O] = \frac{\gamma}{2} (2 J^\dagger O J - J^\dagger J O - O J^\dagger J)$, for a dissipative process with jump operator $J$ and rate $R$.

For our system we have four different dissipative processes with the jump operators $a$, $\sigma^{12}_i$, $\sigma^{13}_i$ and $\sigma^{23}_i$, and the corresponding decay rates $\kappa$, $\Gamma^i_{12}$, $\Gamma^i_{13}$ and $\Gamma^i_{23}$, respectively.

We start by loading the needed packages.

using QuantumCumulants
using ModelingToolkitBase, OrdinaryDiffEqTsit5
using Plots

Then we define the symbolic parameters of the system, the Hilbertspace and the necessary operators. We define an atomic transition operator function $\sigma(i,j,k)$ for the transition from $|j \rangle$ to $|i \rangle$ of atom $k$. Since we only have one FockSpace we do not need to specify the Hilbertspace on which the Destroy operator acts. For the different atomic transitions, however, we need to specify this, since there is more than one NLevelSpace. This information is stored in the .aon field of each operator.

N = 2 # number of atoms
@variables κ g Γ₂₃ Γ₁₃ Γ₁₂ Ω Δc Δ₃

hf = FockSpace(:cavity) # Hilbertspace
ha = ⊗([NLevelSpace(Symbol(:atom, i), 3) for i in 1:N]...)
h = hf ⊗ ha

a = Destroy(h, :a) # Operators
σ(i, j, k) = Transition(h, Symbol("σ_{$k}"), i, j, k + 1)

Now we create the Hamiltonian and the jumps with the corresponding rates of our laser system. We assume here that all atoms are identical.

H =
    -Δc * a'a +
    sum(g * (a' * σ(1, 2, i) + a * σ(2, 1, i)) for i in 1:N) +
    sum(Ω * (σ(3, 1, i) + σ(1, 3, i)) for i in 1:N) - sum(Δ₃ * σ(3, 3, i) for i in 1:N) # Hamiltonian

J = [a; [σ(1, 2, i) for i in 1:N]; [σ(1, 3, i) for i in 1:N]; [σ(2, 3, i) for i in 1:N]] # Jumps

rates = [κ; [Γ₁₂ for i in 1:N]; [Γ₁₃ for i in 1:N]; [Γ₂₃ for i in 1:N]] # Rates

Later we will complete the system automatically, which has the disadvantage that the equations are not ordered. Therefore we define a list of interesting operators, which we want to use later. Note that at least one operator(-product) is needed. We derive the equations for these operators, average them, and automatically complete the system of equations.

ops = [a'a, σ(2, 2, 1), σ(3, 3, 1)] # list of operators

eqs = meanfield(ops, H, J; rates = rates, order = 2) #second order average

complete!(eqs) # automatically complete the system

To calculate the time evolution we create a Julia function which can be used by DifferentialEquations.jl to solve the set of ordinary differential equations.

Build a System out of the MeanfieldEquations

sys = System(eqs; name = :laser)
sys_c = mtkcompile(sys)

Finally, we compute the time evolution after defining an initial state and numerical values for the parameters.

u0 = initial_values(eqs) # initial state

Γ₁₂n = 1.0
Γ₂₃n = 20Γ₁₂n
Γ₁₃n = 2Γ₁₂n
Ωn = 5Γ₁₃n
gn = 2Γ₁₂n
Δcn = 0.0
Δ₃n = 0.0
κn = 0.5Γ₁₂n

ps = (g, Γ₂₃, Γ₁₃, Γ₁₂, Ω, Δc, Δ₃, κ) # list of parameters
p0 = Dict(ps .=> (gn, Γ₂₃n, Γ₁₃n, Γ₁₂n, Ωn, Δcn, Δ₃n, κn))
tend = 10.0 / κn

prob = ODEProblem(sys_c, merge(u0, p0), (0.0, tend))
sol = solve(prob, Tsit5(), reltol = 1.0e-8, abstol = 1.0e-8)

We plot the average photon number and the population inversion of the lasing transition.

n_t = real.(get_solution(sol, a' * a, eqs).(sol.t))
σ22_t = real.(get_solution(sol, σ(2, 2, 1), eqs).(sol.t))
σ33_t = real.(get_solution(sol, σ(3, 3, 1), eqs).(sol.t))
σ22m11_t = 2 .* σ22_t .+ σ33_t .- 1 #σ11 + σ22 + σ33 = 𝟙

p1 = plot(sol.t, n_t, xlabel = "tΓ₁₂", ylabel = "⟨a⁺a⟩", legend = false) # Plot
p2 = plot(sol.t, σ22m11_t, xlabel = "tΓ₁₂", ylabel = "⟨σ22⟩ - ⟨σ11⟩", legend = false)
plot(p1, p2, layout = (1, 2), size = (800, 320), left_margin = 5Plots.mm, bottom_margin = 5Plots.mm, dpi = 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.