Heterodyne detection of emission from atomic ensemble
This example implements the heterodyne detection of the light emitted by an ensemble of atoms, described in H. Yu, et. al., Phys. Rev. Lett 133, 073601 (2024). We describe the stochastic master equation time evolution that governs the measurement backaction on the system.
using QuantumCumulants
using ModelingToolkitBase
using OrdinaryDiffEqLowOrderRK
using StochasticDiffEqLowOrder: SDEProblem, EM, EnsembleProblem, ReturnCode
using DiffEqNoiseProcess: RealWienerProcess
using PlotsThe system we discuss is an ensemble of $N$ equivalent two level systems with transition operators $\sigma^{\alpha \beta}_j$ and frequency $\omega_a$) coupled to a cavity mode $\hat a$, all with the same coupling $g$ and free space decay rate $\gamma$. The cavity has frequency $\omega_c$ and decay rate $\kappa$. We are in the frame rotating with the frequency $\omega_a=\omega_c$.
The system Hamiltonian is
\[H=\omega_c\hat a^\dagger\hat a+\omega_a\sum_j\sigma^{22}_j+\hat a^\dagger\sigma^{12}_j+g\hat a\sigma^{21}_j.\]
@variables N ωₐ γ η χ ωc κ g ξ ωₗ
@register_symbolic pulse(t)
hc = FockSpace(:resonator)
ha = NLevelSpace(:atom, 2)
h = hc ⊗ ha
j = Index(h, :j, N, ha)
k = Index(h, :k, N, ha)
@qnumbers a::Destroy(h, 1)
σ(α, β, k) = IndexedOperator(Transition(h, :σ, α, β, 2), k)
eqs_seed = meanfield(a, ωc * a' * a, [a]; rates = [κ])
t = eqs_seed.iv
H = ωc * a' * a + ωₐ * Σ(σ(2, 2, j), j) + g * a' * Σ(σ(1, 2, j), j) + g * a * Σ(σ(2, 1, j), j);We include decay of the cavity mode, where we choose the decay operator to be $\hat a\exp(\rm{i}\omega_lt)$ for convenience (we will see later why this is reasonable, and you can convince yourself that it gives the same decay term as $\hat a$). The individual decay of the atoms with rate $\gamma$ is given by the decay operator $\sigma^{12}_j$.
We also include an incoherent pump $\sigma^{21}_j$ with amplitude $\eta$, which is a pulse that is on between $t_0$ and $t_0+t_1$. The last Lindblad term is a dephasing term with strength $\chi$ corresponding to the operator $\sigma^{22}_j$.
J = [a * exp(1.0im * ωₗ * t), σ(1, 2, j), σ(2, 1, j), σ(2, 2, j)]
rates = [κ, γ, η * pulse(t), 2 * χ]
pulse(t) = (t > t₀ && t < t₀ + t₁) * 1.0The measurement terms are included by defining measurement efficiencies for all decay channels. A channel with efficiency zero is ignored, i.e. not measured. The operator $\hat a\exp(\rm{i}\omega_lt)$ corresponding to heterodyne detection with local oscillator frequency $\omega_l$ is measured with efficiency $\xi$.
Such measurement terms are then included in the equation of motion for system as
\[\begin{align} d\langle\hat A\rangle=\sqrt{\xi\kappa/2}\langle \hat a^\dagger \rm{e}^{-\rm{i}\omega_lt}\hat A+\hat A\hat a\rm{e}^{\rm{i}\omega_lt} -\hat a\langle\hat a \rm{e}^{\rm{i}\omega_lt}+\hat a^\dagger \rm{e}^{-\rm{i}\omega_lt}\rangle\rangle \end{align}\]
for the operator $\hat a \rm{e}^{\rm{i}\omega_lt}$ corresponding to heterodyne detection.
efficiencies = [ξ, 0, 0, 0]
ops = [a', a' * a, σ(2, 2, k), σ(1, 2, k), a * a]
eqs = meanfield(ops, H, J; rates = rates, efficiencies = efficiencies, direction = Forward(), order = 2, iv = t)\[ \begin{aligned} \partial_{t} \langle a^{\dagger} \rangle &= g i \underset{j}{\overset{N}{\sum}}\langle {\sigma}_{j}^{{21}} \rangle + \langle a^{\dagger} \rangle \left( i \mathtt{{\omega}c} - \frac{1}{2} \cos^{2}\left( t \mathtt{\omega_l} \right) \kappa - \frac{1}{2} \sin^{2}\left( t \mathtt{\omega_l} \right) \kappa \right) + \frac{\mathrm{d}W}{\mathrm{d}t} \left( \langle a^{\dagger}a \rangle \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} + \langle a^{\dagger}a^{\dagger} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} + \langle a^{\dagger} \rangle \left( - \langle a \rangle \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) - \langle a^{\dagger} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \right) \sqrt{\kappa \xi} \right) \\[-0.0em] \partial_{t} \langle a^{\dagger}a \rangle &= g i \underset{j}{\overset{N}{\sum}}\langle a{\sigma}_{j}^{{21}} \rangle - g i \underset{j}{\overset{N}{\sum}}\langle a^{\dagger}{\sigma}_{j}^{{12}} \rangle + \frac{1}{2} \langle a^{\dagger}a \rangle \left( - 2 \cos^{2}\left( t \mathtt{\omega_l} \right) - 2 \sin^{2}\left( t \mathtt{\omega_l} \right) \right) \kappa + \frac{\mathrm{d}W}{\mathrm{d}t} \left( \langle a^{\dagger}a \rangle \left( - \langle a \rangle \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) - \langle a^{\dagger} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \right) \sqrt{\kappa \xi} + \left( \langle a^{\dagger} \rangle \langle aa \rangle + 2 \langle a \rangle \langle a^{\dagger}a \rangle - 2 \langle a \rangle^{2} \langle a^{\dagger} \rangle \right) \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} + \left( \langle a^{\dagger}a^{\dagger} \rangle \langle a \rangle + 2 \langle a^{\dagger} \rangle \langle a^{\dagger}a \rangle - 2 \langle a^{\dagger} \rangle^{2} \langle a \rangle \right) \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} \right) \\[-0.0em] \partial_{t} \langle {\sigma}_{k}^{{22}} \rangle &= \mathrm{pulse}\left( t \right) \eta + \langle a^{\dagger}{\sigma}_{k}^{{12}} \rangle g i - \langle a{\sigma}_{k}^{{21}} \rangle g i + \langle {\sigma}_{k}^{{22}} \rangle \left( - \gamma - \mathrm{pulse}\left( t \right) \eta \right) + \frac{\mathrm{d}W}{\mathrm{d}t} \left( \langle a{\sigma}_{k}^{{22}} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} + \langle a^{\dagger}{\sigma}_{k}^{{22}} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} + \langle {\sigma}_{k}^{{22}} \rangle \left( - \langle a \rangle \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) - \langle a^{\dagger} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \right) \sqrt{\kappa \xi} \right) \\[-0.0em] \partial_{t} \langle {\sigma}_{k}^{{12}} \rangle &= - \langle a \rangle g i + 2 \langle a{\sigma}_{k}^{{22}} \rangle g i + \langle {\sigma}_{k}^{{12}} \rangle \left( - \frac{1}{2} \gamma - \chi - i \mathtt{\omega_a} - \frac{1}{2} \mathrm{pulse}\left( t \right) \eta \right) + \frac{\mathrm{d}W}{\mathrm{d}t} \left( \langle a{\sigma}_{k}^{{12}} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} + \langle a^{\dagger}{\sigma}_{k}^{{12}} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} + \langle {\sigma}_{k}^{{12}} \rangle \left( - \langle a \rangle \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) - \langle a^{\dagger} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \right) \sqrt{\kappa \xi} \right) \\[-0.0em] \partial_{t} \langle aa \rangle &= - 2 g i \underset{j}{\overset{N}{\sum}}\langle a{\sigma}_{j}^{{12}} \rangle + \langle aa \rangle \left( - 2.0 i \mathtt{{\omega}c} - \cos^{2}\left( t \mathtt{\omega_l} \right) \kappa - \sin^{2}\left( t \mathtt{\omega_l} \right) \kappa \right) + \frac{\mathrm{d}W}{\mathrm{d}t} \left( \langle aa \rangle \left( - \langle a \rangle \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) - \langle a^{\dagger} \rangle \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \right) \sqrt{\kappa \xi} + \left( 3 \langle a \rangle \langle aa \rangle - 2 \langle a \rangle^{3} \right) \left( \cos\left( t \mathtt{\omega_l} \right) + \mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} + \left( \langle a^{\dagger} \rangle \langle aa \rangle + 2 \langle a \rangle \langle a^{\dagger}a \rangle - 2 \langle a \rangle^{2} \langle a^{\dagger} \rangle \right) \left( \cos\left( t \mathtt{\omega_l} \right) -\mathit{i} \sin\left( t \mathtt{\omega_l} \right) \right) \sqrt{\kappa \xi} \right) \end{aligned} \]
The completion and scaling as previously discussed in other examples work exactly the same way for equations including noise terms.
eqs_c = complete(eqs; get_adjoints = false)
scaled_eqs = scale(eqs_c)Here we define the actual values for the system parameters. We then show that the deterministic time evolution without the noise terms is still accessible by using the constructor System for the stochastic system of equations and the syntax for the simulation of the time evolution is as usual.
# frequencies are in kHz
ωc_ = 0.0
κ_ = 2π * 1.13e3
ξ_ = 0.12
N_ = 5.0e4
ωₐ_ = 0.0
γ_ = 2π * 0.375
η_ = 2π * 20
χ_ = 0.016
g_ = 2π * 0.73
ωₗ_ = 2π * 1.0e3
t₀ = 0.0
t₁ = 20.0e-3
p = [N, ωₐ, γ, η, χ, ωc, κ, g, ξ, ωₗ]
p0 = [N_, ωₐ_, γ_, η_, χ_, ωc_, κ_, g_, ξ_, ωₗ_]
T_end = 0.1 # 0.1msTo compare the stochastic evolution to the no-measurement case, build the deterministic system by passing MeanfieldEquations(scaled_eqs) to System. The constructor strips the noise drift, which also drops ξ from the compiled parameter list (it appeared only inside the sqrt(ξ κ) factor of the noise term). Because the physical parameter dict p .=> p0 still carries ξ, we filter it through parameter_map(sys, …); the helper keeps only the entries whose key is a live parameter or unknown of the compiled system and silently drops the rest.
sys = mtkcompile(System(MeanfieldEquations(scaled_eqs); name = :sys))
u0 = zeros(ComplexF64, length(scaled_eqs))
dict = parameter_map(sys, merge(Dict(unknowns(sys) .=> u0), Dict(p .=> p0)))
prob = ODEProblem(sys, dict, (0.0, T_end))
sol_det = solve(prob, RK4(), dt = T_end / 2.0e5)
graph = plot(xlabel = "Time (ms)", ylabel = "Cavity photon number")
plot!(graph, sol_det.t, real(get_solution(sol_det, a'a, scaled_eqs)(sol_det.t)), legend = false)
The stochastic time evolution is accessible via the constructor SDESystem, whose syntax is exactly the same as for the System, but with keyword args as defined in the SDE tutorial. We need to provide a noise process for the measurement. If the noise is white the appropriate noise process is a Wiener process. The SDEProblem is constructed just as the ODEProblem, but with an additional noise argument.
We can then make use of the EnsembleProblem, which automatically runs multiple instances of the stochastic equations of motion. The number of trajectories can be set in the solve call. See the tutorial linked above for more details of the function calls here.
sys_st = mtkcompile(System(scaled_eqs; name = :sys_st))
dict_st = parameter_map(sys_st, merge(Dict(unknowns(sys_st) .=> u0), Dict(p .=> p0)))
noise = RealWienerProcess(0.0, 0.0)
prob_st = SDEProblem(sys_st, dict_st, (0.0, T_end); noise = noise)
sol_test = solve(prob_st, EM(); dt = T_end / 2.0e5);
plot(
sol_test.t,
real(get_solution(sol_test, a'a, scaled_eqs)(sol_test.t)),
xlabel = "Time (ms)",
ylabel = "Cavity photon number",
legend = false,
)
eprob = EnsembleProblem(prob_st)
traj = 200
tspan = range(0.0, T_end, length = 201)
sol = solve(
eprob,
EM(),
dt = T_end / 2.0e5,
save_noise = true,
trajectories = traj,
saveat = tspan,
)We plot the average of the cavity photon number for the stochastic and deterministic equation of motion with the trajectories in gray in the background. We can see that the dynamics of the system is indeed modified by the measurement backaction.
n_avg_ = zeros(length(tspan)) # average photon number
a_real_avg_ = zeros(length(tspan)) # average field (real part)
traj_succ = zeros(traj) # trajectories can be numerically unstable
for i in 1:traj
sol_ = sol.u[i]
if sol_.retcode == ReturnCode.Success
n_avg_ .+= real(get_solution(sol_, a'a, scaled_eqs)(tspan))
a_real_avg_ .+= real(get_solution(sol_, a, scaled_eqs)(tspan))
traj_succ[i] = 1
end
end
@show sum(traj_succ)
n_avg = n_avg_ / sum(traj_succ)
a_real_avg = a_real_avg_ / sum(traj_succ)
graph1 = plot(xlabel = "Time (ms)", ylabel = "Cavity photon number")
for i in 1:traj
plot!(
graph1,
sol.u[i].t,
real(get_solution(sol.u[i], a'a, scaled_eqs)(sol.u[i].t)),
color = :grey,
alpha = 0.75,
label = nothing,
)
end
plot!(graph1, tspan, n_avg, color = :red, label = nothing)
For the following cavity field amplitude we see that the ensemble average becomes zero, but the trajectories have finite cavity field amplitude.
graph2 = plot(xlabel = "Time (ms)", ylabel = "Real part cavity field")
for i in 1:traj
plot!(graph2, tspan, real(get_solution(sol.u[i], a', scaled_eqs)(tspan)), color = :grey, alpha = 0.75, label = nothing)
end
plot!(graph2, tspan, a_real_avg, color = :red, label = nothing)
Package versions
These results were obtained using the following versions:
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(
["QuantumCumulants", "ModelingToolkitBase", "OrdinaryDiffEqLowOrderRK", "StochasticDiffEqLowOrder", "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
[77a26b50] DiffEqNoiseProcess v5.36.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.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
[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.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
[19c5a474] StochasticDiffEqCore v2.2.3
[d15fe365] StochasticDiffEqLowOrder v2.0.5
[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
[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.