Mollow Triplet from a Pulsed Drive

The Mollow Triplet example computes the resonance fluorescence spectrum of a constantly driven atom. Here we drive the atom with a time-dependent field that is smoothly switched on, and compute the emission spectrum from the two-time correlation function once the atom has settled into a quasi-steady state. This showcases correlation functions of time-dependent Hamiltonians.

The atom obeys

\[H(t) = \Delta\,\sigma^{ee} + r(t)\,\Omega\,(\sigma^{ge} + \sigma^{eg}),\]

where the dimensionless envelope $r(t)$ ramps the Rabi drive from $0$ to $1$. The emission spectrum is the Fourier transform of the first-order correlation

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

By the quantum regression theorem, with a time-dependent generator the $\tau$-evolution is governed by $H(t_0+\tau)$, not $H(\tau)$: the drive keeps running during the correlation delay, starting from the absolute time $t_0$ at which the original evolution stopped. QuantumCumulants.jl handles this by substituting the original time variable $t \to t_0 + \tau$ in the correlation equations. We pass $t_0$ via the iv0 keyword and supply its value alongside the other parameters.

using QuantumCumulants
using ModelingToolkitBase, OrdinaryDiffEqTsit5, OrdinaryDiffEqLowOrderRK
using QuantumOptics: timecorrelations
using Plots

We register the envelope $r(t)$ as a time-dependent function. As in the Ramsey Spectroscopy example, we seed a meanfield call first so the envelope is built on the system's own independent variable t.

h = NLevelSpace(:atom, (:g, :e))
σ(i, j) = Transition(h, :σ, i, j)

@variables Δ Ω γ
@variables t₀::Real   # absolute time at which the drive evolution stops (t₀)
@register_symbolic r(t)

eqs_seed = meanfield([σ(:e, :e), σ(:e, :g)], Δ * σ(:e, :e), [σ(:g, :e)]; rates = [γ])
t = eqs_seed.iv

H = Δ * σ(:e, :e) + r(t) * Ω * (σ(:g, :e) + σ(:e, :g))
J = [σ(:g, :e)]

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

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

We drive on resonance ($\Delta = 0$) with $\Omega = 3\gamma$, switching the field on smoothly with $r(t) = 1 - e^{-t/t_r}$ so it has effectively saturated long before the stop time $t_0$.

r(t) = 1 - exp(-t / 2.0)

γv, Ωv, Δv = 1.0, 3.0, 0.0
ps = [γ, Ω, Δ]
p0 = [γv, Ωv, Δv]

sys = mtkcompile(System(eqs; name = :atom))
u0 = initial_values(eqs, zeros(ComplexF64, length(eqs)))
tstop = 20.0
prob = ODEProblem(sys, merge(u0, Dict(ps .=> p0)), (0.0, tstop))
sol = solve(prob, Tsit5(); abstol = 1.0e-10, reltol = 1.0e-10)

plot(
    sol.t,
    real.(get_solution(sol, σ(:e, :e), eqs).(sol.t)),
    xlabel = "γt",
    ylabel = "⟨σᵉᵉ⟩",
    label = "excited-state population",
    size = (600, 300),
)

The population settles onto a plateau: at t₀ = tstop the atom is in a quasi-steady state under the (now constant) drive. We build the emission correlation function on this system, passing iv0 = t₀.

c = CorrelationFunction(σ(:e, :g), σ(:g, :e), eqs; iv0 = t₀)

Without iv0, a time-dependent system raises an informative error, since the correlation equations would otherwise contain the orphaned time variable t:

try
    CorrelationFunction(σ(:e, :g), σ(:g, :e), eqs)
catch e
    println(e isa ArgumentError ? "ArgumentError: iv0 is required" : e)
end

ArgumentError: iv0 is required

We solve the $\tau$-evolution. The steady-state initial values are read off the original solution with correlation_u0; the parameters, including t₀, are propagated with correlation_p0.

csys = mtkcompile(System(c; name = :corr))
u0_c = correlation_u0(c, sol.u[end])
p0_c = correlation_p0(c, sol.u[end], [γ => γv, Ω => Ωv, Δ => Δv, t₀ => tstop])

τ_end = 30.0
prob_c = ODEProblem(csys, merge(u0_c, Dict(p0_c)), (0.0, τ_end))
sol_c = solve(prob_c, Tsit5(); abstol = 1.0e-10, reltol = 1.0e-10, save_idxs = 1)

The spectrum is the Fourier transform of the inelastic (connected) part of the correlation, $g(t_0,\tau) - |\langle\sigma^{ge}\rangle|^2$, with the coherent (elastically scattered) contribution removed. We borrow the FFT helper from QuantumOptics.jl.

coh = abs2(get_solution(sol, σ(:g, :e), eqs)(tstop))   # |⟨σᵍᵉ⟩|² at t₀ (elastic peak)
τ = collect(range(0.0, τ_end; length = 2001))          # equidistant grid for the FFT
g_inel = sol_c.(τ) .- coh
ω, S = timecorrelations.correlation2spectrum(τ, g_inel)

To confirm the result, we compare against the textbook constant-drive Mollow triplet: the same atom driven by a time-independent field $\Omega$, evolved to steady state, with its correlation computed the standard way (no iv0).

H_const = Δ * σ(:e, :e) + Ω * (σ(:g, :e) + σ(:e, :g))
eqs_const = meanfield([σ(:e, :e), σ(:e, :g)], H_const, J; rates = [γ])
complete!(eqs_const)

sys_const = mtkcompile(System(eqs_const; name = :atom_const))
sol_const = solve(
    ODEProblem(sys_const, merge(initial_values(eqs_const, zeros(ComplexF64, 2)), Dict(ps .=> p0)), (0.0, 40.0)),
    Tsit5(); abstol = 1.0e-10, reltol = 1.0e-10,
)

c_const = CorrelationFunction(σ(:e, :g), σ(:g, :e), eqs_const)
csys_const = mtkcompile(System(c_const; name = :corr_const))
sol_cc = solve(
    ODEProblem(
        csys_const,
        merge(
            correlation_u0(c_const, sol_const.u[end]),
            Dict(correlation_p0(c_const, sol_const.u[end], [γ => γv, Ω => Ωv, Δ => Δv]))
        ),
        (0.0, τ_end),
    ),
    Tsit5(); abstol = 1.0e-10, reltol = 1.0e-10, save_idxs = 1,
)
coh_const = abs2(get_solution(sol_const, σ(:g, :e), eqs_const)(40.0))
_, S_const = timecorrelations.correlation2spectrum(τ, sol_cc.(τ) .- coh_const)

The two spectra lie on top of each other: the pulsed drive, sampled from a plateau time $t_0$ via iv0, reproduces the steady-state Mollow triplet. The sidebands sit at $\omega \approx \pm 2\Omega$ (the dressed-state splitting at resonance), flanking the central line at $\omega = 0$.

plot(ω, S; xlims = (-12, 12), xlabel = "ω", ylabel = "S(ω)", label = "pulsed drive (iv0 = t₀)", lw = 2, size = (600, 350))
plot!(ω, S_const; xlims = (-12, 12), label = "constant drive (steady state)", ls = :dash, lw = 2)

Package versions

These results were obtained using the following versions:

using InteractiveUtils
versioninfo()

using Pkg
Pkg.status(
    ["QuantumCumulants", "OrdinaryDiffEqTsit5", "ModelingToolkitBase", "QuantumOptics", "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
⌅ [7d9fca2a] Arpack v0.5.3
  [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
  [77a26b50] DiffEqNoiseProcess v5.36.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
  [7a1cc6ca] FFTW v1.10.0
  [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
  [42fd0dbc] IterativeSolvers v0.9.4
  [1019f520] JLFzf v0.1.11
  [682c06a0] JSON v1.10.0
  [ccbc3e58] JumpProcesses v9.33.1
  [0b1a1467] KrylovKit v0.10.4
  [b964fa9f] LaTeXStrings v1.4.1
  [23fbe1c1] Latexify v0.16.12
  [7a12625a] LinearMaps v3.11.4
  [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
  [1344f307] OrdinaryDiffEqLowOrderRK v2.2.6
  [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`
  [6e0679c1] QuantumOptics v1.2.10
  [4f57444f] QuantumOpticsBase v0.5.16
  [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
  [789caeaf] StochasticDiffEq v7.2.0
  [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.