# 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]:
%install_ext http://raw.github.com/goerz/version_information/master/version_information.py
%version_information numpy, matplotlib, sympy

Installed version_information.py. To use it, type:

Out[4]:
SoftwareVersion
Python3.3.5 64bit [GCC 4.2.1 (Apple Inc. build 5577)]
IPython3.0.0
OSDarwin 14.5.0 x86_64 i386 64bit
numpy1.9.2
matplotlib1.4.3
sympy0.7.6
Mon Aug 24 12:15:32 2015 EDT

## 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]:
$$\left(- \frac{i \Delta}{\Omega} \sin{\left (\frac{\Omega t}{2} \right )} + \cos{\left (\frac{\Omega t}{2} \right )}\right) e^{- \frac{i \Delta}{2} t}$$
In [8]:
b

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

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]:
$$0$$
In [10]:
sp.simplify(diff(b, t) + I*(Ω0/2)*a)

Out[10]:
$$0$$

## 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 [11]:
Ω0 = 2*π*0.5 # 2π GHz
Δ = 2*π*0.1 # 2π GHz

In [12]:
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 [13]:
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 [14]:
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='--')
plt.show(fig)

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

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

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

In [18]:
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 [19]:
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 [20]:
dyn = propagate(np.array([0,1]), tgrid=np.linspace(0,10,1000), evals=evals, W=W)

In [21]:
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='--')
plt.show(fig)

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


The maximum error between the analytical and numerical solution is:

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

2.7979669322e-13

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

3.28802706295e-13