An Introduction to Rydberg atoms with ARC

Authors: Nikola Šibalić, Jonathan D. Pritchard, Charles S. Adams, Kevin J. Weatherill

Updated 26/11/2018

In this notebook:

  1. Configure Notebook
  2. Rydberg atom energy levels
  3. Rydberg atom wavefunctions
  4. Matrix elements for atom interaction with electromagnetic radiation
  5. Rydberg lifetimes and black-body radiation induced transitions
  6. Creating highly excited states
  7. Rydberg atom energy level modifications in DC electric fields (Stark shifts)
  8. Interaction between the two atoms excited to Rydberg states

  9. Tuning the interactions with external fields (Förster resonance)

  10. Atomic vapour properties
  11. Advanced Rydberg examples
  12. Advanced use of ARC package: interfacing and expansions

Further reading:

For concise general overview of the current research directions in the Rydberg physics, and explanation of elementary underlying ideas, check new eBook below. It was made specifically for final year undergraduates and a new PhD reserchers and uses new electronic format to bring interactive figures for exploring different regimes and building intuition about processes (made with the help of ARC package):

The eBook above is interactive in web version and EPUB version if used with EPUB3 supporting device (e.g. iBooks on Mac/iPad, Gitden Reader in Android, or Document viewer Xreader on Linux). Under Windows we recommend using web version since the best currently available reader in Windows (Calibre) doesn't have full support for EPUB3.

Configure Notebook

This notebook uses ARC package (Alkali Rydberg Calculator) that can be downloaded from GitHub. Once you download the package to your computer set rootDir for the package below, and run the first cell, before continuing with other examples.

In [3]:
# Configure the matplotlib graphics library and configure it to show 

# show figures inline in the notebook
%matplotlib inline               
import matplotlib.pyplot as plt  # Import library for direct plotting functions
import numpy as np               # Import Numerical Python
from IPython.core.display import display, HTML #Import HTML for formatting output

# NOTE: Uncomment following lines ONLY if you are not using installation via pip
# import sys, os
# rootDir = '/path/to/arc/directory' # e.g. '/Users/Username/Desktop/ARC-Alkali-Rydberg-Calculator'
# sys.path.insert(0,rootDir)

from arc import *                 #Import ARC (Alkali Rydberg Calculator)

Rydberg Atom Energy Levels

Rydberg states are highly excited states of the outer valence electron where the properties can be scaled in terms of the principal quantum number, $n$. Originally observed in the spectral lines of hydrogen, the binding energy of the Rydberg series are given by

\begin{equation} E_{n\ell j} = - \frac{\mathrm{Ry}}{(n-\delta_{n\ell j})^2}, \end{equation}

where $n,\ell,j$ are the quantum numbers, $\mathrm{Ry}$ is the Rydberg constant and $\delta_{n\ell j}$ is the quantum defect. This defect describes the increase in binding energy for an alkali atom with respect to Hydrogen due to penetration and polarisation of the closed inner electron shells, which is most significant for the $\ell=0$ states with highly elliptical orbits, and can be neglected for states with $\ell>3$. The quantum defects are parameterised via

\begin{equation} \delta_{n\ell j} = \delta_0 + \frac{\delta_2}{(n-\delta_0)^2}+\frac{\delta_4}{(n-\delta_0)^4}+ \ldots, \end{equation}

with coefficients $\delta_{0,2,...}$ taken from measured energy levels.

The mass-corrected Rydberg constant is given by $\mathrm{Ry}=\mathrm{Ry}_\infty\times M/(M+m_e)$ where $m_e$ is the electron mass and $M$ is the atomic mass of the nucleus, with

\begin{equation} \mathrm{Ry}_\infty=\frac{e^4m_e}{16\pi^2\epsilon_0^2\hbar^2}. \end{equation}

To demonstrate the effect of quantum defects, the graph below shows the energy levels of Cesium highlighting the effect of the quantum defect between the different $\ell$ manifolds.

In [2]:
#Load parameters for Caesium

nmin=6  #Minimum n
nmax=60 #Maximum n
lmin=0  #Minimum l
lmax=3  #Maxmium l

#Plot Energy Levels of Cesium
levels = LevelPlot(atom)
# plot is interactive when called outside the IPython notebook (e.g. from Python program)
In [3]:
#Plot Quantum Defects of Cs
fig, axes = plt.subplots(1, 1, figsize=(7,5))

axes.set_ylabel('Quantum Defect $\delta_{n\ell j}$')
axes.set_title('Cs Quantum Defects')
Text(0.5,1,'Cs Quantum Defects')

Rydberg atom wavefunctions

Rydberg atom wavefunctions can be obtained from Schrödinger's equation,

\begin{equation} \left[-\frac{1}{2\mu}\left(\frac{\mathrm{d}^2}{\mathrm{d} r^2} +\frac{2}{r}\frac{\mathrm{d}}{\mathrm{d} r}\right) + \frac{\ell(\ell+1)}{2\mu r^2} + V(r)\right]R(r)= E_{n\ell j}R(r), \end{equation}

with reduced mass $\mu=Mm_e/(M+m_e)$ and $V(r)$ a model potential given by

\begin{equation} V(r) = -\frac{Z_{n\ell}(r)}{r}-\frac{\alpha_c}{2r^4}(1-e^{-(r/r_c)^6})+V_{\rm so}, \end{equation}

which at long range gives a $-1/r$ Colomb potential and at short range accounts for the finite size of the core, with radial charge given by

\begin{equation} Z_{n\ell}(r) = 1+(Z-1)e^{-a_1r}-r(a_3+a_4r)e^{-a_2r}. \end{equation}

Wavefunctions are calculated by numerical intergration of the model potential using parameters $a_{1..4},r_c$ and $\alpha_c$ taken from Marinescu Potential includes spin-orbit interaction $V_{\rm so}(r) = \alpha ~\mathbf{L}\cdot\mathbf{S}/(2r^3)$, where $\alpha$ is fine structure constant.

Rydberg atoms are large, with an average size $\langle r \rangle\propto n^2ea_0$, which leads to extremely large dipole matrix elements making them ideal for exploitng strong long range interactions. To illustrate this point, the probabiltiy distribution $\vert rR(r)\vert^2$ for $nS_{1/2}$ atomic wavefunctions are plotted below for increasing $n$ which highlights the dramatic scaling with atoms approaching $1~\mu$m in size for $n=100$ compared to $1~\unicode{x212B}$ for the groundstate!

In [4]:
pqn = [15,30,60] # principal quantum numbers of the states
colors = ["b","r","g"]
l = 0          # S state
j = 0.5        # J = 1/2

plotLegend = []
for i in range(len(pqn)):
    n = pqn[i]
    step = 0.001    
    a1,b1 = atom.radialWavefunction(l,0.5,j,\
                                       atom.getEnergy(n, l, j)/27.211,\
                                       2.0*n*(n+15.0), step)
    legendInfo, = plt.plot(a1,(b1)*(b1),\
                           "-",lw=2,color = colors[i], label = ("n = %d" % n) )
plt.xlabel(r"Distance from nucleus $r$ ($a_0$)")
plt.ylabel(r"$\vert rR(r)\vert^2$")

Matrix Elements

Dipole Matrix Elements

Relevant properties of the Rydberg states can be derived through evaluation of dipole matrix elements. Recalling the separability of the atomic wavefunctions into radial and spherical components $\psi(r,\theta,\phi)=Y_{\ell,m_\ell}(\theta,\phi)R_{n\ell}(r)$ in the uncoupled basis, the dipole matrix elements can be expressed as

\begin{equation} \langle n,\ell,m_\ell\vert r_q \vert n',\ell',m_\ell' \rangle = (-1)^{\ell-m_\ell}\begin{pmatrix}\ell&1&\ell'\\-m_\ell&q&m_\ell'\end{pmatrix}\langle\ell \vert\vert r\vert\vert\ell'\rangle, \end{equation}

where $q=-1,0,1$ corresponds to $\sigma^+,\pi$ and $\sigma^-$ transitions respectively with the reduced matrix element

\begin{equation} \langle\ell\vert\vert r\vert\vert\ell'\rangle = (-1)^\ell\sqrt{(2\ell+1)(2\ell'+1)}\begin{pmatrix}\ell&1&\ell\\0&0&0\end{pmatrix}\mathcal{R}_{n\ell\rightarrow n'\ell'}, \end{equation}

where round braces denote Wigner-$3j$ symbols and the radial matrix element is evaluated from

\begin{equation} \mathcal{R}_{n\ell\rightarrow n'\ell'} = \int_{r_{\rm i}}^{r_{\rm o}} R_{n,\ell}(r)rR_{n',\ell'}(r)r^2\mathrm{d} r. \end{equation}

For highly excited states the hyperfine-structure splitting becomes negligible, however the fine-structure splitting $\hat{\mathcal{H}}=\beta \bf\hat{L}\cdot\hat{S}$ means the relevant basis is then fine structure basis ($j,m_j$) with matrix elements

\begin{equation} \langle n,\ell,j,m_j\vert r_q \vert n',\ell',j',m_j' \rangle = (-1)^{j-m_j}\begin{pmatrix}j&1&j'\\-m_j&q&m_j'\end{pmatrix}\langle j \vert\vert r\vert\vert j'\rangle, \end{equation}

and the reduced matrix element equal to

\begin{equation} \langle j \vert\vert r\vert\vert j'\rangle = (-1)^{\ell+s+j'+1}\delta_{s,s'}\sqrt{(2j+1)(2j'+1)} \begin{Bmatrix}j&1&j'\\\ell'&s&\ell\end{Bmatrix}\langle \ell\vert\vert r\vert\vert\ell'\rangle, \end{equation}

where the curly braces denote a Wigner-$6j$ symbol.

The cell below demonstrates how to extract the relevent matrix elements for different transitions and highlights the extremely large matrix elements of the Rydberg states for transitions to neighbouring states, which scale $\propto n^2$ due to the large electron radius shown above. Thus whilst the alkali $D$ line transitions have elements of order $~e a_0$, the Rydberg states can have $>1000~e a_0$ making them ideally suited for exploiting strong, long range interactions as will be explored below.

Quadrupole Matrix Elements

For accurate interaction potentials at short range, it is necessary to also consider quadrupole matrix elements for the Rydberg atoms. getQuadrupoleMatrixElement

\begin{equation} \mathcal{R}^Q_{n\ell\rightarrow n'\ell'} = \int_{r_{\rm i}}^{r_{\rm o}} R_{n',\ell}(r)r^2R_{n,\ell'}(r)r^2\mathrm{d} r. \end{equation}

In [5]:
#Cs D2 Transition 6S_{1/2}-->6P_{3/2}
print("Cs D2 Transition 6S_{1/2}-->6P_{3/2}")
#Radial Matrix element R_{nlj\rightarrown'l'j'}
print("R_{nlj-n'l'j'} = %.3f ea_0" % atom.getRadialMatrixElement(n1,l1,j1,n2,l2,j2))
#Reduced Matrix Element <l||er||l'>
print("<l||er||l'> = = %.3f ea_0" % atom.getReducedMatrixElementL(n1,l1,j1,n2,l2,j2))
#Reduced Matrix Element <j||er||j'>
print("<j||er||j'> = %.3f ea_0" % atom.getReducedMatrixElementJ(n1,l1,j1,n2,l2,j2))
#Angular Coupling
print("<nljmj|er|n'l'j'mj'> = %.3f ea_0\n" %\

#Cs Rydberg Transition 60S_{1/2}-->60P_{3/2}

print("Cs Rydberg Transition 60S_{1/2}-->60P_{3/2}")
#Radial Matrix element R_{nlj\rightarrown'l'j'}
print("R_{nlj-n'l'j'} = %.3f ea_0" % atom.getRadialMatrixElement(n1,l1,j1,n2,l2,j2))
#Reduced Matrix Element <l||er||l'>
print("<l||er||l'> = = %.3f ea_0" % atom.getReducedMatrixElementL(n1,l1,j1,n2,l2,j2))
#Reduced Matrix Element <j||er||j'>
print("<j||er||j'> = %.3f ea_0" % atom.getReducedMatrixElementJ(n1,l1,j1,n2,l2,j2))
#Angular Coupling
print("<nljmj|er|n'l'j'mj'> = %.3f ea_0\n" % \
Cs D2 Transition 6S_{1/2}-->6P_{3/2}
R_{nlj-n'l'j'} = -5.477 ea_0
<l||er||l'> = = 5.477 ea_0
<j||er||j'> = 6.324 ea_0
<nljmj|er|n'l'j'mj'> = 3.162 ea_0

Cs Rydberg Transition 60S_{1/2}-->60P_{3/2}
R_{nlj-n'l'j'} = 3563.326 ea_0
<l||er||l'> = = -3563.326 ea_0
<j||er||j'> = -4114.575 ea_0
<nljmj|er|n'l'j'mj'> = -2057.287 ea_0

For low-lying states, Numerov integration is usually not accurate enough in estimating dipole matrix element. Package also has a .csv table for each element, that can be manually expanded if needed, with various literature values for dipole matrix elements. These values are either obtained experimentally, or by more advanced theoretical calculations. If value is existing in the literature, package will by default use the literature value with the smallest error estimate.

In [6]:
atom = Caesium()
hasLiteratureValue,dme, info = atom.getLiteratureDME(6,1,1.5,7,0,0.5)
if hasLiteratureValue:
    print("%.2f ea_0" % dme)
    # additional information about the source
    s = "experiment"
    if info[0]==1:
        s = "theory"
    print(" source = %s\n errorEstimate = %.2f\n comment = %s\n reference = %s\n doi = %s" % \
5.61 ea_0
 source = theory
 errorEstimate = 0.01
 comment = scaled, table VI
 reference = Physical Review A, 60 4476 (1999)
 doi = 10.1103/PhysRevA.60.4476

Rydberg Atom Lifetimes

Rydberg atoms have suprisingly long radiative lifetimes. For example lifetimes of $nS$ states are of the order of $10~\mu$s, for $n\sim30$. Lifetimes $\tau$ of $n$ states, with small orbital angular momenutum $\ell$ like $S$ or $P$ states, scales as $\tau\propto n_*^3$ where $n_*=n-\delta$ is the effective principal quantum number. Of all possible $\ell$ states, circular states $\ell=n-1$ have longest lifetimes, that scales as $\tau \propto n_*^5$.

In [7]:
atom = Rubidium()
print("%.2e s" % (atom.getStateLifetime(30,0,0.5)) )
2.43e-05 s

The radiative lifetime is calculated from the Einstein-A coefficients as $1/\tau_0=\sum_{n'\ell'}A_{n\ell\rightarrow n'\ell'}$, with the summation over all dipole coupled states with energy $E_{n'\ell'}<E_{n\ell}$ and the A-coefficient calculated as

\begin{equation} A_{n\ell\rightarrow n'\ell'}=\frac{4\omega_{nn'}^3}{3c^3}\frac{\ell_\mathrm{max}}{2\ell+1}\mathcal{R}^2_{n\ell\rightarrow n'\ell'}. \end{equation}

Here the prefactor of $\omega_{nn'}^3$ means that although the dipole matrix element for transitions to nearby $n$ states is much larger than that for transitions to the ground state (e.g. Cs $\langle 60 S_{1/2}|er| 60 P_{3/2}\rangle = 4114 ~ea_0$ compared to $\langle 6 S_{1/2}|er| 60 P_{3/2}\rangle = 0.0038 ~ea_0$), the A-coefficients for these two transitions are $1.9\mathrm{s}^{-1}$ and $224\mathrm{s}^{-1}$ respectively, resulting in the lowest energy transitions being the dominant decay pathways for radiative decay from the Rydberg states.

At finite temperature, the black-body excitation spectrum causes a reduction in the lifetime due to driving transitions in the GHz and THz region to neighbouring Rydberg states both above and below the energy of the original Rydberg state. The black body loss is given by

\begin{equation} \frac{1}{\tau_\mathrm{BBR}}=\displaystyle\sum_{n'\ell'} \frac{A_{n\ell\rightarrow n'\ell'}}{\exp(\omega_{nn'}/k_\mathrm{B}T)-1}, \end{equation}

with the effective lifetime $1/\tau_\mathrm{eff}=1/\tau_0+1/\tau_\mathrm{BBR}$.

The figure below demonstrates the finite temperature effects by showing the strength of the diffent decay pathways for the $30S_{1/2}$ state in Rb at 300K showing the contribution from radiative decay (red columns) and black-body induced transitions (green). The effective lifetime is reduced by around a factor of 2 at room temperature compared to the natural lifetime.

In [8]:
atom = Rubidium()
pqn = []
y = []
ybb = []
for n in xrange(5,40):
    noBBR = atom.getTransitionRate(30, 0, 0.5, n, 1, 0.5, temperature=0.1)\
            +atom.getTransitionRate(30, 0, 0.5, n, 1, 1.5, temperature=0.1)
    withBBR =  atom.getTransitionRate(30, 0, 0.5, n, 1, 0.5, temperature=300.0)\
            +atom.getTransitionRate(30, 0, 0.5, n, 1, 1.5, temperature=300.0)
y = np.array(y)
ybb = np.array(ybb)

width = 0.4,y,width=width,color="r"),ybb,width=width,color="g")
plt.xlabel("Principal quantum number, $n$")
plt.ylabel(r"Transition rate (s${}^{-1}$)")
plt.title("Transition from 30 $S_{1/2}$ to $n$ $P_{1/2,3/2}$")
plt.legend(("Spontaneous decays","Black-body induced transitions"),fontsize=10)

display(HTML("Lifetime (0 K) &tau;<sub>0</sub> = %.2f &mu;s" % \
             (atom.getStateLifetime(30,0,0.5) *1.e6)))
display(HTML("Lifetime (300 K) &tau;<sub>eff</sub> = %.2f &mu;s" % \
             (atom.getStateLifetime(30,0,0.5,temperature=300.,includeLevelsUpTo=39) *1.e6)))
Lifetime (0 K) τ0 = 24.27 μs
Lifetime (300 K) τeff = 15.54 μs

This plot can be compared with Fig. 1 in Ref. I. I. Beterov et. al, PRA 79, 052504(2009). Similarly, we can reproduce black-body induced state depopulation rates (see Fig. 2 in the same reference):

In [9]:
atom = Rubidium()
stateL = 1; stateJ = 1.5 # P_{1/2} states
BBRdepopulationRate = []
nList = np.arange(6,80)
for n in nList:
    noBBR = 1./atom.getStateLifetime(n,stateL,stateJ,temperature=0)
    withBBR = 1./atom.getStateLifetime(n,stateL,stateJ,temperature=300,includeLevelsUpTo=n+10)
ax = plt.subplot(111)
ax.set_xlabel(r"Principal quantum number, $n$")
ax.set_ylabel(r"BBR induced depopulation rate (s$^{-1}$)")

Comparing lifetimes for different principal quantum numbers $n$ and orbital angular momentum states $l$, we see that lifetimes for $nP$ scale as $\propto n_*^3$, and lifetimes for circular states $n,l=n-1$ scale as $\propto n_*^5$

In [10]:
atom = Rubidium()
pqn = np.arange(10,60,2)

nstar1 = []; nP_lifetime = []
nstar2 = []; circular_state_lifetime = []

for n in pqn:
    nstar1.append(n - atom.getQuantumDefect(n,1,1.0+0.5))

    nstar2.append(n - atom.getQuantumDefect(n,n-1,n-1.5))

nstar1 = np.array(nstar1); nP_lifetime = np.array(nP_lifetime)
nstar2 = np.array(nstar2); circular_state_lifetime = np.array(circular_state_lifetime)
f = plt.figure()
ax = f.add_subplot(1,1,1)

# fitting n^3 dependance for nP state lifetime
def func3(x, a):
    return a+3*x
popt, pcov = curve_fit(func3, np.log(nstar1), np.log(nP_lifetime))

# fitting n^5 dependance for n,l=n-1 streatched state lifetime
def func5(x, a):
    return a+5*x
popt, pcov = curve_fit(func5, np.log(nstar2), np.log(circular_state_lifetime))

ax.set_title("Lifetime scaling")
ax.legend([r"$|n,P,j=3/2\rangle$ states",r"$|n,\ell=n-1,j=n-0.5\rangle$ circular states"],loc=2,fontsize=10)
ax.text(20,1e-3,r"$\propto n_*^5$",fontsize=12)
ax.text(30,3e-5,r"$\propto n_*^3$",fontsize=12)
ax.set_xlabel("Scaled principal quantum number, $n_*$")
ax.set_ylabel("State lifetime (s$^{-1}$)")

Creating highly excited states

Rydberg atoms are highly excited atoms. If the excitation is done in single step necessary energy would typically correspond to UV radition. For example, in Caesium $6 S_{1/2} \rightarrow 60 P_{3/2}$ requred transition wavelength and frequencies are

In [4]:
atom = Caesium()
print("Cs Excitation 6S_{1/2} -> 60P_{3/2}")
n1 = 6; l1 = 0; j1 = 0.5; #Initial State
n2 = 60; l2 = 1; j2 = 1.5; #Final State
print("lambda = %.3f nm" % (atom.getTransitionWavelength(n1,l1,j1,n2,l2,j2)*1e9))
print("omega/2pi = %.3f THz" % (atom.getTransitionFrequency(n1,l1,j1,n2,l2,j2)*1e-12))
Cs Excitation 6S_{1/2} -> 60P_{3/2}
lambda = 318.755 nm
omega/2pi = 940.509 THz

Rabi Frequency

The transition strength is characterised by the Rabi frequency $\Omega=d\cdot\mathcal{E}/\hbar$, where $d$ is the dipole matrix element calcualte above and $\mathcal{E}$ is the electric field of the laser driving the excitation. Recalling that for a Gaussian beam with $1/e^2$ radius of $w_0$ and power $P$, the intensity is given by \begin{equation} I=\frac{2P}{\pi w_0^2}, \end{equation} which is related to the electric field by $I=(1/2)c\epsilon_0\vert \mathcal{E}\vert ^2$, resulting in the Rabi frequency of

\begin{equation} \Omega = \sqrt{\frac{4P}{\epsilon_0c\pi w_0^2}}\frac{\langle n,\ell,j,m_j\vert r_q \vert n',\ell',j',m_j' \rangle}{\hbar}. \end{equation}

Due to the weak overlap between the ground and highly excited Rydberg states (which scales as $\sim n^{*-3/2}$), the resulting matrix elements are relatively small meaning large intensities are required. As an example, consider a 1 W UV laser driving the transition above with a $1/e^2$ waist of 50 $\mu$m. The resulting Rabi frequency is only $\Omega/2\pi=10$ MHz!

In [5]:
#Laser Parameters
waist = 50.e-6 # 50 mu m
P = 1000.e-3 # 500 mW

q=+1; #Light Polarisation (sigma+)

rabiFreq = atom.getRabiFrequency(n1, l1, j1, mj1,
                                 n2, l2, j2,
                                 q, P, waist)
print("rabi Frequency = 2 pi x %.2f MHz" %(rabiFreq/(2*pi)*1e-6))
rabi Frequency = 2 pi x 10.58 MHz

An alternative way to reach these highly excited states is to use a mutli-step excitation scheme with a number of lasers of different wavelength. This enables use of transitions with favourable or convenient wavelengths (for example standard diode/fiber laser sources) with high power or strong transitions. The most commonly used is a two-step excitation via the D-line transition in the alkali's, with the added benefit of enabling electromagnetically induced transparency (EIT) to be engineered for mapping strong Rydberg interactions onto optical photons. In Rubidium, typically the $5 S_{1/2} \rightarrow 5 P_{3/2} \rightarrow 60 S_{1/2}$ transition is favoured, whose properties are

In [13]:
atom = Rubidium()
print("5 S_{1/2} -> 5 P_{3/2}")
print("%.3f nm " % (atom.getTransitionWavelength(5,0,0.5,5,1,1.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(5,0,0.5,5,1,1.5)*1e-12))
print("5 P_{3/2} -> 60 S_{1/2}")
print("%.3f nm " % (atom.getTransitionWavelength(5,1,1.5,60,0,0.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(5,1,1.5,60,0,0.5)*1e-12))
5 S_{1/2} -> 5 P_{3/2}
780.241 nm 
384.230 THz 
5 P_{3/2} -> 60 S_{1/2}
479.839 nm 
624.777 THz 

One can use even more involved schemes, for example in Ceasium three lasers can be used $6 S_{1/2} \rightarrow 6 P_{3/2} \rightarrow 7 S_{1/2} \rightarrow 60 P_{3/2}$. Note that now all the lasers are conveniantly in the near infrared (NIR) range!

In [14]:
atom = Caesium()
print("6 S_{1/2} -> 6 P_{3/2}")
print("%.3f nm " % (atom.getTransitionWavelength(6,0,0.5,6,1,1.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(6,0,0.5,6,1,1.5)*1e-12))
print("6 P_{3/2} -> 7 S_{1/2}")
print("%.3f nm " % (atom.getTransitionWavelength(6,1,1.5,7,0,0.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(6,1,1.5,7,0,0.5)*1e-12))
print("7 S_{1/2} -> 60 P_{3/2}")
print("%.3f nm " % (atom.getTransitionWavelength(7,0,0.5,60,1,1.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(7,0,0.5,60,1,1.5)*1e-12))
6 S_{1/2} -> 6 P_{3/2}
852.347 nm 
351.726 THz 
6 P_{3/2} -> 7 S_{1/2}
1469.892 nm 
203.955 THz 
7 S_{1/2} -> 60 P_{3/2}
779.029 nm 
384.828 THz 

Or, even four lasers can be achieved in order to couple to Rydberg states. For example [reference] $6 S_{1/2} \rightarrow 6P_{3/2} \rightarrow 7S_{1/2} \rightarrow 8 P_{1/2} \rightarrow 52 D_{3/2}$

In [15]:
print("6 S_{1/2} -> 6 P_{3/2}")
print("%.3f nm " % (atom.getTransitionWavelength(6,0,0.5,6,1,1.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(6,0,0.5,6,1,1.5)*1e-12))
print("6 P_{3/2} -> 7 S_{1/2}")
print("%.3f nm " % (atom.getTransitionWavelength(6,1,1.5,7,0,0.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(6,1,1.5,7,0,0.5)*1e-12))
print("7 S_{1/2} -> 8 P_{3/2}")
print("%.3f nm " % (atom.getTransitionWavelength(7,0,0.5,8,1,1.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(7,0,0.5,8,1,1.5)*1e-12))
print("8 P_{1/2} -> 52 D_{3/2}")
print("%.3f nm " % (atom.getTransitionWavelength(8,1,0.5,52,2,1.5)*1e9))
print("%.3f THz " % (atom.getTransitionFrequency(8,1,0.5,52,2,1.5)*1e-12))
6 S_{1/2} -> 6 P_{3/2}
852.347 nm 
351.726 THz 
6 P_{3/2} -> 7 S_{1/2}
1469.892 nm 
203.955 THz 
7 S_{1/2} -> 8 P_{3/2}
1378.174 nm 
217.529 THz 
8 P_{1/2} -> 52 D_{3/2}
1769.006 nm 
169.469 THz 

Now even with moderate laser powers we can reach highly-excited states. For example for 10 mW laser focused down to 50 $\mu$m driving $8~P_{1/2}~m_j=1/2 \rightarrow 52~D_{3/2}~m_j=1/2$ we achieve

In [6]:
waist = 50e-6  # 50 mu m
P = 10.e-3  # 10 mW
q = 0  # laser drives pi transition
rabiFreq = atom.getRabiFrequency(8, 1, 1.5, 0.5, 
                                 52, 2, 1.5,
                                 q, P, waist)
print("rabi Frequency = 2 pi x %.2f MHz" % (rabiFreq / (2*pi) * 1e-6))
rabi Frequency = 2 pi x 3.16 MHz

Quadrupole matrix elements

We can also obtain quadrupole matrix elements. For example, $5 S_{1/2}\rightarrow n D_{5/2}$ transitions, that have been also calculated in D. Tong, PRA 79, 052509 (2009) (see Table I there), we can obtain by calling

In [17]:
atom = Rubidium()
nList = [24,27,29,34,39,44,47,49,54,59]
for n in nList:
    print("%.2e e^2 a_0^4" % atom.getQuadrupoleMatrixElement(5,0,0.5,n,2,2.5)**2)
1.15e-01 e^2 a_0^4
7.94e-02 e^2 a_0^4
6.35e-02 e^2 a_0^4
3.87e-02 e^2 a_0^4
2.52e-02 e^2 a_0^4
1.74e-02 e^2 a_0^4
1.42e-02 e^2 a_0^4
1.25e-02 e^2 a_0^4
9.26e-03 e^2 a_0^4
7.06e-03 e^2 a_0^4

Rydberg Atom Stark Shifts

The large dipole moment of the Rydberg states makes them highly polarisable. For an atom in a d.c. electromagnetic field $\mathcal{E}$ directed along the $z$-axis, the interaction leads to a coupling given by

\begin{equation} \mathcal{H} = \mathcal{H}_0 + \mathcal{E}\hat{z}, \end{equation}

where $\mathcal{H}_0$ is the Hamiltonian for the unperturbed energy levels and $\hat{z}\equiv r_0$ in the spherical basis resulting in a selection rule $\Delta m_j=0$ and a shift dependent on $\vert m_j \vert$. For the low-$\ell$ states, the quantum defects cause the first-order Stark-shift to vanish, resulting in a second-order energy shift $\Delta E = -\frac{1}{2}\alpha_0\mathcal{E}^2$ where the polarizability for state $\vert n \ell j m_j\rangle$ is given by

\begin{equation} \alpha_0 = 2 e^2\displaystyle\sum_{n'\ell'j'\neq n\ell j} \frac{\vert \langle n,\ell,j,m_j\vert r_0 \vert n',\ell',j',m_j' \rangle \vert^2}{E_{n'\ell' j'}-E_{n\ell j}}. \end{equation}

Since the dipole matrix elements scale as $n^{2}$, and the energy difference scales as $n^{-3}$ the scalar polarizability $\alpha_0\propto n^7$. States with $\ell>3$ are degenerate, and therefore experience a linear Stark shift.

The evaluation of polarisability using second-order perturbation theory is limited, and only valid at weak fields. Instead, Stark maps accurate to all orders can be obtained by direct diagonalisatoin of the Hamiltonian $\mathcal{H}$, taking care to include a large enough set of basis states to achieve divergence. As well as revealing the exact energies of the Stark-coupled states, this permits extraction of the mixed electric field eigenstates which are relevant when considering optical excitation of a given Rydberg state which is now mixed into neighbouring energy levels.

The code below provides an example Stark map for the $35S_{1/2},\vert m_j\vert=1/2$ state in Cesium. After identifying the target state, a range of basis states are defined before running the diagonalisation for the range of E-fields specified, returning the desired Stark map. To show the effect of state mixing, the levels are coloured proportional to the fraction of the target state present in each eigenstate. This knowledge can also be used to determine the effective polarisability by extracting the energy level of the eigenstate with a contribution of at least 90% of the original target state.

For high $n$, the range of basis states may need to be increased due to the reduction in energy spacings but this can easily be checked for convergence by running twice with a larger basis set.

In [18]:
#Stark Map Caclulator
#Initialise a Stark-shift Solver for Caesium
calc = StarkMap(Caesium())

#Target state
#Define max/min n values in basis
#Maximum value of l to include (l~20 gives good convergence for states with l<5)

#Initialise Basis States for Solver : progressOutput=True gives verbose output
calc.defineBasis(n0, l0, j0, mj0, nmin, nmax, lmax, progressOutput=True)

Emin=0. #Min E field (V/m)
Emax=20000. #Max E field (V/m)
N=1001 #Number of Points

#Generate Stark Map
calc.diagonalise(np.linspace(Emin,Emax,N), progressOutput=True)
#Show Sark Map
calc.plotLevelDiagram(progressOutput=True,units=1,highlighState = True)
calc.showPlot(interactive = False)
#Return Polarizability of target state    
print("%.5f MHz cm^2 / V^2 " % calc.getPolarizability(showPlot=True, minStateContribution=0.9))
Found  410  states.
Generating matrix...

Finding eigenvectors...

4.15576 MHz cm^2 / V^2 

We can compare values for the polarizability with the measurements in the literature. For example, we can get theoretical predictions for results cited in Ref. M. S. O'Sullivan and B. P. Stoicheff, PRA 31, 2718 (1985) (see Table I there) by calling the following

In [20]:
calc = StarkMap(Rubidium())

stateN = [15,20,25,30,35,40,45,50,55,60,63,65,67,70,75,80]
eFieldRange = np.array([[90,492],[28,188],[19,75],[4,28],[4,14],\
                        [0.003,0.2],[0.1,0.16]])*1.e2 # V/m ranges
polarizabilityList = []
print("State\t\tPolarizability (MHz cm^2 / V^2)\t\tElectric field range (V/cm)")

for i in xrange(len(stateN)):
    n = stateN[i]
    calc.defineBasis(n, 0, 0.5, 0.5, n-5, n+5, 20) #,progressOutput=True)
    minEfield = eFieldRange[i][0]
    maxEfield = eFieldRange[i][1]
    calc.diagonalise(np.linspace(minEfield,maxEfield,100)) # ,progressOutput=True)
    p = calc.getPolarizability( minStateContribution=0.9)
    print("%s\t%.3e\t\t\t\t%.2f-%.2f" % \
        (printStateString(n,0, 0.5), p, minEfield/100.,maxEfield/100.))
State		Polarizability (MHz cm^2 / V^2)		Electric field range (V/cm)
15 S 1/2	7.995e-03				90.00-492.00
20 S 1/2	7.233e-02				28.00-188.00
25 S 1/2	3.755e-01				19.00-75.00
30 S 1/2	1.396e+00				4.00-28.00
35 S 1/2	4.208e+00				4.00-14.00
40 S 1/2	1.067e+01				1.00-5.00
45 S 1/2	2.457e+01				1.40-3.20
50 S 1/2	5.148e+01				0.70-2.00
55 S 1/2	9.957e+01				0.60-1.00
60 S 1/2	1.821e+02				0.15-0.70
63 S 1/2	2.611e+02				0.35-0.70
65 S 1/2	3.251e+02				0.20-0.70
67 S 1/2	3.955e+02				0.15-0.45
70 S 1/2	5.323e+02				0.15-0.30
75 S 1/2	8.532e+02				0.00-0.20
80 S 1/2	1.353e+03				0.10-0.16

This data can be compared also with analytical fit function that they recommend as a polarizability estimate.

In [21]:
ax = plt.subplot(1,1,1)
ax.set_xlabel(r"Principal quantum number, $n$")
ax.set_ylabel(r"Scalar polarizability, $\alpha_0$ (MHz cm$^2$/ V$^2$)")

def functionFromReference_Rb_S_states(atom,n):
    # recommended fit function from M. S. O'Sullivan and B. P. Stoicheff, PRA 31, 2718 (1985) 
    scaledPQN = n - atom.getQuantumDefect(n,0,0.5)
    return 2.202*1.e-9*scaledPQN**6 + 5.53*1.e-11*scaledPQN**7

nValues = np.arange(10,80)
referenceFunction = [functionFromReference_Rb_S_states(calc.atom,n) for n in nValues]
plt.legend(("calculation","reference function"),loc=4,fontsize=10)

We can also compare our calculations with recent calculations and measurements from J. Grimmel, NJP 17, 053005 (2015) for Rubidium $35 S_{1/2}$ and $70 S_{1/2}$ states

In [23]:
calc = StarkMap(Rubidium())

nRange = 13
noOfEigVectors = 500

calc.defineBasis(35, 2, 2.5, 0.5, 35-nRange, 35+nRange, 35,progressOutput=True)
calc.diagonalise(np.linspace(0,235*1.e2,noOfEigVectors),drivingFromState=[5, 1, 1.5, 1.5, -1],\
calc.plotLevelDiagram(units=2,highlighState = True,highlightColour="red")

calc.defineBasis(35, 2, 2.5, 1.5, 35-nRange, 35+nRange, 35,progressOutput=True)
calc.diagonalise(np.linspace(0,235*1.e2,noOfEigVectors),drivingFromState=[5, 1, 1.5, 1.5, 0],\
calc.plotLevelDiagram(units=2,highlighState = True,highlightColour="green", addToExistingPlot=True)

calc.defineBasis(35, 2, 2.5, 2.5, 35-nRange, 35+nRange, 35,progressOutput=True)

calc.diagonalise(np.linspace(0,235*1.e2,noOfEigVectors),drivingFromState=[5, 1, 1.5, 1.5, +1],\
calc.plotLevelDiagram(units=2,highlighState = True,highlightColour="blue", addToExistingPlot=True)

stateEnergy = calc.atom.getEnergy(35,0,0.5)*C_e/(C_h*1.e9),stateEnergy+7)
calc.showPlot(interactive = False)
Found  1636  states.
Generating matrix...

Finding driving field coupling...

Finding eigenvectors...

Found  1584  states.
Generating matrix...

Finding driving field coupling...

Finding eigenvectors...

Found  1532  states.
Generating matrix...

Finding driving field coupling...

Finding eigenvectors...


Rydberg atom interactions

As described above, the Rydberg states have extremely large dipole moments making them ideal candidates for controlled long-range interactions in quantum information processing and for studying effects such as resonant energy transfer or superradiance.

For two atoms separated by distance $\textbf{R}$ the dipole-dipole interaction is given by

\begin{equation} V(\textbf{R}) = \frac{\textbf{$\mu$}_1\cdot\textbf{$\mu$}_2}{R^3} -\frac{3(\textbf{$\mu$}_1\cdot\textbf{R})(\textbf{$\mu$}_2\cdot\textbf{R})}{R^5}, \end{equation}

where the dipole matrix elements $\mathbf{\mu}_{1,2}$ describe transitions from the initial Rydberg states $\vert r \rangle$ to other dipole-coupled states $\vert r'\rangle,\vert r''\rangle$ respectively. The dipole-coupled pair-states have an energy difference $\Delta$ give by

\begin{equation} \Delta = E_{r'}+E_{r''}-2E_{r}, \end{equation}

with the dominant contribution to the atom-atom interaction arising from the pair state of with smallest absolute defect $\vert \Delta \vert$. Reducing this to a simple two-pair problem, the Hamiltonian for the states $\vert r r\rangle,\vert r'r''\rangle$ is given by

\begin{equation} \mathcal{H} = \begin{pmatrix} 0 &V(R) \\ V(R)&\Delta \end{pmatrix}, \end{equation}

which has the eigenvalues

\begin{equation} \lambda_\pm = \frac{\Delta\pm\sqrt{\Delta^2+4V(R)^2}}{2}. \end{equation}

This results in two asympotitic limits for the atom-atom interactions:

i) Long-range ($V(R)\ll\Delta$)

In this limit, also known as the Van der Waals regime, pair-state energy is shifted by $-V(R)^2/\Delta\equiv - C_6/R^6$. The sign of the interaction (attractive/repulsive) is determined by $\Delta$, and the interaction scales as $C_6\propto n^{11}$.

ii) Short-range ($V(R)\gg\Delta$)

This is known as the resonant dipole-dipole regime with $\lambda=\pm C_3/R^3$ where $C_3\propto n^4$, resulting in significant state mixing.

The transition between these two regimes is known as the Van der Waals radius occuring at $V(R_\mathrm{vdW})=\Delta$, with $R_\mathrm{vdW}=\sqrt[6]{C_6/\vert\Delta\vert}$.

Typically $\Delta\neq0$ for the Rydberg states, however for small $\Delta$ an external electric field can be used to Stark-shift the pair states into resonance known as a Förster resonance, resulting in resonant $C_3/R^3$ behaviour at all radii.

In [24]:
# calcualtion of energy defects and c3 terms for atom-pairs
atom = Rubidium()


fig, axes = plt.subplots(1, 1, figsize=(10,6))

l1 = 1 ; j1 = 1.5 #  P_{3/2}
l2 = 0 ; j2 = 0.5 # S_{1/2}
energyDefect = []
for n in nlist:
              |(n+1) %s_{%d/2},n %s_{%d/2}\rangle$"%\

l1 = 1 ; j1 = 0.5 #  P_{3/2}
l2 = 0 ; j2 = 0.5 # S_{1/2}
energyDefect = []
for n in nlist:
              |(n+1) %s_{%d/2},n %s_{%d/2}\rangle$"%\

l1 = 1 ; j1 = 1.5 #  P_{3/2}
l2 = 0 ; j2 = 0.5 # S_{1/2}
energyDefect = []
for n in nlist:
              |(n+2) %s_{%d/2},(n_1) %s_{%d/2}\rangle$"%\

l1 = 1 ; j1 = 0.5 #  P_{3/2}
l2 = 0 ; j2 = 0.5 # S_{1/2}
energyDefect = []
for n in nlist:
              |(n+2) %s_{%d/2},(n_1) %s_{%d/2}\rangle$"%\

axes.set_ylabel("Energy defect, $\Delta$ (GHz)")
axes.set_xlabel("Principal quantum number, $n$")

The obtained plot can be compared with Fig. 8 (for Rb) in Ref. Thad G. Walker and M. Saffman, PRA 77, 032723 (2008). Similar paper provides $C_6$ calculations for individual coulping channels. We can reproduce examples corresponding to Eq. (50a-c) in the same reference simply by calling:

In [25]:
channels = [[70,0,0.5, 70, 1,1.5, 69,1, 1.5],\
            [70,0,0.5, 70, 1,1.5, 69,1, 0.5],\
            [70,0,0.5, 69, 1,1.5, 70,1, 0.5],\
            [70,0,0.5, 70, 1,0.5, 69,1, 0.5]]
print(" = = = Caesium = = = ")
atom = Caesium()
for channel in channels:
    print("%.0f  GHz (mu m)^6" % ( atom.getC6term(*channel)/C_h*1.e27 ))
print("\n = = = Rubidium  = = =")
atom = Rubidium()
for channel in channels:
    print("%.0f  GHz (mu m)^6" % ( atom.getC6term(*channel)/C_h*1.e27 ))
 = = = Caesium = = = 
722  GHz (mu m)^6
316  GHz (mu m)^6
382  GHz (mu m)^6
227  GHz (mu m)^6

 = = = Rubidium  = = =
799  GHz (mu m)^6
543  GHz (mu m)^6
589  GHz (mu m)^6
437  GHz (mu m)^6

Dispersion Coefficients

For the long-range limit, the $C_6$ dispersion coefficient can be calculated using second-order petubation theory valid for $R>R_\mathrm{vdW}$ and taking a sum over all dipole-coupled pair-states using the equation

\begin{equation} C_6 = \displaystyle\sum_{r'r''} \frac{\left| \langle r'r''\vert V(R) \vert rr\rangle\right|^2}{\Delta_{r',r''}}. \end{equation}

As well as the dependence upon the sign of $\Delta$ for attractive/repulsive interactions, the form of the dipole interaction above is anisotropic leading to a strong angular sensitivity. Defining an angle $(\theta,\phi)$ in spherical basis between the quantization axis and the interatomic separation vector $\textbf{R}$, the plots below highlight the difference in interaction strength of the different $\ell$ states. As a general rule, for the alkali atoms $ns$ states are repulsive ($C_6<0$) and relatively isotropic, whilst $nd$ states are attractive ($C_6>0$) and highly anisotropic.

Note: Typically, this sum is constrained to include only pair states up to a maximum value of $\Delta$ to minimise computational overhead, and this can be set to achieve convergence.

In [2]:
#Dipole-Interaction Dispersion Coefficient:60S1/2
#Evaluation of the Cs 60S_1/2 C6 coefficient using perturbation theory (Theta=0,phi=0)
n0=60;l0=0;j0=0.5;mj0=0.5; #Target State
theta=0; #Polar Angle [0-pi]
phi=0; #Azimuthal Angle [0-2pi]
dn = 5; #Range of n to consider (n0-dn:n0+dn)
deltaMax = 25e9 #Max pair-state energy difference [Hz]

#Set target-state and extract value
calculation = PairStateInteractions(Rubidium(), n0,l0,j0,n0,l0,j0, mj0,mj0)
C6 = calculation.getC6perturbatively(theta,phi, dn, deltaMax)
print("C6 [%s] = %.2f GHz (mum)^6" % (printStateString(n0,l0,j0),C6))

#Angular Coupling
#Evaluate C6 as a function of angle
thetaList = np.linspace(0,2*pi,55)
ax = plt.subplot(111, projection='polar')
line = []

c6 = []
for t in thetaList:
    C6 = calculation.getC6perturbatively(t,phi, dn, deltaMax)
# plot results
lineLegend, = plt.plot(thetaList,c6,"-",color="r",label=("mj=%d/2"%int(2*mj0)) )
plt.title("Rb, pairstate  60 $S_{1/2}$, $|C_6|$ in GHz $\mu$m$^{6}$",fontsize=10)
C6 [60 S 1/2] = -138.88 GHz (mum)^6
Note: No saved angular matrix files to be loaded.
(<type 'exceptions.KeyboardInterrupt'>, KeyboardInterrupt(), <traceback object at 0x10434f9e0>)
In [25]:
#Dipole-Interaction Dispersion Coefficient:60SD5/2
n0=60;l0=2;j0=2.5;#Target State
phi=0; #Azimuthal Angle [0-2pi]
dn = 5; #Range of n to consider (n0-dn:n0+dn)
deltaMax = 25e9 #Max pair-state energy difference [Hz]

thetaList = np.linspace(0,2*pi,71)
mj = [0.5,1.5,2.5]
colourList = ["b","g","r"]

ax = plt.subplot(111, projection='polar')
line = []
for i in [0,1,2]:
    calculation = PairStateInteractions(Rubidium(), n0,l0,j0,n0,l0,j0, mj[i], mj[i])
    c6 = []
    for t in thetaList:
        C6=calculation.getC6perturbatively(t,phi, dn, deltaMax);
    # plot results
    lineLegend, = plt.plot(thetaList,c6,"-",color=colourList[i],label=("mj=%d/2"%int(2*mj[i])))
plt.title("Rb, pairstate  60 $D_{5/2}$, $|C_6|$ in GHz $\mu$m$^{6}$",fontsize=10)