# Rabi Cycling in the Two-Level-System — From Scratch¶

This notebook uses sympy to symbolically solve the Schrödinger equation for a two-level system, resulting in the exact analytic expressions for the complex amplitudes and the population dynamics of Rabi cycling.

Having this written out in a computer algebra system allows to easily adapt the calculation to variations in the form of the two-level-system Hamiltonian, or different boundary conditions.

In [1]:
from sympy import Function, symbols, sqrt, Eq, Abs, exp, cos, sin, Matrix, dsolve, solve

In [2]:
from sympy import I as 𝕚

In [3]:
g = Function('g')
e = Function('e')
t = symbols('t', positive=True)
Δ = symbols('Delta', real=True)
Ω0 = symbols('Ω_0', real=True)
Ω = symbols('Omega', real=True)

In [4]:
Ĥ = Matrix([
[  0  , -Ω0/2],
[-Ω0/2,    Δ ]
])
Ĥ

Out[4]:
$\displaystyle \left[\begin{matrix}0 & - \frac{Ω_{0}}{2}\\- \frac{Ω_{0}}{2} & \Delta\end{matrix}\right]$
In [5]:
TDSE_system = [
g(t).diff(t) + 𝕚 * (Ĥ[0,0] * g(t) + Ĥ[0,1] * e(t)),
e(t).diff(t) + 𝕚 * (Ĥ[1,0]*  g(t) + Ĥ[1,1] * e(t)),
]

In [6]:
sols_gen = dsolve(TDSE_system, [g(t), e(t)])


Note that dsolve allows to specify boundary conditions, but the resulting expressions are hard to simply. We do better by solving for the integration constants "manually" later on.

In [7]:
effective_rabi_freq = {
Ω: sqrt(Δ**2 + Ω0**2)
}

In [8]:
Eq(Ω, effective_rabi_freq[Ω])

Out[8]:
$\displaystyle \Omega = \sqrt{\Delta^{2} + Ω_{0}^{2}}$
In [9]:
find_Ω = {
effective_rabi_freq[Ω]: Ω
}

In [10]:
sols_gen[0].subs(find_Ω)

Out[10]:
$\displaystyle g{\left(t \right)} = - \frac{C_{1} Ω_{0} e^{- \frac{t \left(i \Delta - i \Omega\right)}{2}}}{\Delta - \Omega} - \frac{C_{2} Ω_{0} e^{- \frac{t \left(i \Delta + i \Omega\right)}{2}}}{\Delta + \Omega}$
In [11]:
sols_gen[1].subs(find_Ω).expand()

Out[11]:
$\displaystyle e{\left(t \right)} = C_{1} e^{- \frac{i \Delta t}{2}} e^{\frac{i \Omega t}{2}} + C_{2} e^{- \frac{i \Delta t}{2}} e^{- \frac{i \Omega t}{2}}$
In [12]:
boundary_conditions = {
t: 0,
g(t): 1,
e(t) : 0,
}

In [13]:
find_Ω0sq = {
Ω**2 - Δ**2: Ω0**2
}

In [14]:
C1, C2 = symbols("C1, C2")

In [15]:
_integration_constants = solve(
[sol.subs(find_Ω).subs(boundary_conditions) for sol in sols_gen],
[C1, C2]
)
integration_constants = {
k: v.subs(find_Ω0sq)
for (k, v) in _integration_constants.items()
}

In [16]:
Eq(C1, integration_constants[C1])

Out[16]:
$\displaystyle C_{1} = \frac{Ω_{0}}{2 \Omega}$
In [17]:
Eq(C2, integration_constants[C2])

Out[17]:
$\displaystyle C_{2} = - \frac{Ω_{0}}{2 \Omega}$
In [18]:
global_phase = exp(-𝕚 * Δ * t / 2)

In [19]:
def remove_global_phase(eq, global_phase=global_phase):
return Eq(eq.lhs, eq.rhs * global_phase.conjugate())

def restore_global_phase(eq, global_phase=global_phase):
return Eq(eq.lhs, eq.rhs * global_phase)

In [20]:
def collect(eq, term):
return Eq(eq.lhs, eq.rhs.collect(term))

In [21]:
def simplify_sol_g(sol_g):
rhs = (
sol_g
.rhs
.rewrite(sin)
.expand()
.collect(𝕚)
.subs(effective_rabi_freq)
.simplify()
.subs(find_Ω)
)
return Eq(sol_g.lhs, rhs)

In [22]:
sol_g = restore_global_phase(
simplify_sol_g(
collect(
remove_global_phase(
sols_gen[0]
.subs(find_Ω)
.subs(integration_constants)
).expand(),
exp(𝕚 * Ω * t / 2                                                                   )
)
)
)
sol_g

Out[22]:
$\displaystyle g{\left(t \right)} = \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 [23]:
sol_e = restore_global_phase(
remove_global_phase(
sols_gen[1].subs(find_Ω).subs(integration_constants)
).expand().rewrite(sin).expand()
)
sol_e

Out[23]:
$\displaystyle e{\left(t \right)} = \frac{i Ω_{0} e^{- \frac{i \Delta t}{2}} \sin{\left(\frac{\Omega t}{2} \right)}}{\Omega}$
In [24]:
pop_e = (sol_e.rhs * sol_e.rhs.conjugate())
Eq(Abs(e(t))**2, pop_e)

Out[24]:
$\displaystyle \left|{e{\left(t \right)}}\right|^{2} = \frac{Ω_{0}^{2} \sin^{2}{\left(\frac{\Omega t}{2} \right)}}{\Omega^{2}}$