# The Analytical Solution of Rabi Cycling in the Two-Level-System¶

In [1]:
import matplotlib.pylab as plt
import numpy as np
from numpy import pi as π
from numpy import exp, cos, sin, sqrt, abs

In [2]:
%matplotlib notebook

In [3]:
from IPython.display import display

In [4]:
%load_ext version_information
%version_information numpy, matplotlib, sympy

Out[4]:
SoftwareVersion
Python3.7.3 64bit [GCC 7.3.0]
IPython7.8.0
OSLinux 4.14.138+ x86_64 with debian buster sid
numpy1.18.1
matplotlib3.1.2
sympy1.5.1
Tue Jan 28 23:34:30 2020 UTC

## Check the analytical solution to the TLS¶

We consider the solution of the Schrödinger equation for the two level system in the rotating wave approximation $$i \hbar \frac{\partial}{\partial t} \begin{pmatrix} b(t) \\\\ a(t) \end{pmatrix} = \begin{pmatrix} 0 & \frac{\Omega_0}{2} \\\\ \frac{\Omega_0}{2} & \Delta\end{pmatrix} \begin{pmatrix} b(t) \\\\ a(t) \end{pmatrix}$$ with the boundary condition $a(0) = 1$, $b(0) = 0$. Here and everwhere else, we work in energy units of $2 \pi$ GHz and time units of ns, and thus have $\hbar = 1$

In [5]:
import sympy as sp
from sympy import I, diff
sp.init_printing()


We "guess" the analytical solutions (based on D. Tannor, "Introduction to Quantum Mechanics: A Time-dependent Perspective", with some small corrections):

In [6]:
Δ, Ω0, t, Ω = sp.symbols((r'\Delta', r'\Omega_0', 't', r'\Omega'), real=True)
#Ω = sp.sqrt(Δ**2 + Ω0**2)
a = sp.exp(-I*Δ*t/2)*( sp.cos(Ω*t/2) - I*(Δ/Ω)*sp.sin(Ω*t/2) )
b = -I * (Ω0/Ω)*sp.exp(-I*Δ*t/2)*sp.sin(Ω*t/2)

In [7]:
a

Out[7]:
$\displaystyle \left(- \frac{i \Delta \sin{\left(\frac{\Omega t}{2} \right)}}{\Omega} + \cos{\left(\frac{\Omega t}{2} \right)}\right) e^{- \frac{i \Delta t}{2}}$
In [8]:
b

Out[8]:
$\displaystyle - \frac{i \Omega_0 e^{- \frac{i \Delta t}{2}} \sin{\left(\frac{\Omega t}{2} \right)}}{\Omega}$

and check that the left hand side of the Schrödinger Equation equals the right side

In [9]:
sp.simplify((diff(a, t) + I*(Ω0/2)*b + I*Δ*a).subs({Ω:sp.sqrt(Δ**2 + Ω0**2)}))

Out[9]:
$\displaystyle 0$
In [10]:
sp.simplify(diff(b, t) + I*(Ω0/2)*a)

Out[10]:
$\displaystyle 0$

## Ideal π-pulse¶

Population transfer is achieved for $\Omega t = \pi$. The amplitudes in the ideal case of no detuning are:

In [11]:
a.subs({Ω * t: sp.pi, Δ:0})

Out[11]:
$\displaystyle 0$
In [12]:
b.subs({Ω * t: sp.pi, Δ:0, Ω0: Ω})

Out[12]:
$\displaystyle - i$

Note the non-zero phase on $b$: Even a $2 \pi$ pulse (a full cycle) will not reproduce the initial state, but result in a $\pi$-phase shift. Only two full cycles reproduce the exact original state.

For $\Delta \neq 0$, the phases behave in a considerably more complicated way. In particular, the imaginary part of $a$ (i.e., the phase of $\vert0\rangle$) will no longer be zero:

In [13]:
a.subs({Ω * t: sp.pi})

Out[13]:
$\displaystyle - \frac{i \Delta e^{- \frac{i \Delta t}{2}}}{\Omega}$
In [14]:
b.subs({Ω * t: sp.pi})

Out[14]:
$\displaystyle - \frac{i \Omega_0 e^{- \frac{i \Delta t}{2}}}{\Omega}$

## Compare analytical and numerical solution of TLS¶

We now simulate the dynamics of the two-level system numerically, and compare the result to the analytical solution

In [15]:
Ω0 = 2*π*0.5 # 2π GHz
Δ = 2*π*0.1 # 2π GHz

In [16]:
def a(Ω0, Δ, t):
"""Complex amplitude of the TLS solution, for the \ket{1} component"""
Ω = sqrt(Δ**2 + Ω0**2)
return exp(-0.5j*Δ*t)*( cos(0.5*Ω*t) - 1j*(Δ/Ω)*sin(0.5*Ω*t) )

In [17]:
def b(Ω0, Δ, t):
"""Complex amplitude of the TLS solution, for the \ket{0} component"""
Ω = sqrt(Δ**2 + Ω0**2)
return -1j*(Ω0/Ω)*exp(-0.5j*Δ*t)*sin(0.5*Ω*t)

In [18]:
def show_rabi_cycling(Ω0, Δ, tgrid):
fig = plt.figure()
a_s = a(Ω0, Δ, tgrid)
b_s = b(Ω0, Δ, tgrid)
ax.plot(tgrid, abs(a_s)**2, label='pop(1)', color='red')
ax.plot(tgrid, abs(b_s)**2, label='pop(0)', color='blue', ls='--')
fig.show()

In [19]:
show_rabi_cycling(Ω0, Δ, np.linspace(0, 10, 1000))

In [20]:
H = np.matrix([[0, 0.5*Ω0],[0.5*Ω0, Δ]], dtype=np.complex128); H

Out[20]:
matrix([[0.        +0.j, 1.57079633+0.j],
[1.57079633+0.j, 0.62831853+0.j]])
In [21]:
evals, W = np.linalg.eig(H)

In [22]:
def get_U(evals, W, dt):
"""Get the time evolution operator for a single time step"""
D = np.matrix(np.diag(np.exp(-1j*evals*dt)))
W = np.matrix(W)
return np.asarray(W*D*W.H)

In [23]:
def propagate(psi, tgrid, evals, W):
"""Propagate psi along the given time grid, assuming a constant
Hamiltonian with the given eigenvalues and eigenvector matrix.

Returns an array of shape (2,nt) that gives the propagated
amplitude at each point in time."""
nt = len(tgrid)
result = np.zeros((2, nt), dtype=np.complex128)
tgrid = list(tgrid)
dt = tgrid[1] - tgrid[0]
i = 0
U = get_U(evals, W, dt)
while len(tgrid) > 0:
tgrid.pop()
result[:, i] = psi
psi = U.dot(psi)
i += 1
return result

In [24]:
dyn = propagate(np.array([0,1]), tgrid=np.linspace(0,10,1000), evals=evals, W=W)

In [25]:
def show_dyn(dyn, tgrid):
fig = plt.figure()
ax.plot(tgrid, abs(dyn[1,:])**2, label='pop(1)', color='red')
ax.plot(tgrid, abs(dyn[0,:])**2, label='pop(0)', color='blue', ls='--')
fig.show()

In [26]:
show_dyn(dyn, np.linspace(0, 10, 1000))


The maximum error between the analytical and numerical solution is:

In [27]:
print(np.max(np.abs(dyn[0,:] - b(Ω0, Δ, np.linspace(0, 10, 1000)))))

2.29480091309186e-13

In [28]:
print(np.max(np.abs(dyn[1,:] - a(Ω0, Δ, np.linspace(0, 10, 1000)))))

2.7207552755516735e-13