GitHub

This notebook can be found on github

Ramsey spectroscopy of two-level atom ensemble

Consider an ensemble of two-level atoms (spin-1/2 particles) with transition frequency $\omega_{a}$ subject to decay. The Hamiltonian of the system is

$H_{at} = \frac{\Delta}{2}\sum_i\sigma_i^z,$

where $\Delta = \omega_a - \omega_l$ and $\omega_l$ is a reference frequency. We describe the decay with rate $\gamma$ via the Lindblad term

$\mathcal{L}[\rho] = \frac{\gamma}{2}\sum_i\left(2\sigma_i^-\rho\sigma_i^+ - \sigma_i^+\sigma_i^-\rho - \rho\sigma_i^+\sigma_i^-\right).$

Say we want to perform a Ramsey interferometry scheme on this ensemble. This is done in three stages:

  1. A laser pulse (laser frequency $\omega_l$) with amplitude $\eta$ is applied for a time $T$ such that $\eta T=\frac{\pi}{4}$ ($\pi/2$-pulse).

  2. The ensemble is left to evolve freely under the Hamiltonian $H_{at}$ and the Liouvillian $\mathcal{L}[\rho]$.

  3. Finally, another $\pi/2$-pulse is applied and the total population inversion $\sum_i\langle\sigma_i^z\rangle$ is measured.

using QuantumOptics
using PyPlot
N = 2 # Atom number
γ = 0.1 # Decay rate
Δ = 0.05 # Detuning
η = 15.0 # Pulse strength

T = π/4η # Length for pulse
Tsteps = 101
Tlist = collect(range(0, stop=T, length=Tsteps))

τ = γ == 0 ? 1.0 : 0.5/γ # Length of free time evolution; 1.0 if γ=0, else 0.5/γ
τsteps = 201
τlist = collect(range(0, stop=τ, length=τsteps))

b_atom = SpinBasis(1//2)
b_coll = tensor([b_atom for i=1:N]...)
# Single atom operators
sm(i) = embed(b_coll, i, sigmam(b_atom))
sp(i) = embed(b_coll, i, sigmap(b_atom))
sz(i) = embed(b_coll, i, sigmaz(b_atom))

# Collective operators
Sx = 0.5*sum(sm.(1:N) + sp.(1:N))
Sy = 0.5im*sum(sm.(1:N) - sp.(1:N))
Sz = 0.5*sum(sz.(1:N))
H_at = Δ*Sz

The Hamiltonian of the driving laser is $H_l = \eta\sum_i\left(\sigma_i^- + \sigma_i^+\right)=2\eta S_x.$

H_l = 2η*Sx;

J = [sm(i) for i=1:N]
rates = [γ for i=1:N]
ψ₀ = tensor([spindown(b_atom) for i=1:N]...)

A convenient way to illustrate a quantum state is the so-called Bloch sphere, which can be generalized to a collection of atoms. The collective Bloch vector is defined by

$\vec{B} = \left(\langle S_x\rangle, \langle S_y\rangle, \langle S_z\rangle\right)^\mathrm{T}$.

We define a function to calculate this vector.

bloch(ρ) = [real(expect(s, ρ)) for s=[Sx, Sy, Sz]]

The Bloch vector of the initial state (all atoms in the ground state) then points downwards on the sphere.

# Plot a circle with x and y labels (the Bloch sphere from different sides)
function plot_circle(x, y)
    ϕ = collect(range(0, stop=2π, length=201))
    plot(cos.(ϕ), sin.(ϕ))
    
    # Label circle axes
    text(0.9, 0.05, x)
    text(0.05, 0.9, y)
    
    # Draw axes
    plot([0, 0], [0, 1], "k", lw=0.5)
    plot([0, 1], [0, 0], "k", lw=0.5)
end
bloch0 = bloch(dm(ψ₀))

fig = figure(figsize=(9, 3))

ax1 = subplot(131)
plot_circle("x","y")
plot([0, bloch0[1]], [0, bloch0[2]], lw=3)
ax1.axis("off")

ax2 = subplot(132)
plot_circle("x", "z")
plot([0, bloch0[1]], [0, bloch0[3]], lw=3)
ax2.axis("off")

ax3 = subplot(133)
plot_circle("y", "z")
plot([0, bloch0[2]], [0, bloch0[3]], lw=3)
ax3.axis("off")

tight_layout()

png

1) First $\pi/2$-pulse

Tout, ρT = timeevolution.master(Tlist, ψ₀, H_at + H_l, J; rates=rates)

inv_π2 = sum([real(expect(sz(i), ρT)) for i=1:N])
ρπ2 = ρT[end]; # State after first π/2-pulse
bloch_π2 = bloch(ρπ2)

fig = figure(figsize=(9, 3))

ax1 = subplot(131)
plot_circle("x","y")
plot([0, bloch_π2[1]], [0, bloch_π2[2]], lw=3)
ax1.axis("off")

ax2 = subplot(132)
plot_circle("x", "z")
plot([0, bloch_π2[1]], [0, bloch_π2[3]], lw=3)
ax2.axis("off")

ax3 = subplot(133)
plot_circle("y", "z")
plot([0, bloch_π2[2]], [0, bloch_π2[3]], lw=3)
ax3.axis("off")

tight_layout()

png

2) Free time evolution

τout, ρτ = timeevolution.master(τlist, ρπ2, H_at, J; rates=rates);

ρm = ρτ[end]; # State after free time evolution
bloch_m = bloch(ρm)

fig = figure(figsize=(9, 3))

ax1 = subplot(131)
plot_circle("x","y")
plot([0, bloch_m[1]], [0, bloch_m[2]], lw=3)
ax1.axis("off")

ax2 = subplot(132)
plot_circle("x", "z")
plot([0, bloch_m[1]], [0, bloch_m[3]], lw=3)
ax2.axis("off")

ax3 = subplot(133)
plot_circle("y", "z")
plot([0, bloch_m[2]], [0, bloch_m[3]], lw=3)
ax3.axis("off")

tight_layout()

png

3) Second $\pi$-2 pulse

Tout, ρT = timeevolution.master(Tlist, ρm, H_at + H_l, J; rates=rates)
ρπ = ρT[end]; # Final state after second π/2-pulse
bloch_π = bloch(ρπ)

fig = figure(figsize=(9, 3))

ax1 = subplot(131)
plot_circle("x","y")
plot([0, bloch_π[1]], [0, bloch_π[2]], lw=3)
ax1.axis("off")

ax2 = subplot(132)
plot_circle("x", "z")
plot([0, bloch_π[1]], [0, bloch_π[3]], lw=3)
ax2.axis("off")

ax3 = subplot(133)
plot_circle("y", "z")
plot([0, bloch_π[2]], [0, bloch_π[3]], lw=3)
ax3.axis("off")

tight_layout()

png

As you can see, for the chosen parameters we end up in a state near the fully inverted one, with an inversion of

print(bloch_π[3])
0.7475443402382681

Ramsey fringes

Now, in order to obtain the famous Ramsey fringes, we have to scan over the laser frequency, i.e. the detuning. To this end, we need to do the entire above procedure but for different detunings.

Let's write a function for the above procedure but allowing for different detunings that returns the final population inversion.

function ramsey(Δ)
    H_at = Δ*Sz
    Tout, ρ1 = timeevolution.master(Tlist, ψ₀, H_at + H_l, J; rates=rates) # First π/2-pulse
    τout, ρτ = timeevolution.master(τlist, ρ1[end], H_at, J; rates=rates) # Free time evolution
    Tout, ρ2 = timeevolution.master(Tlist, ρτ[end], H_at + H_l, J; rates=rates) # Second π/2-pulse
    real(expect(Sz, ρ2[end])) # Return resulting final inversion
end
Δ_ls = collect(range(-80, stop=80, length=501)) # List of detunings

Sz_exp = ramsey.(Δ_ls)

figure(figsize=(8, 4))
plot(Δ_ls, Sz_exp)
ylabel(L"$\langle S_z\rangle$")
xlabel(L"$\Delta$")
title("Ramsey fringes")

png