Mollow Triplet

In the following example, we show how one can compute the resonance spectrum of a single, coherently driven atom. The Hamiltonian describing the system is given by

\[H = -\Delta\sigma^{ee} + \Omega\left(\sigma^{ge} + \sigma^{eg}\right),\]

where $\Delta = \omega_\ell - \omega_a$ is the detuning between the laser and the atomic resonance frequency and $\Omega$ is the Rabi frequency of the laser. The atom is also subject to decay at rate $\gamma$ with the damping operator $\sigma^{ge} = |g\rangle \langle e |$ projecting it from the excited state $|e\rangle $ to its ground state $|g\rangle$.

using QuantumCumulants
using ModelingToolkitBase, OrdinaryDiffEqLowOrderRK
using Plots

Obtaining the equations of motion for the system is simple. Note that in this case, we are not actually making any assumptions: in the special case of a single atom (even if it has more than two levels), there exists a simple mapping between the equations of motions for averages and the density operator. In our case, the density operator $\rho$ is determined by only two average values, $\langle\sigma^{eg}\rangle = \rho_{eg} = \rho_{ge}^* $ and $\langle \sigma^{ee}\rangle = \rho_{ee} = 1 - \rho_{gg}$. In other words, we are solving the master equation component-wise.

h = NLevelSpace(:atom, (:g, :e)) # Hilbert space

@variables Δ Ω γ # Operators
σ(i, j) = Transition(h, :σ, i, j)
H = Δ * σ(:e, :e) + Ω * (σ(:g, :e) + σ(:e, :g))
J = [σ(:g, :e)]


eqs = meanfield([σ(:e, :g), σ(:e, :e)], H, J; rates = [γ]) # Equations
complete!(eqs)

\[ \begin{aligned} \partial_{t} \langle {\sigma}^{{21}} \rangle &= i \Omega + \langle {\sigma}^{{21}} \rangle \left( - \frac{1}{2} \gamma + i \Delta \right) - 2 \langle {\sigma}^{{22}} \rangle i \Omega \\\\[-0.0em] \partial_{t} \langle {\sigma}^{{22}} \rangle &= - \langle {\sigma}^{{22}} \rangle \gamma + \langle {\sigma}^{{12}} \rangle i \Omega - \langle {\sigma}^{{21}} \rangle i \Omega \end{aligned} \]

In order to compute the correlation function given by

\[g(\tau) = \langle \sigma^{eg}(t_0+\tau)\sigma^{ge}(t_0)\rangle \equiv \langle \sigma^{eg}\sigma^{ge}_0\rangle,\]

where $t_0$ is any time after which the original system has reached steady state. Note, that in the second step we simplified the notation by dropping the temporal arguments. The subscript $0$ indicates operators given at $t_0$.

The correlation function of the system given by eqs can be computed as follows.

c = CorrelationFunction(σ(:e, :g), σ(:g, :e), eqs; steady_state = true) # Correlation Function

Note that the above actually derives another set of equations, which might take some time. Specifically, it is deriving the equation of motion for $g(\tau)$ and automatically completes the system of equation deriving the necessary equations of motion for all average values on which $g(\tau)$ depends. Here, the system of equations reads

c.eqs

\[ \begin{aligned} \partial_{\tau} \langle {\sigma}^{{21}}{\sigma_{\mathrm{0}}}^{{12}} \rangle &= \langle {\sigma}^{{21}}{\sigma_{\mathrm{0}}}^{{12}} \rangle \left( - \frac{1}{2} \gamma + i \Delta \right) + \langle {\sigma_{\mathrm{0}}}^{{12}} \rangle i \Omega - 2 \langle {\sigma}^{{22}}{\sigma_{\mathrm{0}}}^{{12}} \rangle i \Omega \\\\[-0.0em] \partial_{\tau} \langle {\sigma}^{{22}}{\sigma_{\mathrm{0}}}^{{12}} \rangle &= - \langle {\sigma}^{{22}}{\sigma_{\mathrm{0}}}^{{12}} \rangle \gamma + \langle {\sigma}^{{12}}{\sigma_{\mathrm{0}}}^{{12}} \rangle i \Omega - \langle {\sigma}^{{21}}{\sigma_{\mathrm{0}}}^{{12}} \rangle i \Omega \\\\[-0.0em] \partial_{\tau} \langle {\sigma}^{{12}}{\sigma_{\mathrm{0}}}^{{12}} \rangle &= \langle {\sigma}^{{12}}{\sigma_{\mathrm{0}}}^{{12}} \rangle \left( - \frac{1}{2} \gamma - i \Delta \right) - \langle {\sigma_{\mathrm{0}}}^{{12}} \rangle i \Omega + 2 \langle {\sigma}^{{22}}{\sigma_{\mathrm{0}}}^{{12}} \rangle i \Omega \end{aligned} \]

As mentioned above, to compute the time evolution of the system itself, we are effectively solving a master equation. One way to obtain the spectrum would be to solve the system in order to obtain the time evolution of the correlation function, and then take the Fourier transform. However, a computationally more efficient way is to take the Laplace transform of the correlation function directly, which yields a simple linear system of equations. Then we only have to compute a matrix inverse instead of the time evolution to obtain the spectrum.

In the following, we will use the latter approach:

ps = (Δ, Ω, γ)
S = Spectrum(c, ps)

\[\mathcal{F}(\langle {\sigma}^{{21}}{\sigma_{\mathrm{0}}}^{{12}} \rangle)(\omega)\]

Constructing the Spectrum from the correlation function, automatically derives a matrix $A$ and a vector $b$ determining the set of linear equations

\[A(\omega)x(\omega) = b(\omega),\]

where the spectrum is then given by the first entry of the solution vector, $S(\omega)=x_1(\omega)$.

To find the spectrum, we first need to compute the time evolution of the system up to steady state.

sys = System(eqs; name = :sys)
ssys = mtkcompile(sys)

p0 = (0.0, 2.0, 1.0)
u0 = zeros(ComplexF64, 2)
dict = merge(Dict(unknowns(ssys) .=> u0), Dict(ps .=> p0))
prob = ODEProblem(ssys, dict, (0.0, 20.0))
sol = solve(prob, RK4())

plot(
    sol.t,
    real.(get_solution(sol, σ(:e, :e), eqs).(sol.t)),
    xlabel = "γt",
    label = "Excited state population",
)

Now, solving the linear system to obtain the spectrum can simply be done by calling the instance at a range of frequencies, and providing the proper steady-state values and numerical parameters.

ω = range(-6pi, 6pi, length = 1001)
s = S(ω, sol.u[end], p0)
plot(ω, s, xlabel = "ω - ωℓ", label = "Resonance spectrum")

The resulting spectrum shows a prominent peak at the resonance point ($\omega=\omega_\ell=\omega_a$), but also two more peaks around $\omega \approx \pm \Omega^2/\gamma$. These two resonances originate from the dressed states. These three peaks are called Mollow Triplet.

Package versions

These results were obtained using the following versions:

using InteractiveUtils
versioninfo()

using Pkg
Pkg.status(
    ["QuantumCumulants", "OrdinaryDiffEqLowOrderRK", "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 × INTEL(R) XEON(R) PLATINUM 8573C
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, emeraldrapids)
  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.1
  [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.3
  [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.70.0
⌅ [2e0e35c7] Moshi v0.3.9
  [46d2c3a1] MuladdMacro v0.2.7
  [77ba4419] NaNMath v1.1.4
  [6fe1bfb0] OffsetArrays v1.17.0
  [bac558e1] OrderedCollections v2.0.1
  [bbf590c4] OrdinaryDiffEqCore v4.17.1
  [1344f307] OrdinaryDiffEqLowOrderRK v2.2.5
  [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
  [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.53.3
  [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.3
  [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.5
  [0c5d862f] Symbolics v7.39.2
  [ed4db957] TaskLocalValues v0.1.3
  [8ea1fca8] TermInterface v2.0.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.