We study population dynamics and effects of optical pumping on spectral line singal for near-resonant cyclic $D_{2}$ line transitions in Cs and Na atoms with the help of simplified four-level system. A classical laser field of frequency $ω_{L}$ is in resonance between highest energy ground state hyperfine component $\left|2\right>$ highest energy hyperfine state component $\left|4\right>$. Intensity of excitation laser field follows a Guassian pulse shape that is characterize by Rabi frequency coupling strengths $Ω_{24}(t)$ and $Ω_{23}(t)$ between states $\left|2\right> ↔ \left|4\right>$ and $\left|2\right> ↔ \left|3\right>$ respectively. States $\left|1\right>$ is a population trapping ground state (dark state) that does not interact with the excitation laser field.
In [1]:
%%svg
4-level-schematics.svg
Corresponding $D_{2}$ line transition hyperfine level schematics in Na nd Cs atoms.
In [2]:
%%svg
hyperfine-level-shematics.svg
The standard approach for deriving the equation of motion for a system interacting with its environment is to expand the scope of the system to include the environment. The combined quantum system is closed and its evolution is governed by the von Neumann equation \begin{equation} \begin{aligned} \dot{ρ}(t) &=-\frac{i}{ℏ}[H_{tot},ρ_{tot}(t)]\\ H_{tot} &= H_{sys} + H_{env} + H_{int}\text{,} \end{aligned} \end{equation} where the total Hamiltonian $H_{tot}$ includes the original system Hamiltonian $H_{sys}$, the Hamiltonian for the environment $H_{env}$ and the interaction Hamiltonian $H_{int}$ between the system and its environment. To obtain the dynamics of system $H_{sys}$, we can perform a partial trace over the environmental degrees of freedom in von Neumann equation. The most general trace-preserving and completely positive form of this evolution is the Lindblad master equation [[Lindblad(1976)]][communmathphys.48.119], [[Gardiner and Zoller(2004)]][Gardiner_2004], [[Walls and Milburn(2008)]][Walls_2008] \begin{equation} \begin{aligned} \dot{ρ}(t) &=-\frac{i}{ℏ}[H(t),ρ(t)] + \sum_n \frac{1}{2} \left(2 C_n \rho(t) C_n^{\dagger} - ρ(t) C_n^{\dagger} C_n - C_n^{\dagger} C_n ρ(t)\right)\\ H(t) &= H_{sys} + H_{int}\\ C_n &= \sqrt{γ_{n}} A_{n}\text{,} \end{aligned} \end{equation} where $C_n$ are wave function collapse operators, $A_{n}$ are operators through which the environment couples to the system in $H_{int}$ and $γ_{n}$ are the corresponding decay rates.
Let us obtain optical Bloch equation for four-level system. Using dressed state notation given in figure, system Hamiltonian can be written as \begin{equation} \begin{aligned} H_{sys} &= ε_{1} \left|1,n+1\right>\left< 1,n+1\right| + ε_{2} \left|2,n+1\right>\left< 2,n+1\right| + ε_{3} \left|3,n\right>\left< 3,n\right| + ε_{4} \left|4,n\right>\left< 4,n\right|\\ &= \left[\begin{smallmatrix} ε_{1} & 0 & 0 & 0\\0 & ε_{2} & 0 & 0\\0 & 0 & ε_{3} & 0\\0 & 0 & 0 & ε_{4} \end{smallmatrix}\right]\text{.} \end{aligned} \end{equation} And energies of the system Hamiltonain can be written as \begin{equation} \begin{aligned} ε_{1} &= -ℏ ω_{21} & ε_{2} &= 0 & ε_{3} &= Δ_{23} = -ω_{43} - Δ & ε_{4} &= Δ_{24} = - Δ \text{,} \end{aligned} \end{equation} where $ω_{ab} = ΔE_{ab}/ℏ = (E_{a} - E_{b})/ℏ$ and energy level difference in SI unit system using frequency $ν$ and cyclic frequency $ω$ are written as $ΔE = h ν = ℏ ω$.
Corresponding four-level system decay channels are following, $A_{1}$ from state $\left|4,n\right>$ to state $\left|2,n+1\right>$, $A_{2}$ from state $\left|3,n\right>$ to state $\left|2,n+1\right>$, $A_{3}$ from state $\left|3,n\right>$ to state $\left|1,n+1\right>$, such that \begin{equation} \begin{aligned} A_{1} &= \left|2,n+1\right>\left< 4,n\right| = \left[\begin{smallmatrix}0 & 0 & 1 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{smallmatrix}\right] \\ A_{2} &= \left|2,n+1\right>\left< 3,n\right| = \left[\begin{smallmatrix}0 & 0 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{smallmatrix}\right] \\ A_{3} &= \left|1,n+1\right>\left< 3,n\right| = \left[\begin{smallmatrix}0 & 0 & 0 & 0\\0 & 0 & 0 & 1\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{smallmatrix}\right]\text{.} \end{aligned} \end{equation} Corresponding wave function collapse operators are \begin{equation} \begin{aligned} C_{1} &= \sqrt{γ_{1}} A_{1} & C_{2} &= \sqrt{γ_{2}} A_{2} & C_{3} &= \sqrt{γ_{3}} A_{3} \\ γ_{1} &= Π_{42} Γ_{4} & γ_{2} &= Π_{32} Γ_{3} & γ_{3} &= Π_{31} Γ_{3}\\ Π_{42} &= 1 & & & Π_{31} &= 1-Π_{32} \text{,} \end{aligned} \end{equation} where $Γ_{4} = Γ_{3} = Γ = 1/τ_{e}$ is decay rate of excited states, $τ_{e}$ is natural lifetime of excited state and $Π_{F_{e} F_{g}}$ is branching ration for transition $\left|F_{e}\right> → \left|F_{g}\right>$.
The interaction of bound particles with laser light most often originates with the electric-dipole interaction. Therefore, the atom-laser interaction Hamiltonian $H_{int}$ in the dipole approximation is given by interaction energy operator, that is the projection of the electric dipole moment $\vec{d}$ onto the electric field \begin{equation} \begin{aligned} H_{int} &= -\vec{d} · \vec{E}(t)\\ \vec{E}(t) &= \hat{\mathbf{e}} E(t)\cos\left(ω_{L} t - φ\right)\\ E(t) &= E_{0} \exp\left(-2\left(\frac{t}{τ_{tr}}\right)^{2}\right) \end{aligned} \end{equation} and where the atomic dipole operator $\vec{d}$ is given in terms of the atomic electron position $\vec{r}_{e}$ as \begin{equation} \vec{d} = - \boldsymbol{e} \vec{r}_{e} \text{,} \end{equation} where we denote the fundamental charge by $\boldsymbol{e}$, so that the electron charge is $q = - \boldsymbol{e}$. The dipole transition moment between states $\left|ψ_{a}\right>$ and $\left|ψ_{b}\right>$, projected onto the field unit vector $\hat{\mathbf{e}}$, is \begin{equation} d_{ψ_{a}, ψ_{b}} = \left< ψ_{a}\,\middle|\, \vec{d} · \hat{\mathbf{e}}\,\middle|\,ψ_{b}\right> \text{.} \end{equation} For linear polarization we can take the unit vector $\hat{\mathbf{e}}$ to be real and write the interaction Hamiltonain between states $\left|ψ_{a}\right>$ and $\left|ψ_{b}\right>$ as [[Shore(2011)]][Shore_2011] \begin{equation} \begin{aligned} H_{int,ab} &= H_{int,ab} = -d_{ψ_{a}, ψ_{b}} E(t) \cos\left(ω_{L} t - φ\right) ≡ ℏ Ω_{ab}(t) \cos\left(ω_{L}t + φ\right) \text{,} \end{aligned} \end{equation} where we denote the Rabi frequency by $Ω_{ab}(t)$. There are many definitions of Rabi frequency $Ω$ in the literature. The chosen Rabi frequency $Ω_{ab}(t)$ refers to the frequency of population oscillations of resonant two-state system, i.e., when the Rabi frequency remains constant, the populations of resonant two-state system undergo periodic Rabi oscillations at the Rabi frequency $Ω$. The Rabi frequency is defined by \begin{equation} Ω_{ab}(t) ≡ -\frac{d_{ψ_{a}, ψ_{b}} E(t)}{ℏ} \text{.} \end{equation}
In rotating reference frame off-diagonal elements of the atom-laser interaction Hamiltonian $H_{int, ab}$ [[Shore(2011)]][Shore_2011] are \begin{equation} H_{int, ab} = H_{int, ba}^{*} = -d_{ψ_{a}, ψ_{b}} E(t) \cos\left(ω_{L} t - φ\right) e^{-iω_{L}t} = \frac{1}{2} ℏ Ω_{ab}(t) \left(e^{-i φ} + e^{-2i ω_{L}t + iφ}\right) \end{equation} Rotating wave approximation (RWA) of atom-laser interaction is performed by disregarding the terms that vary as $2ω_{L}t$ (counter-rotating terms) when added to constant terms \begin{equation} \cos\left(ω_{L} t - φ\right) e^{-iω_{L}t} = \frac{1}{2} \left(e^{-i φ} + e^{-2i ω_{L}t + iφ}\right) \xrightarrow{RWA} \frac{1}{2} e^{-i φ} \end{equation} Note that the absolute value of the phase $φ$ is not controllable and only when one compares two pulses, both within a coherence time, or pulses affecting two locations, does one need to keep track of phases, therefore, phase $φ$ can be taken as zero at any convenient time \begin{equation} \frac{1}{2} e^{-i φ} → \frac{1}{2}\text{.} \end{equation} Thus the Rabi frequency can be made real $Ω_{ab}(t)$ such that $Ω_{ab}^{*}(t) = Ω_{ab}(t)$. Thus we obtain the atom-laser interaction Hamiltonain in the rotating wave approximation between state $\left|ψ_{a}\right>$ and $\left|ψ_{b}\right>$ as \begin{equation} H_{int, ab} = \frac{1}{2} ℏ Ω_{ab}(t) \text{.} \end{equation} And the atom-laser interaction Hamiltonain in the rotating wave approximation is given by \begin{equation} H_{int} = \frac{ℏ Ω_{23}(t)}{2}\left|2,n+1\right>\left< 3,n\right| + \frac{ℏ Ω_{23}^{*}(t)}{2}\left|3,n\right>\left< 2,n+1\right| + \frac{ℏ Ω_{24}(t)}{2}\left|2,n+1\right>\left< 4,n\right| + \frac{ℏ Ω_{24}^{*}(t)}{2}\left|4,n\right>\left< 2,n+1\right| \end{equation}
When simulating the dynamics of the Hamiltonian we are modeling hyperfine transitions between states $\left|F_{g}\right> → \left|F_{e}\right>$ without taking into account Zeeman sublevels, therefore, for linearly polarized excitation laser, effective reduced dipole moment $d_{eff (F_{g} → F_{e})}$ [[Steck(2010a)]][Steck_Cs_2010], [[Steck(2010b)]][Steck_Na_2010] can be expressed as \begin{equation} \begin{aligned} \left|d_{eff (F_{g} → F_{e})}\right|^{2} &= \left|\left< F_{g}\,\middle|\middle|\, \vec{d} · \hat{\mathbf{e}} \,\middle|\middle|\,F_{e}\right>\right|^{2} = S_{F_{g} F_{e}} \left|\left< J_{g}\,\middle|\middle|\, \vec{d} · \hat{\mathbf{e}}\,\middle|\middle|\,J_{e}\right>\right|^{2} \\ S_{F_{g} F_{e}} &= ( 2 F_{e} + 1)(2 J_{g} + 1 ) \begin{Bmatrix} J_{g} & J_{e} & 1\\ F_{e} & F_{g} & 1 \end{Bmatrix}^{2}\\ \sum_{F_{e}} S_{F_{g} F_{e}} &= 1 \text{,} \end{aligned} \end{equation} where $S_{F_{g} F_{e}}$ are dimensionless relative hyperfine transition-strength factors of each of the $\left|F_{g}\right> → \left|F_{e}\right>$ transitions. Numerical value of the reduced dipole matrix element $\left< J_{g}\,\middle|\middle|\, \vec{d} · \hat{\mathbf{e}}\,\middle|\middle|\,J_{e}\right>$ can be calculated using excited state lifetime from expression [[Loudon(2000)]][Loudon_2000] \begin{equation} \frac{1}{τ_{e}} = Γ_{J_{e} J_{g}} = \frac{ω_{0}^{3}}{3 π \mathit{ε}_{0} ℏ c^3} \frac{2 J_{g} + 1}{2 J_{e} + 1} \left|\left< J_{g}\,\middle|\middle|\, \vec{d}\,\middle|\middle|\,J_{e}\right>\right|^{2} \text{.} \end{equation}
It is convenient to rewrite Rabi frequency $Ω_{ab}(t)$ using reduced Rabi frequency $Ω_{red}(t)$ and relative hyperfine transition-strength factors $S_{F_{g} F_{e}}$, where reduced Raby frequency is defined by \begin{equation} Ω_{red}(t) ≡ \frac{ E(t) \left< J_{g}\,\middle|\middle|\, \vec{d}\,\middle|\middle|\,J_{e}\right>}{ℏ} \text{.} \end{equation} This allows us to obtain Rabi frequency between states $\left|ψ_{a}\right>$ and $\left|ψ_{b}\right>$ as \begin{equation} Ω_{ab}(t) = Ω_{red}(t) \sqrt{S_{F_{a} F_{b}}} \end{equation}
Critical Rabi frequency $Ω_{red,cr}$ is used to estimate at what reduced Rabi frequency value $≈ 63 \%$ of the population will be trapped in a dark state after laser-atom interaction and effects of optical pumping become pronounced in spectral line signal. \begin{equation} \begin{aligned} Ω_{red,cr} = \sqrt{\frac{4 τ_{e} ω_{43}^{2}}{τ_{tr} \Pi_{31} S_{32} \sqrt{\pi}}} \end{aligned} \end{equation}
Critical Rabi frequency $Ω_{red,cr99}$ is used to estimate at what reduced Rabi frequency value $≈ 99 \%$ of the population will be trapped in a dark state after laser-atom interaction, and leading to population depletion in time period smaller than the interaction time $τ_{tr}$. \begin{equation} \begin{aligned} Ω_{red,cr99} = \sqrt{-\ln\left( 1-\frac{99}{100}\right)} Ω_{red,cr} ≈ 2.15 Ω_{red,cr} \end{aligned} \end{equation}
Spectral line signal $J(Δ, Ω_{red})$ is given by a number of photons emitted from an atom during transit time through the excitation laser zone, therefore we can express spectral line signal as \begin{equation} J(Δ, Ω_{red}) = Γ \int\limits_{-∞}^{∞} (ρ_{4,4}(t) + ρ_{3,3}(t) )\, dt \text{.} \end{equation} By using adiabatic aproximation and optical Bloch equations from Lindblad master equation, we can obtain approximate analitic expression of spectral line signal $J(Δ, Ω_{red})$ as \begin{equation} \begin{aligned} J(Δ, Ω_{red}) &= \frac{S_{24} 4 \left(- ω_{43} - Δ \right)^{2} }{S_{23} \left( 4 \left(-Δ \right)^{2} + S_{24} Ω_{red}^{2} + Γ^{2} \right) Π_{31} } \\ &× \left(1- \exp\left( -τ_{tr} \frac{\sqrt{π} Γ Π_{31} S_{23} Ω_{red}^{2}}{8 \left(- ω_{43} - Δ - ½\sqrt{S_{24}} Ω_{red} - ½\sqrt{S_{23}} Ω_{red}\right)^{2}} \right) \right)\text{.} \end{aligned} \end{equation}
In [7]:
las_plot_jored(na_jored_expcase_list_1)
plot4lme(expcase_list[0], expcase_list[0].result)
plot4lme(expcase_list[1], expcase_list[1].result)
In [8]:
las_plot_jored(cs_jored_expcase_list_1)
plot4lme(expcase_list[2], expcase_list[2].result)
plot4lme(expcase_list[3], expcase_list[3].result)
In [9]:
las_plot_jdelta(na_jdelta_expcase_list_1)
las_plot_jdelta(na_jdelta_expcase_list_2)
las_plot_jdelta(cs_jdelta_expcase_list_1)
las_plot_jdelta(cs_jdelta_expcase_list_2)
Please note that when simulating the dynamics of the total Hamiltonian:
We set the reduced Plank constant value from QuTiP backend, i.e., it is set to \begin{equation} ℏ = 1 \end{equation}
For laser-atom interaction Hamiltonian we are using a precomputed product of all of the coefficients before time-depenent exponent function \begin{equation} hcf_{ab} = \frac{ℏ Ω_{red} \sqrt{S_{F_{a} F_{b}}} }{2} \end{equation}
Compute population dinamics for $Ω_{red} = Ω_{red,cr}$ and $Ω_{red} = Ω_{red,cr99}$.
Compute $J(Ω_{red})$ for $Ω_{red}$ from $0.25$ MHz to $50$ MHz.
Compute $J(Δ)$ for $Δ$ from $-40$ MHz to $40$ MHz.
Compare $J(Δ,Ω_{red})$ evaluated numerically from density matrix with approximate analytic expression.
In [3]:
%pylab inline
%config InlineBackend.figure_formats = {'svg',}
import seaborn as sns
flatui = ["#9b59b6", "#3498db", "#95a5a6", "#e74c3c", "#34495e", "#2ecc71"]
# rc('axes', prop_cycle=cycler('color', flatui))
# sns.palplot(sns.color_palette(flatui))
sns.set(palette=flatui)
rc('svg', fonttype='none')
from qutip import *
from fractions import Fraction
from IPython.display import Markdown
In [4]:
def solve4lme(args):
'''
Solve RWA Hamiltonian for 4-level system.
RWA Hamiltonian states |1,n+1⟩, |2,n+1⟩, |3,n⟩, |4,n⟩.
RWA Hamiltonian energies ε1=-w21, ε2=0, ε3=-w43-Δ, ε4=-Δ.
A classical laser field is in resonance between states |2⟩↔|4⟩.
Excitation schematics (bare atomic states picture)
# |4⟩
# -----------E4
# ↓Δ
# ----------- |3⟩
# / ↓Γ*Π42 ---------------E3
# Ω24 ↕ ↓ / ↓Π32*Γ ↓(1-Π32)*Γ
# / ↓ Ω23↕ ↓ ↓
# -----------------------E2 ↓
# |2⟩ -----------E1
# |1⟩
Usage::
>>> args = {'delta': 0, # detuning Δ from excited state |4⟩ in s^-1
'gamma': 1, # excited state decay rate Γ
'tint': 5920, # laser-atom interaction time, np.float32
'w21': 1756, # hyperfile splitting 1⇥2
'w43': 48, # hyperfile splitting 3⇥4
'hcf23': 3, # interaction coeff 2↔3, np.float32
'hcf24': 5, # interaction coeff 2↔4, np.float32
'cbr1': 1, # branching ratio Π42 for |2⟩→|4⟩
'cbr2': 7/12, # branching ratio Π32 for |3⟩→|2⟩
'cbr3': 5/12, # branching ratio Π31=(1-Π32) for |3⟩→|1⟩
'nsteps': 50000} # Max. number of internal ode steps/call.
>>> results = solve4lme(args)
:param args: A dictionary of the form ``(parameter, value)``.
:return: QuTiP Result object with mesolve data.
:rtype: A :class:`qutip.solver.Result`
'''
# Define atomic states |1>, |2>, |3>, |4>
st1, st2, st3, st4 = map(lambda st: basis(4, st), range(4))
# Operators for the diagonal elements of atomic Hamiltonian
# outer product of |1><1|, |2><2|, |3><3|, |4><4|
sig11, sig22, sig33, sig44 = map(ket2dm, [st1, st2, st3, st4])
# Operators for the off-diagonal elements of Hamiltonian
# Interaction Hamiltonian |2>↔|3> and |2>↔|4>
# Decay modes (collapse operators) for |4>→|2>, |3>→|2>, |3>→|1>
sig24 = st2 * st4.dag() # |2><4|
sig23 = st2 * st3.dag() # |2><3|
sig13 = st1 * st3.dag() # |1><3|
# Define time independent collapse operators
# collapse operator for 4->2, sqrt(Π42*Γ) * |2><4|
C1 = np.sqrt(args['cbr1'] * args['gamma']) * sig24
# collapse operator for 3->2, sqrt(Π32*Γ) * |2><3|
C2 = np.sqrt(args['cbr2'] * args['gamma']) * sig23
# collapse operator for 3->1, sqrt(Π31*Γ) * |1><3|
C3 = np.sqrt(args['cbr3'] * args['gamma']) * sig13
# Define list of collapse operators
c_op_list = [C1, C2, C3]
# Define time vector
t = linspace(-args['tint']*2, args['tint']*2, 101)
# Set up the time independent system Hamiltonians
# ε1=-w21, ε2 = 0, ε3=-w43-Δ, ε4=-Δ.
# state1 energy level position ε1=-w21
HS1 = -args['w21'] * sig11
# state3 energy level position ε3=-w43-Δ
HS3 = (-args['w43'] - args['delta']) * sig33
# state4 energy level position ε4=-Δ
HS4 = -args['delta'] * sig44
# Set up operators for the time varying Hamiltonians
# for laser-atom interaction Ω23(t) and Ω24(t)
HI1 = sig23.dag() + sig23
HI2 = sig24.dag() + sig24
# Set up the time varying RWA Hamiltonian with time dependant
# coefficients based on QuTip Cython string functions format
HRWA = [HS1, HS3, HS4,
[HI1, 'hcf23 * exp(-2*(t / tint) ** 2)'],
[HI2, 'hcf24 * exp(-2*(t / tint) ** 2)']]
# Define initial state as state |2>
psi0 = st2
# Define ODE solver options
opts=Odeoptions()
opts.nsteps = args['nsteps']
# Workaround for QuTip version 3.1.0 error with Cython string functions
# Convert Cython string function variable values to numpy.float32
# # error: two or more data types in declaration specifiers
# # typedef npy_double _Complex __pyx_t_npy_double_complex;
args['hcf23'] = np.float32(args['hcf23'])
args['hcf24'] = np.float32(args['hcf24'])
# Solve RWA Hamiltonian
output = mesolve(HRWA, psi0, t, c_op_list, [sig11, sig22, sig33, sig44],
args=args, options=opts, progress_bar=True)
# return mesolve
return output
def plot4lme(atomdata, results=None, saveplot=None):
'''Plot QuTiP Result object.'''
if results != None:
rs = results
else:
rs = qload(atomdata.filename)
rho11, rho22, rho33, rho44 = rs.expect
timescale = 1e6
t = rs.times * timescale
# Define pump strength as a function of time for plotting
wp = lambda t, tw, A: A * np.exp(-2*(t / tw) ** 2)
# Plot the results
fig = figure()
subplot(211)
if atomdata != None:
# colors #a6cee3 #2078b4 #afdd8a #35a12e #fa9897 #e31a1c
detuning = atomdata.delta
detuning_MHz = detuning / (timescale * 2*pi)
detuning_MHz_str = str(round(detuning_MHz,2))
fig.suptitle(atomdata.symbol + ', $τ_{int} = %s$ (s), ' % str(atomdata.tint) \
+ '$Δ = %s$ (MHz)' % detuning_MHz_str)
ored_sat_24_MHz = atomdata.osat / (timescale * 2 * pi * np.sqrt(atomdata.st24))
ored_sat_24_MHz_str = str(round(ored_sat_24_MHz,2))
axhline(y=ored_sat_24_MHz, color='#a6cee3',
label='$Ω_{red,sat,24} ≈ %s$ MHz' % ored_sat_24_MHz_str)
ored_cr_MHz = atomdata.__class__(atomdata.tint).ocr \
/ (timescale * 2 * pi)
ored_cr_MHz_str = str(round(ored_cr_MHz,2))
axhline(y=ored_cr_MHz, color='#afdd8a',
label='$Ω_{red,cr} ≈ %s$ MHz' % ored_cr_MHz_str)
ored_cr_pop_99_MHz = atomdata.__class__(atomdata.tint, pop=99).ocr \
/ (timescale * 2 * pi)
ored_cr_pop_99_MHz_str = str(round(ored_cr_pop_99_MHz,2))
axhline(y=ored_cr_pop_99_MHz, color='#fa9897',
label='$Ω_{red,cr}|_{pop=99} ≈ %s$ MHz' % ored_cr_pop_99_MHz_str)
plot(t, wp(t, atomdata.tint * timescale,
atomdata.hcf23 * 2 / (timescale * 2*pi)), '-', label='$Ω_{23}$')
plot(t, wp(t, atomdata.tint * timescale,
atomdata.hcf24 * 2 / (timescale * 2*pi)), '-', label='$Ω_{24}$')
ylabel('Coupling strength $Ω$ (MHz)')
lg_1 = legend()
lg_1.draw_frame(True)
subplot(212)
plot(t, rho11, '-', label='$ρ_{1,1}$')
plot(t, rho22, '-', label='$ρ_{2,2}$')
plot(t, rho33, '-', label='$ρ_{3,3}$')
plot(t, rho44, '-', label='$ρ_{4,4}$')
ylabel('Population $ρ$ (rel. units)')
xlabel('Time ($\mu s$)')
lg_2 = legend()
lg_2.draw_frame(True)
# check if filename for exporting figure is provided
if saveplot != None:
savefig(saveplot)
show()
class AtomData4Levels:
'''Setup parameters for 4-level system RWA Hamiltonian
Usage::
>>> class Atom4levelD2line(AtomData4Levels):
gamma = 1 # excited state decay rate Γ in s^-1
taue = 1 # exctied state lifetime in s
w21 = 1756 # hyperfile splitting 1⇤⇥2 in s^-1
w43 = 48 # hyperfile splitting 3⇤⇥4 in s^-1
cbr1 = 1 # branching ratio Π42 for |2⟩→|4⟩
cbr2 = 1/2 # branching ratio Π32 for |3⟩→|2⟩
cbr3 = 1/2 # branching ratio Π31=(1-Π32) for |3⟩→|1⟩
st23 = 1 # rel. HF trans. strength factor
st24 = 1 # rel. HF trans. strength factor
>>> tint = 5920 # laser-atom interaction time in s
>>> csd2 = Atom4levelD2line(tint, delta)
:param tint : laser-atom interaction time in s
:param delta: detuning Δ from excited state |4⟩ in s^-1
:param gamma: excited state decay rate Γ in s^-1
:param taue : exctied state lifetime in s
:param w21 : hyperfile splitting 1⇤⇥2 in s^-1
:param w43 : hyperfile splitting 3⇤⇥4 in s^-1
:param cbr1 : branching ratio Π42 for |2⟩→|4⟩
:param cbr2 : branching ratio Π32 for |3⟩→|2⟩
:param cbr3 : branching ratio Π31=(1-Π32) for |3⟩→|1⟩
:param st23 : HF transition strength factors divided by
:param st24 : square of reduced dipole matrix element
:return: Parameters for 4-level system RWA Hamiltonian.
:rtype: A :class:`AtomData4Levels`
'''
listargs = ('delta', 'tint', 'gamma', 'taue', 'w21', 'w43', 'hcf23', 'hcf24',
'cbr1', 'cbr2', 'cbr3', 'nsteps')
def __init__(self, tint, delta=0, ored=None, pop=None, nsteps=50000):
self.tint = tint
self.delta = delta
self.ored = ored
self.pop = pop
self.nsteps = nsteps
self.update()
def update(self):
'''Update parameters'''
self.osat = self.omega_saturation(self.gamma)
if self.pop == None:
self.pop = (1-np.exp(-1))*100
self.pop_coef = 1
else:
self.pop_coef = - np.log(1-self.pop/100)
self.ocr = self.omega_critical(self.tint, self.taue, self.w43-self.delta,
self.st23, self.cbr3, self.pop_coef)
if self.ored == None:
self.ored = self.ocr
self.hcf23 = self.hcf(self.ored, self.st23)
self.hcf24 = self.hcf(self.ored, self.st24)
self.argsaslist = (self.delta, self.tint, self.gamma, self.taue, self.w21,
self.w43, self.hcf23, self.hcf24, self.cbr1, self.cbr2,
self.cbr3, self.nsteps)
self.args = self.gendict(self.listargs, self.argsaslist)
self.filename = self.args_to_filename(self.listargs, **self.args)
def args_to_filename(self, listargs, **kwargs):
'''Return filename from list of args'''
return ''.join(list(map(lambda i: i + ':' + str(kwargs[i]), listargs)))
def omega_saturation(self, gamma):
'''Return aturation Rabi frequency value for reduced Rabi frequency'''
return gamma / np.sqrt(2)
def omega_critical(self, ti, te, w43, s23, p3, pop_coef):
'''Return critical Rabi frequency value for reduced Rabi frequency'''
return np.sqrt(4 * pop_coef * te * (w43 ** 2) / (ti * p3 * s23))
def hcf(self, ored, strength):
'''Return interaction Hamiltonian coeff'''
return ored * np.sqrt(strength) / 2
def gendict(self, listargs, args):
'''Return Cesium D2 Line Data as dictionary'''
return dict(zip(listargs, args))
class CsD2data:
'''Store cesium D2 line data.
Data obtained from Daniel A. Steck, Cesium D Line Data (2010).
Available online at http://steck.us/alkalidata
'''
symbol = 'Cs'
nspin = Fraction(7,2) # nuclear spin
jg = Fraction(1,2) # ground state J
je = Fraction(3,2) # excited state J
# relative hyperfine transition strength factors [D. Steck, Alkali D Line Data]
# S_{F_{g} F_{e}}
# F_{g} = 3
s32, s33, s34 = Fraction(20,56), Fraction(21,56), Fraction(15,56)
# F_{g} = 4
s43, s44, s45 = Fraction(7,72), Fraction(21,72), Fraction(44,72)
skeys = ((3,2), (3,3), (3,4), (4,3), (4,4), (4,5))
svals = (s32, s33, s34, s43, s44, s45)
sge = dict(zip(skeys, svals))
gamma = 32.889e+6 # decay rate Γ 32.889(84)e+6 in s^-1
taue = 30.405e-9 # lifetime 30.405(77) in s
wg43 = 2*pi * 9192.631770e+6 # F_{g} hfs 3⇤⇥4 9192.631770e+6 in s^-1
we54 = 2*pi * 251.0916e+6 # F_{e} hfs 4⇤⇥5 251.0916(20)e+6 in s^-1
we43 = 2*pi * 201.2871e+6 # F_{e} hfs 3⇤⇥4 201.2871(11)e+6 in s^-1
we32 = 2*pi * 151.2247e+6 # F_{e} hfs 2⇤⇥3 151.2247(16)e+6 in s^-1
cbr54 = 1 # branching ratio Π54 for |F_{5}⟩→|F_{4}⟩
cbr43 = Fraction(5,12) # branching ratio Π43 for |F_{4}⟩→|F_{3}⟩
cbr44 = Fraction(7,12) # branching ratio Π44 for |F_{4}⟩→|F_{4}⟩
cbr33 = Fraction(3,4) # branching ratio Π33 for |F_{3}⟩→|F_{3}⟩
cbr34 = Fraction(1,4) # branching ratio Π34 for |F_{3}⟩→|F_{4}⟩
cbr23 = 1 # branching ratio Π23 for |F_{2}⟩→|F_{3}⟩
class Cs4L(AtomData4Levels,CsD2data):
'''Store cesium D2 line data parameters for 4-level system RWA Hamiltonian'''
w21 = np.float32(CsD2data.wg43)
w43 = np.float32(CsD2data.we54)
cbr1 = np.float32(CsD2data.cbr54)
cbr2 = np.float32(CsD2data.cbr44)
cbr3 = np.float32(CsD2data.cbr43)
st24 = np.float32(CsD2data.sge[4,5])
st23 = np.float32(CsD2data.sge[4,4])
class NaD2data:
'''Store sodium D2 line data
Data obtained from Daniel A. Steck, Sodium D Line Data (2010).
Available online at http://steck.us/alkalidata
'''
symbol = 'Na'
nspin = Fraction(3,2) # nuclear spin
jg = Fraction(1,2) # ground state J
je = Fraction(3,2) # excited state J
# relative hyperfine transition strength factors [D. Steck, Alkali D Line Data]
# S_{F_{g} F_{e}}
# F_{g} = 1
s10 = Fraction(2,12)
s11 = Fraction(5,12)
s12 = Fraction(5,12)
# F_{g} = 2
s21 = Fraction(1,20)
s22 = Fraction(5,20)
s23 = Fraction(14,20)
skeys = ((1,0), (1,1), (1,2), (2,1), (2,2), (2,3))
svals = (s10, s11, s12, s21, s22, s23)
sge = dict(zip(skeys, svals))
gamma = 61.542e+6 # decay rate Γ 61.542(29)e+6 in s^-1
taue = 16.2492e-9 # lifetime 16.2492(77) in s
wg21 = 2*pi * 1771.6261288e+6 # F_{g} hfs 1⇤⇥2 1771.6261288(10)e+6 in s^-1
we32 = 2*pi * 58.326e+6 # F_{e} hfs 2⇤⇥3 58.326(43)e+6 in s^-1
we21 = 2*pi * 34.344e+6 # F_{e} hfs 1⇤⇥2 34.344(49)e+6 in s^-1
we10 = 2*pi * 15.810e+6 # F_{e} hfs 0⇤⇥1 15.810(80)e+6 in s^-1
cbr32 = 1 # branching ratio Π32 for |F_{3}⟩→|F_{2}⟩
cbr21 = Fraction(1,2) # branching ratio Π21 for |F_{2}⟩→|F_{1}⟩
cbr22 = Fraction(1,2) # branching ratio Π22 for |F_{2}⟩→|F_{2}⟩
cbr11 = Fraction(5,6) # branching ratio Π11 for |F_{1}⟩→|F_{1}⟩
cbr12 = Fraction(1,6) # branching ratio Π12 for |F_{1}⟩→|F_{2}⟩
cbr01 = 1 # branching ratio Π01 for |F_{0}⟩→|F_{1}⟩
class Na4L(AtomData4Levels,NaD2data):
'''Store sodium D2 line data parameters for 4-level system RWA Hamiltonian'''
w21 = np.float32(NaD2data.wg21);
w43 = np.float32(NaD2data.we32)
cbr1 = np.float32(NaD2data.cbr32)
cbr2 = np.float32(NaD2data.cbr22)
cbr3 = np.float32(NaD2data.cbr21)
st24 = np.float32(NaD2data.sge[2,3])
st23 = np.float32(NaD2data.sge[2,2])
def lorentz_norm(delta, deltapeak, gamma):
'''Return Lorentzian profile of natural linewidth norlalized to 1 at peak.'''
g, w, w0 = gamma, delta, deltapeak
return g ** 2 / (4 * ((w - w0)**2 + (g ** 2)/4))
def las_analytic(delta, omegared, atom):
'''Return laser absorption signal J using approximate analytic expression.'''
delta23 = - atom.w43 - delta
delta24 = - delta
delta23eff = delta23 - omegared * (np.sqrt(atom.st23) + np.sqrt(atom.st24)) / 2
delta24effsq = (delta24 ** 2) + (atom.st24 * omegared ** 2)/4
tau_pump = (atom.taue * 8 * (delta23eff ** 2)) \
/ (np.sqrt(pi) * atom.cbr3 * atom.st23 * omegared ** 2)
J = (atom.st24 * 4 * (delta23 ** 2) * (1 - np.exp(-atom.tint/tau_pump))) \
/ (atom.st23 * (4 * (delta24effsq) + (atom.gamma) ** 2) * atom.cbr3)
return J
def las_numeric_single(results, gamma):
'''Return laser absorption signal J by integrating populations
rho33, rho44 from QuTiP result object using Trapezoidal rule
and multiplying it by decay rate Γ.'''
rho11, rho22, rho33, rho44 = results.expect
times = results.times
return (numpy.trapz(rho33, times) + numpy.trapz(rho44, times)) * gamma
def solve4lme_list(atomlist, precomputed=True):
'''Compute and add RWA Hamiltonian QuTiP result object
to AtomData4Levels object in atomlist.'''
for atom in atomlist:
if precomputed:
atom.result = qload(atom.filename)
else:
atom.result = solve4lme(atom.args)
qsave(atom.result, atom.filename)
atom.jnumval = las_numeric_single(atom.result, atom.gamma)
def las_plot_jored(jored_expcase_list):
'''Plot J(Ω_{red}) for fixed Δ'''
# colors #a6cee3 #2078b4 #afdd8a #35a12e #fa9897 #e31a1c
atom = jored_expcase_list[0]
delta = atom.delta
delta_MHz = delta * AFtoMHz
delta_MHz_str = str(round(delta_MHz,2))
jored_list = np.array(list(map(lambda atom: atom.ored, jored_expcase_list)))
jored_list_vals = np.array(list(map(lambda atom: atom.jnumval, jored_expcase_list)))
jored_num_list = np.linspace(jored_list[0], jored_list[-1], 201)
ored_sat_24 = atom.osat / np.sqrt(atom.st24)
ored_sat_24_MHz = ored_sat_24 * AFtoMHz
ored_sat_24_MHz_str = str(round(ored_sat_24_MHz,2))
orec_cr = atom.__class__(atom.tint).ocr
ored_cr_MHz = orec_cr * AFtoMHz
ored_cr_MHz_str = str(round(float(ored_cr_MHz),2))
orec_cr_pop_99 = atom.__class__(atom.tint,pop=99).ocr
ored_cr_pop_99_MHz = orec_cr_pop_99 * AFtoMHz
ored_cr_pop_99_MHz_str = str(round(float(ored_cr_pop_99_MHz),2))
figure()
axvline(x=ored_sat_24_MHz, color='#a6cee3',
label='$Ω_{red,sat,24} ≈$ %s MHz' % ored_sat_24_MHz_str)
axvline(x=ored_cr_MHz, color='#afdd8a',
label='$Ω_{red,cr} ≈$ %s MHz' % ored_cr_MHz_str)
axvline(x=ored_cr_pop_99_MHz, color='#fa9897',
label='$Ω_{red,cr}|_{pop=0.01} ≈$ %s MHz' % ored_cr_pop_99_MHz_str)
plot(jored_list * AFtoMHz,jored_list_vals, '-', label='numeric calc.')
plot(jored_num_list * AFtoMHz, las_analytic(delta, jored_num_list, atom),
'--', label='analytic approx.')
ylabel('$J(Ω_{red})$ (arb. units)')
xlabel('$Ω_{red}$ (MHz)')
title('%s, $Δ$ = %s (MHz), $τ_{int}$ = %s (s)' \
% (atom.symbol, delta_MHz_str, str(atom.tint)))
lg = legend()
lg.draw_frame(True)
show()
def las_plot_jdelta(jdelta_expcase_list):
'''Plot J(Δ) for fixed Ω_{red}'''
# colors #a6cee3 #2078b4 #afdd8a #35a12e #fa9897 #e31a1c
atom = jdelta_expcase_list[0]
ored = atom.ored
ored_MHz = ored * AFtoMHz
ored_MHz_str = str(round(ored_MHz,2))
gamma = atom.gamma
gamma_MHz = gamma * AFtoMHz
gamma_MHz_str = str(round(gamma_MHz,2))
jdelta_list = np.array(list(map(lambda atom: atom.delta, jdelta_expcase_list)))
jdelta_list_vals = np.array(list(map(lambda atom: atom.jnumval, jdelta_expcase_list)))
jdelta_num_list = np.linspace(jdelta_list[0], jdelta_list[-1], 201)
figure()
axvline(x=-gamma_MHz/2, color='#afdd8a',
label='±$Γ/2$, $Γ≈$ %s MHz (FWHM)' % gamma_MHz_str)
axvline(x=gamma_MHz/2, color='#afdd8a')
plot(jdelta_num_list * AFtoMHz,
lorentz_norm(jdelta_num_list, 0, gamma) * jdelta_list_vals.max(),
'-', label='natural linewidth')
plot(jdelta_list * AFtoMHz, jdelta_list_vals, '-', label='numeric calc.')
plot(jdelta_num_list * AFtoMHz, las_analytic(jdelta_num_list, ored, atom),
'--', label='analytic approx.')
ylabel('$J(Δ)$ (arb. units)')
xlabel('$Δ$ (MHz)')
title('%s, $Ω_{red}$ = %s (MHz), $τ_{int}$ = %s (s)' \
% (atom.symbol, ored_MHz_str, str(atom.tint)))
lg = legend()
lg.draw_frame(True)
show()
# converting from/to frequency in MHz to/from angular frequency s^{-1}
MHztoAF = 2 * pi * 1e+6
AFtoMHz = 1/(MHztoAF)
In [5]:
%%capture
expcase_list = Na4L(0.0002), Na4L(0.0002, pop=99), Cs4L(0.0002), Cs4L(0.0002, pop=99)
# set precomputed to False to avoid loading precomputed and saved qutip results
solve4lme_list(expcase_list, precomputed=True)
In [6]:
%%capture
# Calculate numerical laser absorption signal J(Ω_{red})|Δ=0
# laser-atom interaction time 0.0002 s,
# Ω_{red} from 0.5 MHz to 50 MHz with vairable step size
na_jored_list = np.linspace(0.25,2.25, 9)
na_jored_list = np.concatenate((na_jored_list, np.linspace(2.5,14.5, 25)), axis=0)
na_jored_list = np.concatenate((na_jored_list, np.linspace(15,50, 15)), axis=0)
na_jored_MHz_list_1 = na_jored_list
na_jored_expcase_list_1 = list(map(lambda omegaMHz:
Na4L(0.0002, ored = omegaMHz * MHztoAF),
na_jored_MHz_list_1))
# set precomputed to False to avoid loading precomputed and saved qutip results
solve4lme_list(na_jored_expcase_list_1, precomputed=True)
# Calculate numerical laser absorption signal J(Ω_{red})|Δ=0
# laser-atom interaction time 0.0002 s,
# Ω_{red} from 0.5 MHz to 50 MHz with vairable step size
cs_jored_list = np.linspace(0.25,2.25, 9)
cs_jored_list = np.concatenate((cs_jored_list, np.linspace(2.5,14.5, 13)), axis=0)
cs_jored_list = np.concatenate((cs_jored_list, np.linspace(15,50, 15)), axis=0)
cs_jored_MHz_list_1 = cs_jored_list
cs_jored_expcase_list_1 = list(map(lambda omegaMHz:
Cs4L(0.0002, ored = omegaMHz * MHztoAF),
cs_jored_MHz_list_1))
# set precomputed to False to avoid loading precomputed and saved qutip results
solve4lme_list(cs_jored_expcase_list_1, precomputed=True)
# Calculate numerical laser absorption signal J(Δ)|Ω_{red}=Ω_{red,cr}|Δ=0
# laser-atom interaction time 0.0002 s,
# Δ from -40 MHz to 40 MHz with vairable step size
na_jdelta_list = np.linspace(-40,-5, 15)
na_jdelta_list = np.concatenate((na_jdelta_list, np.linspace(-4,-1, 4)), axis=0)
na_jdelta_list = np.concatenate((na_jdelta_list, np.linspace(-0.5,0.5, 5)), axis=0)
na_jdelta_list = np.concatenate((na_jdelta_list, np.linspace(1,4, 4)), axis=0)
na_jdelta_list = np.concatenate((na_jdelta_list, np.linspace(5,40, 15)), axis=0)
na_jdelta_MHz_list_1 = na_jdelta_list
na_jdelta_expcase_list_1 = list(map(lambda deltaMHz:
Na4L(0.0002, delta = deltaMHz * MHztoAF,
ored = Na4L(0.0002).ocr),
na_jdelta_MHz_list_1))
# set precomputed to False to avoid loading precomputed and saved qutip results
solve4lme_list(na_jdelta_expcase_list_1, precomputed=True)
# Calculate numerical laser absorption signal J(Δ)|Ω_{red}=Ω_{red,cr}|Δ=0
# laser-atom interaction time 0.0002 s,
# Δ from -40 MHz to 40 MHz with vairable step size
na_jdelta_list_2 = np.linspace(-40,-5, 15)
na_jdelta_list_2 = np.concatenate((na_jdelta_list_2, np.linspace(-4,-1, 4)), axis=0)
na_jdelta_list_2 = np.concatenate((na_jdelta_list_2, np.linspace(-0.5,0.5, 5)), axis=0)
na_jdelta_list_2 = np.concatenate((na_jdelta_list_2, np.linspace(1,4, 4)), axis=0)
na_jdelta_list_2 = np.concatenate((na_jdelta_list_2, np.linspace(5,40, 15)), axis=0)
na_jdelta_MHz_list_2 = na_jdelta_list_2
na_jdelta_expcase_list_2 = list(map(lambda deltaMHz:
Na4L(0.0002, delta = deltaMHz * MHztoAF,
ored = Na4L(0.0002, pop=99).ocr),
na_jdelta_MHz_list_2))
# set precomputed to False to avoid loading precomputed and saved qutip results
solve4lme_list(na_jdelta_expcase_list_2, precomputed=True)
# Calculate numerical laser absorption signal J(Δ)|Ω_{red}=Ω_{red,cr}|Δ=0
# laser-atom interaction time 0.0002 s,
# Δ from -40 MHz to 40 MHz with vairable step size
cs_jdelta_list = np.linspace(-40,-5, 15)
cs_jdelta_list = np.concatenate((cs_jdelta_list, np.linspace(-4,-1, 4)), axis=0)
cs_jdelta_list = np.concatenate((cs_jdelta_list, np.linspace(-0.5,0.5, 5)), axis=0)
cs_jdelta_list = np.concatenate((cs_jdelta_list, np.linspace(1,4, 4)), axis=0)
cs_jdelta_list = np.concatenate((cs_jdelta_list, np.linspace(5,40, 15)), axis=0)
cs_jdelta_MHz_list_1 = cs_jdelta_list
cs_jdelta_expcase_list_1 = list(map(lambda deltaMHz:
Cs4L(0.0002, delta = deltaMHz * MHztoAF,
ored = Cs4L(0.0002).ocr,
nsteps=100000),
cs_jdelta_MHz_list_1))
# set precomputed to False to avoid loading precomputed and saved qutip results
solve4lme_list(cs_jdelta_expcase_list_1, precomputed=True)
# Calculate numerical laser absorption signal J(Δ)|Ω_{red}=Ω_{red,cr}|Δ=0
# laser-atom interaction time 0.0002 s,
# Δ from -40 MHz to 40 MHz with vairable step size
cs_jdelta_list_2 = np.linspace(-40,-5, 15)
cs_jdelta_list_2 = np.concatenate((cs_jdelta_list_2, np.linspace(-4,-1, 4)), axis=0)
cs_jdelta_list_2 = np.concatenate((cs_jdelta_list_2, np.linspace(-0.5,0.5, 5)), axis=0)
cs_jdelta_list_2 = np.concatenate((cs_jdelta_list_2, np.linspace(1,4, 4)), axis=0)
cs_jdelta_list_2 = np.concatenate((cs_jdelta_list_2, np.linspace(5,40, 15)), axis=0)
cs_jdelta_MHz_list_2 = na_jdelta_list_2
cs_jdelta_expcase_list_2 = list(map(lambda deltaMHz:
Cs4L(0.0002, delta = deltaMHz * MHztoAF,
ored = Cs4L(0.0002, pop=99).ocr,
nsteps=100000),
cs_jdelta_MHz_list_2))
# set precomputed to False to avoid loading precomputed and saved qutip results
solve4lme_list(cs_jdelta_expcase_list_2, precomputed=True)