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
%reload_ext version_information
%version_information numpy, matplotlib, sympy
Installed version_information.py. To use it, type:
  %load_ext version_information
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 \begin{equation} 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} \end{equation} 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()
    ax = fig.add_subplot(111)
    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 = fig.add_subplot(111)
    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