Cavity cooling

This notebook can be found on github

Cavity cooling of a two-level atom

Cavity cooling relies on spontaneous emission of an atom that is coupled to a field that has a certain energy mismatch. The atomic motion then compensates for the missing energy.

Consider a two-level atom moving in a coherently driven cavity. The Hamiltonian is

$H = -\Delta_c a^\dagger a - \Delta_a\sigma^+\sigma^- + \eta(a + a^\dagger) + g\cos(x)(a^\dagger \sigma^- + a\sigma^+),$

where $\Delta_c = \omega_l - \omega_c$, $\Delta_a=\omega_l - \omega_a$ and $\omega_c$ is the cavity frequency, $\omega_a$ the atomic transition frequency and $\omega_l$ the frequency of the laser driving the cavity. The pump strength is $\eta$, while the atom couples to the cavity with $g$ and the cavity has a mode function $\cos(x)$. The atomic position $x$ is in units of the inverse cavity wave number. Note that this is a number rather than an operator since we treat the atomic motion classically.

Decay and cavity damping are modeled by the Liouvillians

$\mathcal{L}_c[\rho] = \kappa(2a\rho a^\dagger - a^\dagger a \rho - \rho a^\dagger a),$

$\mathcal{L}_a[\rho] = \gamma(2\sigma^-\rho\sigma^+ - \sigma^+\sigma^-\rho - \rho\sigma^+\sigma^-)$.

In addition we have the classical differential equations of the atomic motion obtained by making use of the Eherenfest theorem,

$\dot{x} = 2\omega_r p$,

$\dot{p} = -\partial_x\langle H\rangle = 2\sin(x)~\Re\left\{\langle a^\dagger\sigma^-\rangle\right\},$

where $p$ is the momentum in units of the cavity wave number and $\omega_r=k_c^2/(2m)$ is the recoil frequency.

Note, that in order for the cavity to cool the atom, a cavity photon must have less energy than a photon on resonance with the atom, i.e. the cavity must be red-detuned from the atom ($\omega_c < \omega_a$). We will now use the implemented semi-classical master equation to solve the above dynamics.

First, we define the necessary parameters, our Hilbert space and corresponding operators.

using QuantumOptics

# Parameters
Nc = 16
γ = 1.
g = γ/2.
κ = 0.5γ
ωr = .15γ
Δc = -γ
Δa = -2γ
η = γ
tmax = 400
tsteps = 10*tmax
dt = tmax/tsteps
tlist = [0:dt:tmax;]

# Hilbert space
bc = FockBasis(Nc)
ba = SpinBasis(1//2)

# Operators
a = destroy(bc) ⊗ one(ba)
sm = one(bc) ⊗ sigmam(ba)
sp = dagger(sm)

We now define the Hamiltonian as two separate parts, one that is position dependent and one that is independent of $x$. This is for reasons of performance, since it is more efficient to simply add the two parts where one is multiplied by $\cos(x)$ than creating the entire Hamiltonian in every time-step.

# Hamiltonian
Hx = g*(a*sp + ad*sm) # ∝ cos(x)

# Jump operators
J = [sqrt(2κ)*a, sqrt(2γ)*sm]
Jdagger = map(dagger, J)

The semi-classical time evolution requires two functions as arguments: one function that returns the updated Hamiltonian and jump operators at every step and one that updates the vector containing the time derivatives of the classical variables. Note, that the efficiency of the operations performed inside these function is very important for overall perfomance, each function is called at every time-step.

function fquantum(t, psi, u) # psi is the quantum, u the classical part
x = u[1]
return H0 + Hx*cos(x), J, Jdagger
end

adsm = ad*sm # Define to avoid doing a multiplication at every step
function fclassical(t, psi, u, du) # du is a vector containing the increments of the classical variables
du[1] = 2ωr*u[2]
end

Before we calculate the time evolution, we need to define the initial state which needs to be a semi-classical one. It consists of a quantum part (a ket or density operator) and a vector of classical variables. Note: the vector containing the classical variables needs to be complex.

x0 = sqrt(2) # Some arbitrary initial position
p0 = 7.0 # Some arbitrary initial momentum
u0 = complex.([x0, p0])

ψ0 = fockstate(bc, 0) ⊗ spindown(ba)

ψsc0 = semiclassical.State(ψ0, u0)

Finally, we can calculate the time evolution by calling dynamic master equation solver from the semiclassical module.

tout, ρt = semiclassical.master_dynamic(tlist, ψsc0, fquantum, fclassical)

We can now calculate expectation values as with any other density operator and retrieve the classical variables.

x = []
p = []
for r=ρt
push!(x, real(r.classical[1]))
push!(p, real(r.classical[2]))
end

pe = real(expect(sp*sm, ρt))
using PyPlot

figure(figsize=(10, 5))
subplot(221)
plot(tout, x)
xlabel(L"$t$")
ylabel(L"$x$")

subplot(222)
plot(tout, p.^2)
xlabel(L"$t$")
ylabel(L"$E_{kin}$")

subplot(223)
plot(tout, n)
xlabel(L"$t$")
ylabel(L"$n$")

subplot(224)
plot(tout, pe)
xlabel(L"$t$")
ylabel(L"$\langle \sigma^+\sigma^-\rangle$")

tight_layout()