Optomechanical Cooling

In this example, we show how to implement a cooling scheme based on radiation pressure coupling of light to a mechanical oscillator, such as a membrane. The oscillator is placed inside an optical cavity. The cavity is driven by a laser and the resulting radiation pressure of the cavity field effectively couples the photons in the cavity mode to the vibrational phonons of the mechanical oscillator mode. This model is based on the one studied in C. Genes, et. al., Phys. Rev. A 77, 033804 (2008), and the Hamiltonian reads

\[H = -\hbar\Delta a^\dagger a + \hbar\omega_m b^\dagger b + \hbar Ga^\dagger a \left(b + b^\dagger\right) + \hbar E \left(a + a^\dagger\right),\]

where $\Delta = \omega_\ell - \omega_c$ is the detuning between the driving laser ($\omega_\ell$) and the cavity ($\omega_c$). The amplitude of the laser is denoted by $E$, the resonance frequency of the mechanical oscillator by $\omega_m$, and the radiation pressure coupling is given by $G$. Additionally, photons leak out of the cavity at a rate $\kappa$. We start by loading the needed packages and specifying the model.

using QuantumCumulants
using OrdinaryDiffEqLowOrderRK, ModelingToolkitBase
using Plots

hc = FockSpace(:cavity) # Hilbertspace
hm = FockSpace(:motion)
h = hc ⊗ hm

@qnumbers a::Destroy(h, 1) b::Destroy(h, 2) # Operators


@variables Δ ωₘ E G κ # Parameters


H = -Δ * a' * a + ωₘ * b' * b + G * a' * a * (b + b') + E * (a + a') # Hamiltonian


J = [a] # Jump operators & rates
rates = [κ]

We are specifically interested in the average number of photons $\langle a^\dagger a \rangle$ and phonons $\langle b^\dagger b \rangle$. Thus, we first derive the equations for these two averages. We restrict our description to a second order cumulant expansion.

ops = [a' * a, b' * b] # Derive equations
eqs = meanfield(ops, H, J; rates = rates, order = 2)

\[ \begin{aligned} \partial_{t} \langle a^{\dagger}a \rangle &= - \langle a^{\dagger}a \rangle \kappa + E \langle a \rangle i - E \langle a^{\dagger} \rangle i \\\\[-0.0em] \partial_{t} \langle b^{\dagger}b \rangle &= G \left( \langle a^{\dagger} \rangle \langle ab \rangle + \langle b \rangle \langle a^{\dagger}a \rangle + \langle a^{\dagger}b \rangle \langle a \rangle - 2 \langle a^{\dagger} \rangle \langle a \rangle \langle b \rangle \right) i - G \left( \langle a \rangle \langle a^{\dagger}b^{\dagger} \rangle + \langle b^{\dagger} \rangle \langle a^{\dagger}a \rangle + \langle ab^{\dagger} \rangle \langle a^{\dagger} \rangle - 2 \langle a^{\dagger} \rangle \langle a \rangle \langle b^{\dagger} \rangle \right) i \end{aligned} \]

To get a closed set of equations we automatically complete the system.

eqs_completed = complete!(deepcopy(eqs)) # Complete equations

\[ \begin{aligned} \partial_{t} \langle a^{\dagger}a \rangle &= - \langle a^{\dagger}a \rangle \kappa + E \langle a \rangle i - E \langle a^{\dagger} \rangle i \\\\[-0.0em] \partial_{t} \langle b^{\dagger}b \rangle &= G \left( \langle a^{\dagger} \rangle \langle ab \rangle + \langle b \rangle \langle a^{\dagger}a \rangle + \langle a^{\dagger}b \rangle \langle a \rangle - 2 \langle a^{\dagger} \rangle \langle a \rangle \langle b \rangle \right) i - G \left( \langle a \rangle \langle a^{\dagger}b^{\dagger} \rangle + \langle b^{\dagger} \rangle \langle a^{\dagger}a \rangle + \langle ab^{\dagger} \rangle \langle a^{\dagger} \rangle - 2 \langle a^{\dagger} \rangle \langle a \rangle \langle b^{\dagger} \rangle \right) i \\\\[-0.0em] \partial_{t} \langle a \rangle &= - E i - G \langle ab \rangle i - G \langle ab^{\dagger} \rangle i + \langle a \rangle \left( - \frac{1}{2} \kappa + i \Delta \right) \\\\[-0.0em] \partial_{t} \langle b^{\dagger} \rangle &= G \langle a^{\dagger}a \rangle i + \langle b^{\dagger} \rangle i \mathtt{\omega_m} \\\\[-0.0em] \partial_{t} \langle a^{\dagger}b^{\dagger} \rangle &= E \langle b^{\dagger} \rangle i + G \langle a^{\dagger} \rangle i + \langle a^{\dagger}b^{\dagger} \rangle \left( - \frac{1}{2} \kappa + i \left( - \Delta + \mathtt{\omega_m} \right) \right) + G \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) i + G \left( \langle b^{\dagger}b^{\dagger} \rangle \langle a^{\dagger} \rangle + 2 \langle b^{\dagger} \rangle \langle a^{\dagger}b^{\dagger} \rangle - 2 \langle b^{\dagger} \rangle^{2} \langle a^{\dagger} \rangle \right) i + G \left( \langle a^{\dagger} \rangle \langle b^{\dagger}b \rangle + \langle b \rangle \langle a^{\dagger}b^{\dagger} \rangle + \langle a^{\dagger}b \rangle \langle b^{\dagger} \rangle - 2 \langle a^{\dagger} \rangle \langle b^{\dagger} \rangle \langle b \rangle \right) i \\\\[-0.0em] \partial_{t} \langle ab^{\dagger} \rangle &= - E \langle b^{\dagger} \rangle i + \langle ab^{\dagger} \rangle \left( - \frac{1}{2} \kappa + i \left( \Delta + \mathtt{\omega_m} \right) \right) + G \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) i - G \left( \langle b^{\dagger}b^{\dagger} \rangle \langle a \rangle + 2 \langle ab^{\dagger} \rangle \langle b^{\dagger} \rangle - 2 \langle b^{\dagger} \rangle^{2} \langle a \rangle \right) i - G \left( \langle a \rangle \langle b^{\dagger}b \rangle + \langle ab \rangle \langle b^{\dagger} \rangle + \langle ab^{\dagger} \rangle \langle b \rangle - 2 \langle a \rangle \langle b^{\dagger} \rangle \langle b \rangle \right) i \\\\[-0.0em] \partial_{t} \langle a^{\dagger}a^{\dagger} \rangle &= 2 E \langle a^{\dagger} \rangle i + \langle a^{\dagger}a^{\dagger} \rangle \left( - \kappa - 2 i \Delta \right) + 2 G \left( \langle a^{\dagger}a^{\dagger} \rangle \langle b \rangle + 2 \langle a^{\dagger} \rangle \langle a^{\dagger}b \rangle - 2 \langle a^{\dagger} \rangle^{2} \langle b \rangle \right) i + 2 G \left( \langle a^{\dagger}a^{\dagger} \rangle \langle b^{\dagger} \rangle + 2 \langle a^{\dagger} \rangle \langle a^{\dagger}b^{\dagger} \rangle - 2 \langle a^{\dagger} \rangle^{2} \langle b^{\dagger} \rangle \right) i \\\\[-0.0em] \partial_{t} \langle b^{\dagger}b^{\dagger} \rangle &= 2 \langle b^{\dagger}b^{\dagger} \rangle i \mathtt{\omega_m} + 2 G \left( \langle a \rangle \langle a^{\dagger}b^{\dagger} \rangle + \langle b^{\dagger} \rangle \langle a^{\dagger}a \rangle + \langle ab^{\dagger} \rangle \langle a^{\dagger} \rangle - 2 \langle a^{\dagger} \rangle \langle a \rangle \langle b^{\dagger} \rangle \right) i \end{aligned} \]

To calculate the dynamics we create a system of ordinary differential equations, which can be used by DifferentialEquations.jl.

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

Finally, we need to define the numerical parameters and the initial state of the system. We will consider the membrane at room temperature. Its vibrational mode is in a thermal state with an average number of phonons that can be estimated from $k_B T = n_\mathrm{vib}\hbar \omega_m$. If the resonator has a resonance frequency of $\omega_m = 10\mathrm{MHz}$, then the number of phonons at room temperature ($T\approx 300K$) is approximately $n_\mathrm{vib} \approx 4\times 10^6$.

u0 = initial_values(eqs_completed; defaults = Dict(average(b' * b) => 4.0e6 + 0im)) # Initial state (4e6 phonons)

p0 = Dict{Num, ComplexF64}(Δ => -10.0 + 0im, ωₘ => 1.0 + 0im, E => 200.0 + 0im, G => 0.0125 + 0im, κ => 20.0 + 0im) # System parameters
prob = ODEProblem(sys_c, merge(u0, p0), (0.0, 60000.0))
sol = solve(prob, RK4())

t = real.(sol.t) # Plot results
phonons = real.(get_solution(sol, b'b, eqs_completed).(sol.t))
T = 7.5e-5 * phonons
photons = real.(get_solution(sol, a'a, eqs_completed).(sol.t))

p1 = plot(t, T, ylabel = "T in K", legend = false)
p2 = plot(t, photons, xlabel = "t⋅ωₘ", ylabel = "⟨a⁺a⟩", legend = false)
plot(p1, p2, layout = (2, 1), size = (650, 400))

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 × 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
  [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
  [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.