Adapted from QuantumOptics' ramsey.ipynb
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:
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).
The ensemble is left to evolve freely under the Hamiltonian $H_{at}$ and the Liouvillian $\mathcal{L}[\rho]$.
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: matplotlib
┌ Info: Installing pyqt package to avoid buggy tkagg backend. └ @ PyPlot /Users/goerz/.julia/packages/PyPlot/XaELc/src/init.jl:118
using PyCall
qutip = pyimport("qutip");
N = 2 # Atom number
γ = 0.0# 0.1 # Decay rate
Δ = 0.2#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.
function show_bloch(bloch)
b = qutip.Bloch()
b.add_vectors(bloch)
b.zlabel = ["e", "g"]
b.view = [-60, 30]
b.render()
end
show_bloch (generic function with 1 method)
bloch0 = bloch(ψ₀)
3-element Vector{Float64}: 0.0 0.0 -1.0
show_bloch(bloch0)
┌ Warning: `vendor()` is deprecated, use `BLAS.get_config()` and inspect the output instead │ caller = npyinitialize() at numpy.jl:67 └ @ PyCall /Users/goerz/.julia/packages/PyCall/3fwVL/src/numpy.jl:67
1) First $\pi/2$-pulse
H_at.data
4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries: 0.2+0.0im ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -0.2+0.0im
H_l.data
4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 8 stored entries: ⋅ 15.0+0.0im 15.0+0.0im ⋅ 15.0+0.0im ⋅ ⋅ 15.0+0.0im 15.0+0.0im ⋅ ⋅ 15.0+0.0im ⋅ 15.0+0.0im 15.0+0.0im ⋅
ψ₀
Ket(dim=4) basis: [Spin(1/2) ⊗ Spin(1/2)] 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.0 + 0.0im
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)
3-element Vector{Float64}: -0.006666603090060102 0.9999777714501946 -9.536490984507484e-6
norm(bloch_π2)
0.9999999935411014
show_bloch(bloch_π2)
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)
3-element Vector{Float64}: -0.2051986295311271 0.9787203426613192 -9.536490984507484e-6
show_bloch(bloch_m)
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(ρπ)
3-element Vector{Float64}: -0.21171423046768093 -0.0013925901987427288 0.9773306039920862
show_bloch(bloch_π)
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.9773306039920862
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);
function show_ramsey_fringes(Δ_ls, Sz_exp)
fig, ax = matplotlib.pyplot.subplots(figsize = (8, 4))
ax.plot(Δ_ls, Sz_exp)
for y in (-1, 0, 1)
ax.axhline(y, color="black", lw=0.5)
end
ax.axvline(0, color="black", lw=0.5)
ax.set_ylabel(raw"$\langle S_z\rangle$")
ax.set_xlabel(raw"$\Delta$")
ax.set_title("Ramsey fringes");
end;
show_ramsey_fringes(Δ_ls, Sz_exp);