The Jaynes-Cummings Hamiltonian models the dynamics of a qubit, or two-level system, with a harmonic oscillator. We can write it using three terms: one for the harmonic oscillator (resonator), of for the spin-1/2 (qubit) and an interaction term:
$$ H = \hbar \omega_r a^\dagger a - \frac{1}{2} \hbar \omega_q \sigma_z + \hbar g (a + a^\dagger)(\sigma_+ + \sigma_-) $$Note that this equation is valid before the rotating wave approximation.
In this tutorial, we explore the ground states of the Hamiltonian above when varying the strenght of the coupling g. We start by loading the JaynesCummings package, and a plotting package.
using JaynesCummings, PlotlyJS
Plotly javascript loaded.
To load again call
init_notebook(true)
We now define a few parameters. Because the Hilbert space of a quantum harmonic oscillator is infinite dimensional, we must truncate it to a finite value. You should always set this value to a number higher than the highest anticipated resonator state. In this case, 50 is enough. ω_q
and ω_r
are the qubit and resonator frequencies. Their value is not important in this case.
We also define a range of coupling to determine the effect of varying $g$.
N = 120
ω_q = 1.0
ω_r = 1.0
w_samples = 100
gnumber = 251
couplings = linspace(0,2.5,gnumber);
We now find the groundstate of the Hamiltonian for each given coupling. This involves generating the Hamiltonian, and finding its eigenvalues and eigenvectors. Because the eig
function returns the vectors sorted from lowest to highest eigenvalue (energy), the groundstate is simply the first vector.
For each groundstate ψ
, we calculate the expectation value of $a^\dagger a$ and $\sigma_+ \sigma_-$. These quantities represent the average number of excitations in either system.
resonator_excitation = Array{Float64}(gnumber)
qubit_excitation = Array{Float64}(gnumber)
groundstate_energy = Array{Float64}(gnumber)
adaga = qeye(2)⊗(a_dagger(N)*a(N))
σpσm = (σ_plus*σ_minus)⊗qeye(N)
for (i,g) in enumerate(couplings)
H = gen_hamiltonian(N,ω_q,ω_r,g,rwa=false)
evals, ekets = eig(H)
groundstate_energy[i] = evals[1]
ψ = ekets[:,1]
resonator_excitation[i] = real(ψ'*adaga*ψ)[1]
qubit_excitation[i] = real(ψ'*σpσm*ψ)[1]
end
We can now plot the average resonator and qubit occupation as a function of the coupling g
.
plot([scatter(;x=couplings,y=resonator_excitation,name="Resonator Occupation"),
scatter(;x=couplings,y=qubit_excitation,name="Qubit Occupation")],
Layout(;xaxis=attr(title="Coupling")))
We can see that at low interaction strenght, the qubit and resonator are both in their respective ground states. As the coupling is increased however, there is an energetic advantage to having a non-zero occupation number. We can notice this by plotting the groundstate energy as a function of coupling:
plot(scatter(;x=couplings,y=groundstate_energy/real(ħ)),Layout(;xaxis=attr(title="Coupling"),yaxis=attr(title="Energy (ħ)")))
The energy drops monotonically, starting from $-\hbar/2$, the energy of the qubit in the ground state.
We can examine the density matrix of the reduced resonator system at the highest coupling by tracing out the qubit from the full density matrix:
H = gen_hamiltonian(N,ω_q,ω_r,couplings[end],rwa=false)
_, ekets = eig(H)
ρ = ekets[:,1]*ekets[:,1]'
ρ_resonator = ptrace(ρ,[2,N],1)
plot_densitymatrix(ρ_resonator)
W = calc_wignerfunction_resonator(ρ_resonator,w_samples,maxdisp=7.5);
plot_wignerfunction3d(W,maxdisp=7.5)