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 \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]:
$\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()
    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='--')
    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 = 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='--')
    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