ARC Additions: Notebook of new functions and testing / benchmarking to verify parameters and sources#

Goal: Transfer hyperfine atomic physics functions from internal Matlab toolbox to ARC.

Jonathan Pritchard, University of Strathclyde (2020)

[8]:
#Configure Environment

from arc import *

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

cs = Cesium()
rb87 = Rubidium87()

Check ARC Parameters for D1/2 and inverted (n+1)P_{1/2,3/2}#

  • New function A, B = self.getHFSCoefficients(self, n, l, j, s=None) added to return HFS constants in Hz. Refernce [5] and [6] have useful reference info

Hyperfine Structure Functions#

  • A, B = atom.getHFSCoefficients(n,l,j)

Returns HFS constants in Hz. This is currently using hard-coded values for Cs and Rb, but has been written as a place-holder for adding in a table. Refernce [5] and [6] have useful reference info.

  • getHFSEnergyShift(j, f, A, B)

Returns HFS Energy shift from centre of mass.

\begin{equation} \Delta E_\mathrm{hfs} = \frac{A}{2}K+B\frac{\frac{3}{2}K(K+1)-2I(I+1)J(J+1)}{2I(2I-1)2J(2J-1)}, \end{equation} where \(K=F(F+1)-I(I+1)-J(J+1)\).

References#

[1] Steck

[2] M. Ortiz and J. Campos, J. Quant. Spectrosc. Radiat. Transfer 26, 107 (1981)

[3] Vasilyev, A. A. and Savukov, I. M. and Safronova, M. S. and Berry, H. G., PRA 66, 020101 (2002)

[4] M. S. Safronova, Carl J. Williams, and Charles W. Clark, PRA 69, 022509 (2004)

[5] D. Feiertag, A. Sahm, and G. zu Putlitz, Z. Phys. 255, 93 (1972)

[6] E. Arrimondo, M. Inguscio and P. Violino, Rev. Mod. Phys. 49, 31 (1977)

[2]:
print("Cs Numbers")
print("==========")
print("State Lifetimes:")
# Lifetime Comparison -->
print("\t[1] : Tau6P1/2 = 34.894(44) ns, ARC = %2.3f ns (%2.2f %% error)"
      % (cs.getStateLifetime(6, 1, 0.5) * 1e9,
         (cs.getStateLifetime(6, 1, 0.5) * 1e9 - 34.894) / 34.894 * 100)
         )
print("\t[1] : Tau6P3/2 = 30.473(44) ns, ARC = %2.3f ns (%2.2f %% error)"
      % (cs.getStateLifetime(6, 1, 1.5) * 1e9,
         (cs.getStateLifetime(6, 1, 1.5) * 1e9 - 30.473) / 30.473 * 100)
        )
print("\t[2] : Tau7P1/2 = 155 ns, ARC = %2.3f ns (%2.2f %% error)"
      % (cs.getStateLifetime(7, 1, 0.5) * 1e9,
         (cs.getStateLifetime(7, 1, 0.5) * 1e9 - 155) / 155 * 100)
        )
print("\t[2] : Tau7P3/2 = 133 ns, ARC = %2.3f ns (%2.2f %% error)"
      % (cs.getStateLifetime(7, 1, 1.5) * 1e9,
         (cs.getStateLifetime(7, 1, 1.5) * 1e9 - 133) / 133 * 100)
         )
# Matrix Elements
print("Matrix Elements")
r = 4.5003; rm = cs.getReducedMatrixElementJ(6, 0, 0.5, 6, 1, 0.5)
print("\t[1] : <6S1/2||er||6P1/2> = %1.4f, ARC = %2.4f ns (%2.2f %% error)"
      % (r, rm, (rm - r) / r * 100)
     )
r = 6.3338; rm = cs.getReducedMatrixElementJ(6, 0, 0.5, 6, 1, 1.5)
print("\t[1] : <6S1/2||er||6P3/2> = %1.4f, ARC = %2.4f ns (%2.2f %% error)"
      % (r, rm, (rm - r) / r * 100)
      )
r = 0.2757; rm = cs.getReducedMatrixElementJ(6, 0, 0.5, 7,1,0.5)
print("\t[3] : <6S1/2||er||7P1/2> = %1.4f, ARC = %2.4f ns (%2.2f %% error)"
      % (r, rm, (rm - r) / r * 100)
      )
r = 0.5856; rm = cs.getReducedMatrixElementJ(6, 0, 0.5, 7, 1, 1.5)
print("\t[3] : <6S1/2||er||7P3/2> = %1.4f, ARC = %2.4f ns (%2.2f %% error)"
      % (r, rm, (rm - r) / r * 100)
      )
# Hyperfine Splitting Coefficients
print("Hyperfine Constants")
A, B = cs.getHFSCoefficients(6, 0, 0.5)
omegaQ = 4 * A
print("\t[1] : 6S1/2 A = %1.4f MHz,B = %1.4f MHz" % (A * 1e-6, B * 1e-6))
A, B = cs.getHFSCoefficients(6, 1, 0.5)
print("\t[1] : 6P1/2 A = %1.4f MHz,B = %1.4f MHz" % (A * 1e-6, B * 1e-6))
A, B = cs.getHFSCoefficients(6, 1, 1.5)
print("\t[1] : 6P3/2 A = %1.4f MHz,B = %1.4f MHz" % (A * 1e-6, B * 1e-6))
A, B = cs.getHFSCoefficients(7, 1, 0.5)
print("\t[5] : 7P1/2 A = %1.4f MHz,B = %1.4f MHz" % (A * 1e-6, B * 1e-6))
A, B = cs.getHFSCoefficients(7, 1, 1.5)
print("\t[6] : 7P3/2 A = %1.4f MHz,B = %1.4f MHz" % (A * 1e-6, B * 1e-6))

print("\nVerify getHFSCoefficients(j,f,A,B) using Cs 6P3/2")
A, B = cs.getHFSCoefficients(6, 1, 1.5)
Ehfs = cs.getHFSEnergyShift(1.5, 5, A, B)
print("\t[1] : 6P3/2 F=5: Ehfs = +263.81 MHz --> ARC Ehfs = %f MHz" % (Ehfs * 1e-6))
Ehfs = cs.getHFSEnergyShift(1.5, 4, A, B)
print("\t[1] : 6P3/2 F=4: Ehfs =  +12.815 MHz --> ARC Ehfs = %f MHz" % (Ehfs * 1e-6))
Ehfs = cs.getHFSEnergyShift(1.5, 3, A, B)
print("\t[1] : 6P3/2 F=3: Ehfs = -188.44 MHz --> ARC Ehfs = %f MHz" % (Ehfs * 1e-6))
Ehfs = cs.getHFSEnergyShift(1.5, 2, A, B)
print("\t[1] : 6P3/2 F=2: Ehfs = -339.64 MHz --> ARC Ehfs = %f MHz" % (Ehfs * 1e-6))


# Clock Frequency
print("\n\t[1] : 6S1/2 Clock Frequency = 9.192631770 GHz --> ARC = %f GHz (%2.2f %% error)"
      % (omegaQ * 1e-9, (omegaQ - 9.192631770e9) / 9.192631770e9 * 100)
     )

print("\nRb Numbers")
print("==========")
print("State Lifetimes:")

# Lifetime Comparison -->
print("\t[1] : Tau5P1/2 = 27.70(4) ns, ARC = %2.3f ns (%2.2f %% error)"
      % (rb87.getStateLifetime(5, 1, 0.5) * 1e9,
         (rb87.getStateLifetime(5, 1, 0.5) * 1e9 - 27.70) / 27.70 * 100)
      )
print("\t[1] : Tau5P3/2 = 26.24(4) ns, ARC = %2.3f ns (%2.2f %% error)"
      % (rb87.getStateLifetime(5, 1, 1.5) * 1e9,(rb87.getStateLifetime(5, 1, 1.5) * 1e9 - 26.24) / 26.24 * 100))
print("\t[4] : Tau6P1/2 = 125 ns, ARC = %2.3f ns (%2.2f %% error)"
       % (rb87.getStateLifetime(6, 1, 0.5) * 1e9,
         (rb87.getStateLifetime(6, 1, 0.5) * 1e9 - 125) / 125 * 100)
      )
print("\t[4] : Tau6P3/2 = 112 ns, ARC = %2.3f ns (%2.2f %% error)" %
      (rb87.getStateLifetime(6, 1, 1.5) * 1e9,
       (rb87.getStateLifetime(6, 1, 1.5) * 1e9 - 112) / 112 * 100)
      )

# Matrix Elements
print("Matrix Elements")
r = 4.2313
rm = rb87.getReducedMatrixElementJ(5, 0, 0.5, 5, 1, 0.5)
print("\t[4] : <5S1/2||er||5P1/2> = %1.4f, ARC = %2.4f ns (%2.2f %% error)"
      % (r, rm, (rm - r) / r * 100)
     )
r = 5.9780
rm = rb87.getReducedMatrixElementJ(5, 0, 0.5, 5, 1, 1.5)
print("\t[4] : <5S1/2||er||5P3/2> = %1.4f, ARC = %2.4f ns (%2.2f %% error)"
      % (r, rm, (rm - r) / r * 100)
     )
r = 0.333
rm = rb87.getReducedMatrixElementJ(5, 0, 0.5, 6, 1, 0.5)
print("\t[4] : <5S1/2||er||6P1/2> = %1.4f, ARC = %2.4f ns (%2.2f %% error)"
      % (r, rm, (rm - r) / r * 100)
     )
r = 0.541
rm = rb87.getReducedMatrixElementJ(5, 0, 0.5, 6, 1, 1.5)
print("\t[4] : <5S1/2||er||6P3/2> = %1.4f, ARC = %2.4f ns (%2.2f %% error)"
      % (r, rm, (rm - r) / r * 100)
     )
print("Hyperfine Constants")
A, B = rb87.getHFSCoefficients(5, 0, 0.5)
omegaQ = (0.5 + rb87.I) * A
print("\t[1] : 5S1/2 A = %1.4f MHz,B = %1.4f MHz"
      % (A * 1e-6, B * 1e-6))
A, B = rb87.getHFSCoefficients(5, 1, 0.5)
print("\t[1] : 5P1/2 A = %1.4f MHz,B = %1.4f MHz"
      % (A * 1e-6, B * 1e-6))
A, B = rb87.getHFSCoefficients(5, 1, 1.5)
print("\t[1] : 5P3/2 A = %1.4f MHz,B = %1.4f MHz"
      % (A * 1e-6, B * 1e-6)
     )
A, B = rb87.getHFSCoefficients(6, 1, 0.5)
print("\t[5] : 6P1/2 A = %1.4f MHz,B = %1.4f MHz"
      % (A * 1e-6, B * 1e-6)
     )
A, B = rb87.getHFSCoefficients(6, 1, 1.5)
print("\t[6] : 6P3/2 A = %1.4f MHz,B = %1.4f MHz"
      % (A * 1e-6, B * 1e-6)
     )

# Clock Frequency
print("\n\t[1] : 5S1/2 Clock Frequency =  6.83468261090429 GHz --> ARC = %f GHz (%2.2f %% error)"
      % (omegaQ * 1e-9, (omegaQ - 6.83468261090429e9) / 6.83468261090429e9 * 100)
      )
Cs Numbers
==========
State Lifetimes:
        [1] : Tau6P1/2 = 34.894(44) ns, ARC = 34.894 ns (0.00 % error)
        [1] : Tau6P3/2 = 30.473(44) ns, ARC = 30.474 ns (0.00 % error)
        [2] : Tau7P1/2 = 155 ns, ARC = 154.368 ns (-0.41 % error)
        [2] : Tau7P3/2 = 133 ns, ARC = 129.775 ns (-2.42 % error)
Matrix Elements
        [1] : <6S1/2||er||6P1/2> = 4.5003, ARC = 4.5003 ns (0.00 % error)
        [1] : <6S1/2||er||6P3/2> = 6.3338, ARC = 6.3337 ns (-0.00 % error)
        [3] : <6S1/2||er||7P1/2> = 0.2757, ARC = 0.2757 ns (0.00 % error)
        [3] : <6S1/2||er||7P3/2> = 0.5856, ARC = 0.5856 ns (0.00 % error)
Hyperfine Constants
        [1] : 6S1/2 A = 2298.1579 MHz,B = 0.0000 MHz
        [1] : 6P1/2 A = 291.9201 MHz,B = 0.0000 MHz
        [1] : 6P3/2 A = 50.2750 MHz,B = -0.5300 MHz
        [5] : 7P1/2 A = 94.3500 MHz,B = 0.0000 MHz
        [6] : 7P3/2 A = 16.6050 MHz,B = -0.1500 MHz

Verify getHFSCoefficients(j,f,A,B) using Cs 6P3/2
        [1] : 6P3/2 F=5: Ehfs = +263.81 MHz --> ARC Ehfs = 263.811250 MHz
        [1] : 6P3/2 F=4: Ehfs =  +12.815 MHz --> ARC Ehfs = 12.814821 MHz
        [1] : 6P3/2 F=3: Ehfs = -188.44 MHz --> ARC Ehfs = -188.436607 MHz
        [1] : 6P3/2 F=2: Ehfs = -339.64 MHz --> ARC Ehfs = -339.640179 MHz

        [1] : 6S1/2 Clock Frequency = 9.192631770 GHz --> ARC = 9.192632 GHz (0.00 % error)

Rb Numbers
==========
State Lifetimes:
        [1] : Tau5P1/2 = 27.70(4) ns, ARC = 27.690 ns (-0.04 % error)
        [1] : Tau5P3/2 = 26.24(4) ns, ARC = 26.238 ns (-0.01 % error)
        [4] : Tau6P1/2 = 125 ns, ARC = 129.295 ns (3.44 % error)
        [4] : Tau6P3/2 = 112 ns, ARC = 118.143 ns (5.48 % error)
Matrix Elements
        [4] : <5S1/2||er||5P1/2> = 4.2313, ARC = 4.2321 ns (0.02 % error)
        [4] : <5S1/2||er||5P3/2> = 5.9780, ARC = 5.9783 ns (0.01 % error)
        [4] : <5S1/2||er||6P1/2> = 0.3330, ARC = 0.3330 ns (0.00 % error)
        [4] : <5S1/2||er||6P3/2> = 0.5410, ARC = 0.5410 ns (0.00 % error)
Hyperfine Constants
        [1] : 5S1/2 A = 3417.3413 MHz,B = 0.0000 MHz
        [1] : 5P1/2 A = 408.3280 MHz,B = 0.0000 MHz
        [1] : 5P3/2 A = 84.7185 MHz,B = 12.4965 MHz
        [5] : 6P1/2 A = 132.5600 MHz,B = 0.0000 MHz
        [6] : 6P3/2 A = 27.6740 MHz,B = 3.9450 MHz

        [1] : 5S1/2 Clock Frequency =  6.83468261090429 GHz --> ARC = 6.834683 GHz (-0.00 % error)

Angular Momentum Functions#

  • getSphericalDipoleMatrixElement(j1, mj1, j2, mj2, q)

    Angular part of spherical matrix element \(A\) in units of reduced matrix element \(\langle j\vert\vert\mu\vert\vert j'\rangle\) (Wigner-Eckart decomposition)

\begin{equation} \langle j m_j\vert \mu_q \vert j' m_j' \rangle = (-1)^{j-m_j}\begin{pmatrix}j&1&j'\\-m_j&-q&m_j'\end{pmatrix}\langle j\vert\vert\mu\vert\vert j'\rangle = A\langle j\vert\vert\mu\vert\vert j'\rangle, \end{equation}

  • _reducedMatrixElementFJ(j1, f1, j2, f2)

Reduced dipole matrix element \(\langle f \vert\vert\mu\vert\vert f'\rangle\) in units of \(\langle j\vert\vert\mu\vert\vert j'\rangle\)

\begin{equation} \langle f \vert\vert\mu\vert\vert f'\rangle = (-1)^{j+I+f'+1}\sqrt{(2f+1)(2f'+1)} \begin{Bmatrix}f&1&f'\\j'&I&j\end{Bmatrix}\langle j\vert\vert\mu\vert\vert j'\rangle, \label{eq:frf} \end{equation}

  • getSphericalMatrixElementHFStoFS(j1, f1, mf1, j2, mj2, q)

Angular coefficient of matrix element \(\langle jfm_f \vert er_q\vert j'f'm_f'\rangle\) for coupling to state without resolved HFS, in units of \(\langle j\vert\vert\mu\vert\vert j'\rangle\)

\begin{align} \langle f,m_f \vert\mu_q\vert j',m_j'\rangle &=\displaystyle\sum_{f'}C_{j'm_j'I(m_f+q-m_j')}^{f'(m_f+q)}(-1)^{f-m_f}\begin{pmatrix}f&1&f'\\-m_f&-q&m_f+q\end{pmatrix}\\& \times(-1)^{j+I+f'+1}\sqrt{(2f+1)(2f'+1)} \begin{Bmatrix}f&1&f'\\j'&I&j\end{Bmatrix}\langle j\vert\vert\mu\vert\vert j'\rangle \end{align}

  • getMagneticDipoleMatrixElementHFS(l, j, f1, mf1, m2, mf2, q)

Returns spherical component of magnetic dipole transition from \(\vert f_1,m_{f_1}\rangle\rightarrow\vert f_2,m_{f_2}\rangle\) in units of \(\mu_B B_q\) corresponding to relative strength of microwave transitions between hyperfine states of the same \(n\ell j\) state driven by magnetic field with polarisation \(q\).

\begin{align} \langle f_1,m_{f_1}\vert \mu_q \vert f_2,m_{f_2}\rangle = g_J \mu_B B_q (-1)^{f_2+j+I+1+f_1-m_{f_1}}\sqrt{(2f_1+1)(2f_2+1)j(j+1)(2j+1)} \begin{pmatrix}f_1&1&f_2\\-m_{f_1} & -q & m_{f_2}\end{pmatrix}\begin{Bmatrix}f_1&1&f_2\\j & I & j\end{Bmatrix} \end{align}

  • getDipoleMatrixElementHFStoFS(n1, l1, j1, f1, mf1, n2, l2, j2, mj2, q)

Full dipole matrix element \(\langle jfm_f \vert er_q\vert j'f'm_f'\rangle\) for coupling to state without resolved HFS, in units of \(ea_0\)

\begin{align} \langle f,m_f \vert\mu_q\vert j',m_j'\rangle &=\displaystyle\sum_{f'}C_{j'm_j'I(m_f+q-m_j')}^{f'(m_f+q)}(-1)^{f-m_f}\begin{pmatrix}f&1&f'\\-m_f&-q&m_f+q\end{pmatrix}\\& \times(-1)^{j+I+f'+1}\sqrt{(2f+1)(2f'+1)} \begin{Bmatrix}f&1&f'\\j'&I&j\end{Bmatrix}\langle j\vert\vert\mu\vert\vert j'\rangle \end{align}

  • getBranchingRatio(jg, fg, mfg, je, fe, mfe)

Branching ratio \(b\) for fraction of spontaneous decay from \(\vert{j_e,f_e,m_{f_e}}\rangle\rightarrow\vert{j_g,f_g,m_{f_g}}\rangle\).

\begin{align} b = \displaystyle\sum_q (2j_e+1)\left(\begin{matrix}f_g & 1 & f_e \\ -m_{f_g} & q & m_{f_e}\end{matrix}\right)^2 \vert \langle j_e,f_e\vert \vert er \vert\vert j_g,f_g\rangle\vert^2 \end{align}

  • getSaturationIntensity(ng, lg, jg, fg, mfg, ne, fe, je, fe, mfe)

Returns saturation intenstiy \(I_\mathrm{sat}\) for transition from \(\vert{j_g,f_g,m_{f_g}}\rangle\rightarrow\vert{j_e,f_e,m_{f_e}}\rangle\) in units of W/m\(^2\).

\begin{align} I_\mathrm{sat} = \frac{c\epsilon_0\Gamma^2\hbar^2}{4\vert \epsilon_q\cdot\mathrm{d}\vert^2} \end{align}

  • getSaturationIntensityIsotropic(ng, lg, jg, fg, ne, fe, je, fe)

Returns saturation intenstiy \(I_\mathrm{sat}\) averaged over transiton from \(f_g\rightarrow f_e\) for isotropic polarisation in units of W/m\(^2\).

\begin{align} I_\mathrm{sat,iso} = \frac{c\epsilon_0\Gamma^2\hbar^2}{4\vert \mathrm{d_{iso}}\vert^2}, \end{align} where \(\vert \mathrm{d_{iso}}\vert^2\) is the isotropic dipole moment obtained by averaging over all possible transitions and polarisations. \begin{align} \vert d_{iso} \vert^2 = \frac{1}{3(2fg+1)}\displaystyle\sum_{q,m_{f_g}}\vert\langle ng,lg,jg,fg,mfg|er|ne,le,je,fe,mfg+q\rangle\vert^2 \end{align}

[3]:
# Test sph=atm.getSphericalDipoleMatrixElement(j1,mj1,j2,mj2,q)
# ============================================================
print("\nTest sph=atm.getSphericalDipoleMatrixElement(j1,mj1,j2,mj2,q):")
print("==============================================================")
print("Spherical matrix element in units of reduced matrix element <j'||er||j> for j and  <f'||er||f> for f")
sph = cs.getSphericalDipoleMatrixElement(0.5, 0.5, 1.5, 1.5, 1)
print("|1/2,1/2> to |3/2,3/2> for q=1 (sigma+) : sph = 1/2 = %2.5f, sph = %2.5f"
      % (1 / 2, sph)
     )
sph = cs.getSphericalDipoleMatrixElement(4, 4, 5, 5, 1)
print("|4,4> to |5,5> for q=1 (sigma+) : sph = sqrt{1/11} = %2.5f, sph = %2.5f"
      % (np.sqrt(1 / 11), sph)
     )
sph = cs.getSphericalDipoleMatrixElement(4, 4, 4, 4, 0)
print("|4,4> to |4,4> for q=0 (pi) : sph = sqrt{4/45} = %2.5f, sph = %2.5f"
      % (np.sqrt(4 / 45), sph))
sph = cs.getSphericalDipoleMatrixElement(4, 4, 3, 3, -1)
print("|4,4> to |3,3> for q=-1 (sigma-) : sph = 1/3 = %2.5f, sph = %2.5f"
      % (1 / 3, sph)
     )

# Test sph=atm.getSphericalMatrixElementHFStoFS(j1,f1,mf1,j2,mj2,q)
# =================================================================
print("\nTest sph=atm.getSphericalMatrixElementHFStoFS(j1,f1,mf1,j2,mj2,q):")
print("==================================================================")
print("Spherical matrix element in units of reduced matrix element <j'||er||j'>")
sph = cs.getSphericalMatrixElementHFStoFS(1.5, 5, 5, 0.5, 0.5, -1)
print("Cs P3/2 |5,5> to S1/2 |1/2,1/2> for q=-1 (sigma-) : sph = 1/2 = %2.5f, sph = %2.5f"
       % (1 / 2, sph)
     )
sph = cs.getSphericalMatrixElementHFStoFS(1.5, 5, 1, 0.5, 0.5, 1)
print("Cs P3/2 |5,1> to S1/2 |1/2,1/2> for q=1 (sigma+) : sph = sqrt(1/40) = %2.5f, sph = %2.5f"
       % (np.sqrt(1 / 40), sph)
     )
sph = cs.getSphericalMatrixElementHFStoFS(1.5, 5, 1, 0.5, 0.5, 0)
print("Cs P3/2 |5,1> to S1/2 |1/2,1/2> for q=0 (pi) : sph = sqrt(1/12) = %2.5f, sph = %2.5f"
       % (np.sqrt(1 / 12), sph)
     )
sph = cs.getSphericalMatrixElementHFStoFS(1.5, 5, 1, 0.5, 0.5, -1)
print("Cs P3/2 |5,1> to S1/2 |1/2,1/2> for q=-1 (sigma-) : sph = sqrt(1/24) = %2.5f, sph = %2.5f"
      % (np.sqrt(1 / 24), sph))

# Test sph=atm.getDipoleMatrixElementHFStoFS(j1,f1,mf1,j2,mj2,q)
# ==============================================================
print("\nTest rme=atm.getDipoleMatrixElementHFStoFS(n1,l1,j1,f1,mf1,n2,l2,j2,mj2,q):")
print("=============================================================================")
print("Spherical matrix element in ea0")
rmeJ = cs.getReducedMatrixElementJ(6, 1, 1.5, 80, 0, 0.5)
print("Cs <6P3/2||er||80S1/2> = 0.0086 = %f" % (rmeJ))
rme = cs.getDipoleMatrixElementHFStoFS(6, 1, 1.5, 5, 5, 80, 0, 0.5, 0.5, -1)
print("Cs 6P3/2 |5,5> to 80S1/2 |1/2,1/2> for q=-1 (sigma-) : 0.0043 ea0, rme = %f ea0"
      % (rme))
rme = cs.getDipoleMatrixElementHFStoFS(6, 1, 1.5, 5, 1, 80, 0, 0.5, 0.5, -1)
print("Cs 6P3/2 |5,1> to 80S1/2 |1/2,1/2> for q=-1 (sigma-) : 0.0018 ea0, rme = %f ea0"
      % (rme))

rmeJ = cs.getReducedMatrixElementJ(6, 1, 1.5, 80, 2, 2.5)
print("Cs <6P3/2||er||80D5/2> = -0.0259 = %f"
      % (rmeJ))
rme = cs.getDipoleMatrixElementHFStoFS(6, 1, 1.5, 5, 5, 80, 2, 2.5, 2.5, 1)
print("Cs 6P3/2 |5,5> to 80D5/2 |5/2,5/2> for q=+1 (sigma+) : -0.0106 ea0, rme = %f ea0"
      % (rme))

sph = cs.getSphericalMatrixElementHFStoFS(1.5, 5, 5, 2.5, 2.5, 1)
print("Cs P3/2 |5,5> to D5/2 |5/2,5/2> for q=+1 (sigma+) : sph = sqrt(1/6) = %2.5f, sph = %2.5f"
       % (np.sqrt(1 / 6), sph))


# Test sph=atm.getMagneticDipoleMatrixElementHFS(l,j,f1,mf1,f2,mf2,q)
#==================================================================
print("\nTest sphB=atm.getMagneticDipoleMatrixElementHFS(l,j,f1,mf1,f2,mf2,q):")
print("=====================================================================")
print("Spherical matrix element in uB*B")
sphB = cs.getMagneticDipoleMatrixElementHFS(0, 0.5, 4, 0, 3, 0, 0)
print("Cs 6S1/2 <4,0|B_0|3,0> = 1.000 = %f muB"
      % (sphB))
sphB = cs.getMagneticDipoleMatrixElementHFS(0, 0.5, 4, 3, 3, 3, 0)
print("Cs 6S1/2 <4,3|B_0|3,3> = sqrt(7/16) = %f = %f muB"
      % (np.sqrt(7 / 16), sphB))
sphB = cs.getMagneticDipoleMatrixElementHFS(0, 0.5, 4, 4, 3, 3, -1)
print("Cs 6S1/2 <4,4|B_1|3,3> = sqrt(7/4) = %f = %f muB"
      % (np.sqrt(7 / 4), sphB))
sphB = rb87.getMagneticDipoleMatrixElementHFS(0, 0.5, 2, 0, 1, 0, 0)
print("Rb 5S1/2 <2,0|B_0|1,0> = 1.000 = %f muB"
       % (sphB))
sphB = rb87.getMagneticDipoleMatrixElementHFS(0, 0.5, 2, 1, 1, 1, 0)
print("Rb 5S1/2 <2,1|B_0|1,1> = sqrt(3/4) = %f = %f muB"
      % (np.sqrt(3 / 4), sphB))
sphB = rb87.getMagneticDipoleMatrixElementHFS(0, 0.5, 2, 2, 1, 1, -1)
print("Rb 5S1/2 <2,2|B_1|1,1> = sqrt(3/2) = %f = %f muB"
      % (np.sqrt(3 / 2), sphB))

# Test getBranchingRatio(jg,fg,mfg,je,fe,mfe)
# ===========================================
print("\nTest br=atm.getBranchingRatio(jg,fg,mfg,je,fe,mfe):")
print("===================================================")
# Cs P3/2 |5,5> to S1/2|4,4> = 1
br = cs.getBranchingRatio(0.5, 4, 4, 1.5, 5, 5)
print("Cs P3/2 |5,5> to S1/2|4,4> = 1, br = %2.5f"
      % (br))
# Cs P3/2 |5,4> to S1/2|4,4> = 1/5
br = cs.getBranchingRatio(0.5, 4, 4, 1.5, 5, 4)
print("Cs P3/2 |5,4> to S1/2|4,4> = 1/5 = %2.5f, br = %2.5f"
      % (1 / 5, br))
# Cs P3/2 |5,4> to S1/2|4,3> = 4/5
br = cs.getBranchingRatio(0.5, 4, 3, 1.5, 5, 4)
print("Cs P3/2 |5,4> to S1/2|4,3> = 4/5 = %2.5f, br = %2.5f"
      % (4 / 5, br))
# Cs P3/2 |4,4> to S1/2|4,4> = 7/15
br = cs.getBranchingRatio(0.5, 4, 4, 1.5, 4, 4)
print("Cs P3/2 |4,4> to S1/2|4,4> = 7/15 = %2.5f, br = %2.5f"
      % (7 / 15, br))
# Cs P3/2 |4,4> to S1/2|4,3> = 7/60
br = cs.getBranchingRatio(0.5, 4, 3, 1.5, 4, 4)
print("Cs P3/2 |4,4> to S1/2|4,3> = 7/60 = %2.5f, br = %2.5f"
      % (7 / 60, br))
# Cs P3/2 |4,4> to S1/2|3,3> = 7/60
br = cs.getBranchingRatio(0.5, 3, 3, 1.5, 4, 4)
print("Cs P3/2 |4,4> to S1/2|3,3> = 5/12 = %2.5f, br = %2.5f"
      % (5 / 12, br))


# Test getSaturationIntensity(ng,lg,jg,fg,mfg,ne,le,je,fe,mfe)
#============================================================
print("\nTest Is=atm.getSaturationIntensity(ng,lg,jg,fg,mfg,ne,le,je,fe,mfe):")
print("====================================================================")
# Cs P3/2 |5,5> to S1/2|4,4> = 11 W/m^2
Is = cs.getSaturationIntensity(6, 0, 0.5, 4, 4, 6, 1, 1.5, 5, 5)
print("Cs 6S1/2|4,4> to 6P3/2 |5,5> Is = 11 W/m^2, br = %2.5f W/m^2"
      %(Is))
Is = cs.getSaturationIntensity(6, 0, 0.5, 4, 0, 6, 1, 1.5, 5, 0)
print("Cs 6S1/2|4,0> to 6P3/2 |5,0> Is = 19.8 W/m^2, br = %2.5f W/m^2"
      %(Is))

# Test getSaturationIntensityIsotropic(ng,lg,jg,fg,ne,le,je,fe)
# ============================================================
print("\nTest Is=atm.getSaturationIntensityIsotropic(ng,lg,jg,fg,ne,le,je,fe):")
print("======================================================================")
# Cs P3/2 |5,5> to S1/2|4,4> = 11 W/m^2
Is = cs.getSaturationIntensityIsotropic(6, 0, 0.5, 4, 6, 1, 1.5, 5)
print("Cs 6S1/2 F=4 to 6P3/2 F=5 Is_ISO = 27.059 W/m^2 [1], Is_ISO = %2.5f W/m^2"
      %(Is))
Is = cs.getSaturationIntensityIsotropic(6, 0, 0.5, 4, 6, 1, 1.5, 4)
print("Cs 6S1/2 F=4 to 6P3/2 F=4 Is_ISO = 56.69 W/m^2, Is_ISO = %2.5f W/m^2"
      %(Is))


Test sph=atm.getSphericalDipoleMatrixElement(j1,mj1,j2,mj2,q):
==============================================================
Spherical matrix element in units of reduced matrix element <j'||er||j> for j and  <f'||er||f> for f
|1/2,1/2> to |3/2,3/2> for q=1 (sigma+) : sph = 1/2 = 0.50000, sph = 0.50000
|4,4> to |5,5> for q=1 (sigma+) : sph = sqrt{1/11} = 0.30151, sph = 0.30151
|4,4> to |4,4> for q=0 (pi) : sph = sqrt{4/45} = 0.29814, sph = 0.29814
|4,4> to |3,3> for q=-1 (sigma-) : sph = 1/3 = 0.33333, sph = 0.33333

Test sph=atm.getSphericalMatrixElementHFStoFS(j1,f1,mf1,j2,mj2,q):
==================================================================
Spherical matrix element in units of reduced matrix element <j'||er||j'>
Cs P3/2 |5,5> to S1/2 |1/2,1/2> for q=-1 (sigma-) : sph = 1/2 = 0.50000, sph = 0.50000
Cs P3/2 |5,1> to S1/2 |1/2,1/2> for q=1 (sigma+) : sph = sqrt(1/40) = 0.15811, sph = 0.15811
Cs P3/2 |5,1> to S1/2 |1/2,1/2> for q=0 (pi) : sph = sqrt(1/12) = 0.28868, sph = 0.28868
Cs P3/2 |5,1> to S1/2 |1/2,1/2> for q=-1 (sigma-) : sph = sqrt(1/24) = 0.20412, sph = 0.20412

Test rme=atm.getDipoleMatrixElementHFStoFS(n1,l1,j1,f1,mf1,n2,l2,j2,mj2,q):
=============================================================================
Spherical matrix element in ea0
Cs <6P3/2||er||80S1/2> = 0.0086 = -0.008715
Cs 6P3/2 |5,5> to 80S1/2 |1/2,1/2> for q=-1 (sigma-) : 0.0043 ea0, rme = -0.004358 ea0
Cs 6P3/2 |5,1> to 80S1/2 |1/2,1/2> for q=-1 (sigma-) : 0.0018 ea0, rme = -0.001779 ea0
Cs <6P3/2||er||80D5/2> = -0.0259 = -0.025645
Cs 6P3/2 |5,5> to 80D5/2 |5/2,5/2> for q=+1 (sigma+) : -0.0106 ea0, rme = -0.010469 ea0
Cs P3/2 |5,5> to D5/2 |5/2,5/2> for q=+1 (sigma+) : sph = sqrt(1/6) = 0.40825, sph = 0.40825

Test sphB=atm.getMagneticDipoleMatrixElementHFS(l,j,f1,mf1,f2,mf2,q):
=====================================================================
Spherical matrix element in uB*B
Cs 6S1/2 <4,0|B_0|3,0> = 1.000 = 1.000000 muB
Cs 6S1/2 <4,3|B_0|3,3> = sqrt(7/16) = 0.661438 = 0.661438 muB
Cs 6S1/2 <4,4|B_1|3,3> = sqrt(7/4) = 1.322876 = 1.322876 muB
Rb 5S1/2 <2,0|B_0|1,0> = 1.000 = 1.000000 muB
Rb 5S1/2 <2,1|B_0|1,1> = sqrt(3/4) = 0.866025 = 0.866025 muB
Rb 5S1/2 <2,2|B_1|1,1> = sqrt(3/2) = 1.224745 = 1.224745 muB

Test br=atm.getBranchingRatio(jg,fg,mfg,je,fe,mfe):
===================================================
Cs P3/2 |5,5> to S1/2|4,4> = 1, br = 1.00000
Cs P3/2 |5,4> to S1/2|4,4> = 1/5 = 0.20000, br = 0.20000
Cs P3/2 |5,4> to S1/2|4,3> = 4/5 = 0.80000, br = 0.80000
Cs P3/2 |4,4> to S1/2|4,4> = 7/15 = 0.46667, br = 0.46667
Cs P3/2 |4,4> to S1/2|4,3> = 7/60 = 0.11667, br = 0.11667
Cs P3/2 |4,4> to S1/2|3,3> = 5/12 = 0.41667, br = 0.41667

Test Is=atm.getSaturationIntensity(ng,lg,jg,fg,mfg,ne,le,je,fe,mfe):
====================================================================
Cs 6S1/2|4,4> to 6P3/2 |5,5> Is = 11 W/m^2, br = 11.02370 W/m^2
Cs 6S1/2|4,0> to 6P3/2 |5,0> Is = 19.8 W/m^2, br = 19.84266 W/m^2

Test Is=atm.getSaturationIntensityIsotropic(ng,lg,jg,fg,ne,le,je,fe):
======================================================================
Cs 6S1/2 F=4 to 6P3/2 F=5 Is_ISO = 27.059 W/m^2 [1], Is_ISO = 27.05817 W/m^2
Cs 6S1/2 F=4 to 6P3/2 F=4 Is_ISO = 56.69 W/m^2, Is_ISO = 56.69332 W/m^2

Magnetic Dipole Matrix Element#

  • getMagneticDipoleMatrixElementHFS(l, j, f1, mf1, m2, mf2, q)

Returns spherical component of magnetic dipole transition from \(\vert f_1,m_{f_1}\rangle\rightarrow\vert f_2,m_{f_2}\rangle\) in units of \(\mu_B B_q\) corresponding to relative strength of microwave transitions between hyperfine states of the same \(n\ell j\) state driven by magnetic field with polarisation \(q\).

\begin{align} \langle f_1,m_{f_1}\vert \mu_q \vert f_2,m_{f_2}\rangle = g_J \mu_B B_q (-1)^{f_2+j+I+1+f_1-m_{f_1}}\sqrt{(2f_1+1)(2f_2+1)j(j+1)(2j+1)} \begin{pmatrix}f_1&1&f_2\\-m_{f_1} & -q & m_{f_2}\end{pmatrix}\begin{Bmatrix}f_1&1&f_2\\j & I & j\end{Bmatrix} \end{align}

Ground-State Raman Transitions#

0c3f702ff72a4325a1dc1ef10e406b29

We consider here the case of two-photon ground-state Raman transitions from \(\vert 0 \rangle = \vert f_0,m_{f_0}\rangle\) to \(\vert 1 \rangle = \vert f_1,m_{f_1}\rangle\) via an intermediate excited state \(\vert e\rangle = \vert f_e,m_{f_e}\rangle\). In the rotating wave approximation, the Hamiltonian can be written as

\begin{equation} \mathcal{H} = \hbar\begin{pmatrix} 0 & \Omega^a_{0\rightarrow f_e}/2 & 0 \\\Omega^a_{0\rightarrow f_e}/2 & -\Delta & \Omega^b_{1\rightarrow f_e}/2 \\ 0 & \Omega^b_{1\rightarrow f_e}/2 & -\delta\end{pmatrix}, \end{equation}

where \(\Omega^a_{0\rightarrow f_e}\) is the Rabi frequency from \(\vert 0\rangle\rightarrow\vert e\rangle\) due to laser A and \(\Omega^b_{1\rightarrow f_e}\) is the Rabi frequency from \(\vert 1\rangle\rightarrow\vert e\rangle\) due to laser B, \(\Delta=(\Delta_a+\Delta_a)/2\) is the average detuning from the intermediate state and \(\delta=\Delta_a-\Delta_b\) is the two-photon detuning.

In the limit \(\Delta \gg \Omega^a,\Omega^b,\Gamma_e\) the intermediate excited state can be adiabatically elliminiated leading to an effective two-photon Rabi frequency \(\Omega_R = \Omega^a\Omega^b/2\Delta\) and an effective detuning \(\delta_\mathrm{eff}=\delta+\Delta_\mathrm{AC}\) where \(\Delta_\mathrm{AC}\) is the differential Stark shift given by

\begin{equation} \Delta_\mathrm{AC} = \frac{\vert\Omega^b_{0\rightarrow f_e}\vert^2}{4\Delta} - \frac{\vert\Omega^a_{1\rightarrow f_e}\vert^2}{4\Delta}. \end{equation}

On the two-photon transition (\(\delta_\mathrm{eff}=0\)), the probability of spontanteous emission during a \(\pi\)-pulse can be evaluated from \(P_\mathrm{sc}=\Gamma_e\rho_e\tau_\pi\), where \(\tau_\pi=\pi/\Omega_R\) is the \(\pi\)-pulse duration and \(\rho_e\) is the excited state population given by

\begin{equation} \rho_e = \frac{1}{2}\left(\frac{\vert\Omega^a_{0\rightarrow f_e}\vert^2}{2\Delta^2}+\frac{\vert\Omega^b_{1\rightarrow f_e}\vert^2}{2\Delta^2}\right). \end{equation}

For a multi-level atom, the terms describing Raman transition, diffrential AC Stark shift and scattering probability must be calculated by summing over the \(\vert f_e,m_{f_e}\rangle\) states of the intermediate excited state. Additionally, for transitions between hyperfine ground states separated by frequency \(\omega_{01}\) we must also account for the off-resonant AC Stark shift and scattering terms arising from laser A acting on \(\vert 1 \rangle\) with detuning \(\Delta'=\Delta+\omega_{01}\) and laser B acting on \(\vert 0 \rangle\) with detuning \(\Delta'=\Delta-\omega_{01}\).

\begin{align} \Omega_R&=\displaystyle\sum_{f_e,m_{f_e}}\frac{\Omega^a_{0\rightarrow f_e}\Omega^b_{1\rightarrow f_e}}{2(\Delta-\Delta_{f_e})}\\ \Delta_{\mathrm{AC}} &= \displaystyle\sum_{f_e,m_{f_e}}\left[\frac{\vert\Omega^a_{0\rightarrow f_e}\vert^2-\vert\Omega^b_{1\rightarrow f_e}\vert^2}{4(\Delta-\Delta_{f_e})}+\frac{\vert\Omega^a_{1\rightarrow f_e}\vert^2}{4(\Delta+\omega_{01}-\Delta_{f_e})}-\frac{\vert\Omega^b_{0\rightarrow f_e}\vert^2}{4(\Delta-\omega_{01}-\Delta_{f_e})}\right],\\ P_\mathrm{sc} &= \frac{\Gamma_et_\pi}{2}\displaystyle\sum_{f_e,m_{f_e}}\left[\frac{\vert\Omega^a_{0\rightarrow f_e}\vert^2}{2(\Delta-\Delta_{f_e})^2}+\frac{\vert\Omega^b_{1\rightarrow f_e}\vert^2}{2(\Delta-\Delta_{f_e})^2}+\frac{\vert\Omega^a_{1\rightarrow f_e}\vert^2}{4(\Delta+\omega_{01}-\Delta_{f_e})^2}+\frac{\vert\Omega^b_{0\rightarrow f_e}\vert^2}{4(\Delta-\omega_{01}-\Delta_{f_e})^2}\right], \end{align}

  • OmegaR, AC, Psc = groundStateRamanTransition(Pa, wa, qa, Pb, wb, qb, Delta, f0, mf0, f1, mf1, ne, le, je)

Returns two-photon Rabi frequency \(\Omega_R\), differential AC Starkshift \(\Delta_\mathrm{AC}\) and scattering probability \(P_\mathrm{sc}\) for transition from \(\vert f_0,m_{f_0}\rangle\rightarrow\vert f_1,m_{f_1}\rangle\) via intermediate excited state with quantum numbers \(n_e,\ell_e,j_e\) with detuning \(\Delta\) from excited state centre of mass using laser \(a/b\) with power \(P\), waist \(w\) and polarisation \(q\) driving transitions from \(f_0,f_1\) respectively.

[4]:
# Test Raman Transition Calculation
#==================================

#Laser Parameters
Pa = 0.1e-3; wa = 5e-6; qa = -1;
Pb = 0.1e-3; wb = 5e-6; qb = -1;

#Detuning Array
Delta0 = -50e9 * 2 * np.pi
Delta = np.linspace(-150, 150, 10000) * 1e9 * 2 * np.pi

#D1 Line
ne = 6; le = 1; je = 0.5; f1 = 4; mf1 = 0; f0 = 3; mf0 = 0;
OmR, AC, Psc = cs.groundStateRamanTransition(Pa, wa, qa,
                                             Pb, wb, qb,
                                             Delta0,
                                             f0, mf0, f1, mf1,
                                             ne, le, je)
print('\nRaman Transitions |%d,%d> to |%d,%d> via %dP_%d/2'
      % (f1, mf1, f0, mf0, ne, je*2))
print('\tParameters:\tDelta/2pi = %2.f GHz, Pa = %1.2f mW, wa = %2.1f um, Pa = %1.2f mW, wa = %2.1f um' \
      % (Delta0 / 2.0 / np.pi * 1e-9, Pa * 1e3, wa * 1e6, Pb * 1e3, wb * 1e6))
print('\tResults:\tOmR/2pi = %2.5f MHz, AC/2pi = %2.5f MHz, Psc = %f'
      % (OmR / 2.0 / np.pi * 1e-6, AC / 2.0 / np.pi * 1e-6, Psc))

#D2 Line
ne = 6; le = 1; je = 1.5;
OmR, AC, Psc = cs.groundStateRamanTransition(Pa, wa, qa,
                                             Pb, wb, qb,
                                             Delta0,
                                             f0, mf0,
                                             f1, mf1,
                                             ne,le,je)
print('\nRaman Transitions |%d,%d> to |%d,%d> via %dP_%d/2' % (f1, mf1, f0, mf0, ne, je*2))
print('\tParameters:\tDelta/2pi = %2.f GHz, Pa = %1.2f mW, wa = %2.1f um, Pa = %1.2f mW, wa = %2.1f um' \
      % (Delta0 / 2.0 / np.pi * 1e-9, Pa * 1e3, wa * 1e6, Pb * 1e3, wb * 1e6))
print('\tResults:\tOmR/2pi = %2.5f MHz, AC/2pi = %2.5f MHz, Psc = %f'
      % (OmR / 2.0 / np.pi * 1e-6, AC / 2.0 / np.pi * 1e-6, Psc))

#7P1/2
ne = 7; le = 1; je = 0.5; Pa = 10e-3; Pb = 10e-3;
OmR, AC, Psc = cs.groundStateRamanTransition(Pa, wa, qa,
                                             Pb, wb, qb,
                                             Delta0,
                                             f0, mf0,
                                             f1, mf1,
                                             ne, le, je)
print('\nRaman Transitions |%d,%d> to |%d,%d> via %dP_%d/2'
       % (f1, mf1, f0, mf0, ne, je*2))
print('\tParameters:\tDelta/2pi = %2.f GHz, Pa = %1.2f mW, wa = %2.1f um, Pa = %1.2f mW, wa = %2.1f um' \
      % (Delta0 / 2.0 / np.pi * 1e-9, Pa * 1e3, wa * 1e6, Pb * 1e3, wb * 1e6)
     )
print('\tResults:\tOmR/2pi = %2.5f MHz, AC/2pi = %2.5f MHz, Psc = %f'
      % (OmR / 2.0 / np.pi * 1e-6, AC / 2.0 / np.pi * 1e-6, Psc))

#7P3/2
ne = 7; le = 1; je = 1.5; Pa = 10e-3; Pb = 10e-3;
OmR, AC, Psc = cs.groundStateRamanTransition(Pa, wa, qa,
                                             Pb, wb, qb,
                                             Delta0,
                                             f0, mf0, f1,
                                             mf1, ne, le, je)
print('\nRaman Transitions |%d,%d> to |%d,%d> via %dP_%d/2'
      % (f1, mf1, f0, mf0, ne, je * 2))
print('\tParameters:\tDelta/2pi = %2.f GHz, Pa = %1.2f mW, wa = %2.1f um, Pa = %1.2f mW, wa = %2.1f um' \
      % (Delta0 / 2.0 / np.pi * 1e-9, Pa * 1e3, wa * 1e6, Pb * 1e3, wb * 1e6))
print('\tResults:\tOmR/2pi = %2.5f MHz, AC/2pi = %2.5f MHz, Psc = %f'
      % (OmR / 2.0 / np.pi * 1e-6, AC / 2.0 / np.pi * 1e-6, Psc))

#Prepare Output Figure
f, ax = plt.subplots(1, 3, figsize=(12, 5))
ax[0].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[0].set_ylabel(r'$\Omega/2\pi$ (MHz)')
ax[0].set_ylim([-20, 20])
ax[1].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[1].set_ylabel(r'$\Delta_\mathrm{AC}/2\pi$ (MHz)')
ax[1].set_ylim([-20, 20])
ax[2].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[2].set_ylabel(r'$P_\mathrm{sc}$')
ax[2].set_ylim([0, 0.01])

#Reset to D2 Line
ne = 6; le = 1; je = 1.5; f1 = 4; mf1 = 0; f0 = 3; mf0 = 0;
Pa = 0.1e-3; wa = 5e-6; qa = -1;
Pb = 0.1e-3; wb = 5e-6; qb = -1;

for mf in range(-3,4):
    OmR, AC, Psc = cs.groundStateRamanTransition(Pa, wa, qa,
                                                 Pb, wb, qb,
                                                 Delta,
                                                 f1, mf, f0, mf + qa - qb,
                                                 ne, le, je)
    ax[0].plot(Delta / (2e9 * np.pi), OmR / (2e6 * np.pi))
    ax[1].plot(Delta / (2e9 * np.pi), AC / (2e6 * np.pi))
    ax[2].plot(Delta / (2e9 * np.pi), Psc,
               label=format("$\\vert%d,%d\\rangle\\rightarrow\\vert%d,%d\\rangle$"
                            % (f1,mf,f0,mf)))
ax[2].legend(loc=4)
ax[1].set_title('D2 Raman Transitions')
plt.show()

Raman Transitions |4,0> to |3,0> via 6P_1/2
        Parameters:     Delta/2pi = -50 GHz, Pa = 0.10 mW, wa = 5.0 um, Pa = 0.10 mW, wa = 5.0 um
        Results:        OmR/2pi = -10.58926 MHz, AC/2pi = 2.01241 MHz, Psc = 0.000603

Raman Transitions |4,0> to |3,0> via 6P_3/2
        Parameters:     Delta/2pi = -50 GHz, Pa = 0.10 mW, wa = 5.0 um, Pa = 0.10 mW, wa = 5.0 um
        Results:        OmR/2pi = 10.51457 MHz, AC/2pi = 3.86727 MHz, Psc = 0.001380

Raman Transitions |4,0> to |3,0> via 7P_1/2
        Parameters:     Delta/2pi = -50 GHz, Pa = 10.00 mW, wa = 5.0 um, Pa = 10.00 mW, wa = 5.0 um
        Results:        OmR/2pi = -3.97773 MHz, AC/2pi = 0.75654 MHz, Psc = 0.000136

Raman Transitions |4,0> to |3,0> via 7P_3/2
        Parameters:     Delta/2pi = -50 GHz, Pa = 10.00 mW, wa = 5.0 um, Pa = 10.00 mW, wa = 5.0 um
        Results:        OmR/2pi = 8.98075 MHz, AC/2pi = 3.37987 MHz, Psc = 0.000325
_images/ARC_3_1_additions_8_1.png

Two-photon Rydberg excitation#

3b51fadc003e4d0d9d17f51946e180b0

We consider here the case of two-photon ground-state Raman transitions from \(\vert g \rangle = \vert f_g,m_{f_g}\rangle\) to \(\vert r \rangle = \vert n_r,l_r,j_r,m_{j_r}\rangle\) via an intermediate excited state \(\vert e\rangle = \vert f_e,m_{f_e}\rangle\). In the rotating wave approximation, the Hamiltonian can be written as

\begin{equation} \mathcal{H} = \hbar\begin{pmatrix} 0 & \Omega_p/2 & 0 \\\Omega_p/2 & -\Delta_p & \Omega_c/2 \\ 0 & \Omega_c/2 & -\delta\end{pmatrix}, \end{equation}

where \(\Omega_p\) is the Rabi frequency from \(\vert 0\rangle\rightarrow\vert e\rangle\) due to the probe laser and \(\Omega_c\) is the Rabi frequency from \(\vert e\rangle\rightarrow\vert r\rangle\) due to the coupling laser, \(\delta=(\Delta_p+\Delta_c)\) is the two-photon detuning.

In the limit \(\Delta_p \gg \Omega_p,\Omega_c,\Gamma_e\) the intermediate excited state can be adiabatically elliminiated leading to an effective two-photon Rabi frequency \(\Omega_R = \Omega_p\Omega_c/2\Delta\) and an effective detuning \(\delta_\mathrm{eff}=\delta+\Delta_\mathrm{AC}\) where \(\Delta_\mathrm{AC}\) is the differential Stark shift given by

\begin{equation} \Delta_\mathrm{AC} = \frac{\vert\Omega_p\vert^2}{4\Delta_p} - \frac{\vert\Omega_c\vert^2}{4\Delta_c}. \end{equation}

On the two-photon transition (\(\delta_\mathrm{eff}=0\)), the probability of spontanteous emission during a \(\pi\)-pulse can be evaluated from \(P_\mathrm{sc}=\Gamma_e\rho_e\tau_\pi\), where \(\tau_\pi=\pi/\Omega_R\) is the \(\pi\)-pulse duration and \(\rho_e\) is the excited state population given by

\begin{equation} \rho_e = \frac{1}{2}\left(\frac{\vert\Omega_p\vert^2}{2\Delta_p^2}+\frac{\vert\Omega_c\vert^2}{2\Delta_c^2}\right). \end{equation}

For a multi-level atom, the terms describing Raman transition, diffrential AC Stark shift and scattering probability must be calculated by summing over the \(\vert f_e,m_{f_e}\rangle\) states of the intermediate excited state. Additionally, for transitions between hyperfine ground states separated by frequency \(\omega_{01}\) we must also account for the off-resonant AC Stark shift and scattering terms arising from laser A acting on \(\vert 1 \rangle\) with detuning \(\Delta'=\Delta+\omega_{01}\) and laser B acting on \(\vert 0 \rangle\) with detuning \(\Delta'=\Delta-\omega_{01}\).

\begin{align} \Omega_R&=\displaystyle\sum_{f_e,m_{f_e}}\frac{\Omega_p^{g\rightarrow f_e}\Omega_c^{f_e\rightarrow r}}{2(\Delta-\Delta_{f_e})}\\ \Delta_{\mathrm{AC}_g} &= \displaystyle\sum_{f_e,m_{f_e}}\frac{\vert\Omega_p^{g\rightarrow f_e}\vert^2}{4(\Delta-\Delta_{f_e})}\\ \Delta_{\mathrm{AC}_r} &= \displaystyle\sum_{f_e,m_{f_e}}\frac{\vert\Omega_c^{f_e\rightarrow r}\vert^2}{4(\Delta-\Delta_{f_e})}\\ P_\mathrm{sc} &= \frac{\Gamma_et_\pi}{2}\displaystyle\sum_{f_e,m_{f_e}}\left[\frac{\vert\Omega_p^{g\rightarrow f_e}\vert^2}{2(\Delta-\Delta_{f_e})^2}+\frac{\vert\Omega_c^{f_e\rightarrow r}\vert^2}{2(\Delta-\Delta_{f_e})^2}\right] \end{align}

  • OmegaR, ACg, ACr, Psc = twoPhotonRydbergExcitation(Pp, wp, qp, Pc, wc, qc, Delta, fg, mfg, ne, le, je, nr, lr, jr, mjr)

Returns two-photon Rabi frequency \(\Omega_R\), ground-state AC Starkshift \(\Delta_\mathrm{AC_g}\), Rydberg-state AC Starkshift \(\Delta_\mathrm{AC_r}\) and scattering probability \(P_\mathrm{sc}\) for transition from \(\vert f_g,m_{f_g}\rangle\rightarrow n_rL_{j_r}\vert j_r,m_{j_r}\rangle\) via intermediate excited state with quantum numbers \(n_e,\ell_e,j_e\) with detuning \(\Delta\) from excited state centre of mass driven by probe/coupling lasers with power \(P\), waist \(w\) and polarisation \(q\) driving transitions from \(\vert g\rangle,\vert r\rangle\) respectively.

[5]:
#Rydberg Excitation vs 6P3/2
Pp=500e-9;wp=3e-6;qp=1;
Pc=200e-3;wc=5e-6;qc=-1;
#Ground State
fg=4;mfg=0;
#Excited state (6P3/2)
ne=6;le=1;je=1.5;
#Rydberg state (82S1/2)
nr=80;lr=0;jr=0.5;mjr=-0.5;

SPD='SPD'

#Detuning Array
Delta0=1.1e9*2*np.pi
Delta=np.linspace(-2.5,2.5,1001)*1e9*2*np.pi

[OmR,ACg,ACr,Psc]=cs.twoPhotonRydbergExcitation(Pp,wp,qp,Pc,wc,qc,Delta0,fg,mfg,ne,le,je,nr,lr,jr,mjr)
print('\nRydberg Excitation 6S1/2 |%d,%d> to %d%c |%d/2,%d/2> via %dP_%d/2' % (fg,mfg,nr,SPD[lr],2*jr,2*mjr,ne,je*2))
print('\n\tParameters:\tDelta/2pi = %2.2f GHz, Pp = %1.2f uW, wp = %2.1f um, qp=%d, Pc = %1.2f mW, wc = %2.1f um, qc=%d' \
      % (Delta0/2.0/np.pi*1e-9,Pp*1e6,wp*1e6,qp,Pc*1e3,wc*1e6,qc))
print('\n\tResults:\tOmR/2pi = %2.5f MHz\n\t\t\tACg/2pi = %2.5f MHz\n\t\t\tACr/2pi = %2.5f MHz\n\t\t\tDelta_AC/2pi = %2.5f MHz\n\t\t\tPsc = %f' \
      % (OmR/2.0/np.pi*1e-6,ACg/2.0/np.pi*1e-6,ACr/2.0/np.pi*1e-6,(ACg-ACr)/2.0/np.pi*1e-6,Psc))

#Prepare Output Figure
[f,ax]=plt.subplots(1,3,figsize=(16, 5))
ax[0].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[0].set_ylabel(r'$\Omega/2\pi$ (MHz)')
ax[0].set_ylim([-10,10])
ax[1].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[1].set_ylabel(r'$\Delta_\mathrm{AC}/2\pi$ (MHz)')
ax[1].set_ylim([-20,20])
ax[2].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[2].set_ylabel(r'$P_\mathrm{sc}$')
ax[2].set_ylim([0,0.25])
ax[0].plot(Delta/(2e9*np.pi),0*Delta,'k--')
ax[1].plot(Delta/(2e9*np.pi),0*Delta,'k--')

mj=-jr-1
for m in range(int(2*jr+1)):
    mj+=1
    for mfg in [-fg,0,fg]:
        [OmR,ACg,ACr,Psc]=cs.twoPhotonRydbergExcitation(Pp,wp,qp,Pc,wc,qc,Delta,fg,mfg,ne,le,je,nr,lr,jr,mj)
        if(np.sum(np.abs(OmR))>0.1):
            ax[0].plot(Delta/(2e9*np.pi),OmR/(2e6*np.pi))
            ax[1].plot(Delta/(2e9*np.pi),(ACg-ACr)/(2e6*np.pi))
            ax[2].plot(Delta/(2e9*np.pi),Psc,label=format("$\\vert%d,%d\\rangle\\rightarrow%d%c\\vert%1.0f/2,%1.0f/2\\rangle$" \
                                                      % (fg,mfg,nr,SPD[lr],2*jr,2*mj)))
ax[2].legend(loc=1)
ax[1].set_title('Rydberg excitation 6S_{1/2} -> %dP_{%1.f/2} -> %d%c_{%1.0f/2}'%(ne,2*je,nr,SPD[lr],2*jr))


plt.show()

Rydberg Excitation 6S1/2 |4,0> to 80S |1/2,-1/2> via 6P_3/2

        Parameters:     Delta/2pi = 1.10 GHz, Pp = 0.50 uW, wp = 3.0 um, qp=1, Pc = 200.00 mW, wc = 5.0 um, qc=-1

        Results:        OmR/2pi = -3.31688 MHz
                        ACg/2pi = 7.64903 MHz
                        ACr/2pi = 7.54096 MHz
                        Delta_AC/2pi = 0.10807 MHz
                        Psc = 0.077362
/home/nikola/Documents/Desktop/github/newARC3/ARC-Alkali-Rydberg-Calculator/doc/../arc/alkali_atom_functions.py:2616: RuntimeWarning: divide by zero encountered in true_divide
  tau_pi = pi / abs(OmegaR)
_images/ARC_3_1_additions_10_2.png
[6]:
#Rydberg Excitation vs 7P1/2
Pp=1.535e-3;wp=5e-6;qp=1;
Pc=20e-3;wc=3e-6;qc=-1;
#Ground State
fg=4;mfg=0;
#Excited state (7P1/2)
ne=7;le=1;je=0.5;
#Rydberg state (82S1/2)
nr=80;lr=0;jr=0.5;mjr=-0.5;

#Detuning Array
Delta0=1e9*2*np.pi
Delta=np.linspace(-2.5,2.5,1001)*1e9*2*np.pi

[OmR,ACg,ACr,Psc]=cs.twoPhotonRydbergExcitation(Pp,wp,qp,Pc,wc,qc,Delta0,fg,mfg,ne,le,je,nr,lr,jr,mjr)
print('\nRydberg Excitation 6S1/2 |%d,%d> to %d%c |%d/2,%d/2> via %dP_%d/2' % (fg,mfg,nr,SPD[lr],2*jr,2*mjr,ne,je*2))
print('\n\tParameters:\tDelta/2pi = %2.2f GHz, Pp = %1.2f uW, wp = %2.1f um, qp=%d, Pc = %1.2f mW, wc = %2.1f um, qc=%d' \
      % (Delta0/2.0/np.pi*1e-9,Pp*1e6,wp*1e6,qp,Pc*1e3,wc*1e6,qc))
print('\n\tResults:\tOmR/2pi = %2.5f MHz\n\t\t\tACg/2pi = %2.5f MHz\n\t\t\tACr/2pi = %2.5f MHz\n\t\t\tDelta_AC/2pi = %2.5f MHz\n\t\t\tPsc = %f' \
      % (OmR/2.0/np.pi*1e-6,ACg/2.0/np.pi*1e-6,ACr/2.0/np.pi*1e-6,(ACg-ACr)/2.0/np.pi*1e-6,Psc))

#Prepare Output Figure
[f,ax]=plt.subplots(1,3,figsize=(16, 5))
ax[0].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[0].set_ylabel(r'$\Omega/2\pi$ (MHz)')
ax[0].set_ylim([-10,10])
ax[1].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[1].set_ylabel(r'$\Delta_\mathrm{AC}/2\pi$ (MHz)')
ax[1].set_ylim([-20,20])
ax[2].set_xlabel(r'$\Delta/2\pi$ (GHz)')
ax[2].set_ylabel(r'$P_\mathrm{sc}$')
ax[2].set_ylim([0,0.25])
ax[0].plot(Delta/(2e9*np.pi),0*Delta,'k--')
ax[1].plot(Delta/(2e9*np.pi),0*Delta,'k--')

SPD='SPD'

mj=-jr-1
for m in range(int(2*jr+1)):
    mj+=1
    for mfg in [-fg,0,fg]:
        [OmR,ACg,ACr,Psc]=cs.twoPhotonRydbergExcitation(Pp,wp,qp,Pc,wc,qc,Delta,fg,mfg,ne,le,je,nr,lr,jr,mj)
        if(np.sum(np.abs(OmR))>0.1):
            ax[0].plot(Delta/(2e9*np.pi),OmR/(2e6*np.pi))
            ax[1].plot(Delta/(2e9*np.pi),(ACg-ACr)/(2e6*np.pi))
            ax[2].plot(Delta/(2e9*np.pi),Psc,label=format("$\\vert%d,%d\\rangle\\rightarrow%d%c\\vert%1.0f/2,%1.0f/2\\rangle$" \
                                                      % (fg,mfg,nr,SPD[lr],2*jr,2*mj)))
ax[2].legend(loc=1)
ax[1].set_title('Rydberg excitation 6S_{1/2} -> %dP_{%1.f/2} -> %d%c_{%1.0f/2}'%(ne,2*je,nr,SPD[lr],2*jr))


plt.show()

Rydberg Excitation 6S1/2 |4,0> to 80S |1/2,-1/2> via 7P_1/2

        Parameters:     Delta/2pi = 1.00 GHz, Pp = 1535.00 uW, wp = 5.0 um, qp=1, Pc = 20.00 mW, wc = 3.0 um, qc=-1

        Results:        OmR/2pi = 11.70666 MHz
                        ACg/2pi = 16.15630 MHz
                        ACr/2pi = 16.59130 MHz
                        Delta_AC/2pi = -0.43500 MHz
                        Psc = 0.009768
/home/nikola/Documents/Desktop/github/newARC3/ARC-Alkali-Rydberg-Calculator/doc/../arc/alkali_atom_functions.py:2616: RuntimeWarning: divide by zero encountered in true_divide
  tau_pi = pi / abs(OmegaR)
_images/ARC_3_1_additions_11_2.png

Magnetic Functions#

  • getLandegj(l, j, s=0.5)

    Returns fine-structure Landé g-factor \(g_J\)

\begin{equation} g_J\simeq 1+\frac{j(j+1)+s(s+1)-l(l+1)}{2j(j+1)} \end{equation}

  • getLandegjExact(l, j, s=0.5)

    Returns fine-structure Landé g-factor \(g_J\) using precise \(g_L,\) \(g_S\) values

\begin{equation} g_J=g_L\frac{j(j+1)-s(s+1)+l(l+1)}{2j(j+1)}+g_S\frac{j(j+1)+s(s+1)-l(l+1)}{2j(j+1)} \end{equation}

  • getLandegf(l, j, f, s=0.5)

    Returns hyperfine-structure Landé g-factor \(g_F\)

\begin{equation} g_F\simeq g_J\frac{f(f+1)-I(I+1)+j(j+1)}{2f(f+1)} \end{equation}

  • getLandegfExact(l, j, f, s=0.5)

    Returns hyperfine-structure Landé g-factor \(g_F\) using precise \(g_L,\) \(g_S\), \(g_I\) values

\begin{equation} g_F=g_J\frac{f(f+1)-I(I+1)+j(j+1)}{2f(f+1)}+g_I\frac{f(f+1)+I(I+1)-j(j+1)}{2f(f+1)} \end{equation}

  • Ez = breitRabi(l, j, B, Ahfs, Bhfs)

Return exact Zeman energy levels for the l,j manifold via diagonalisation of \(\mathcal{H}=\mathcal{H}_\mathrm{hfs}+\mathcal{H}_\mathrm{Z}\) where \(\mathcal{H}_\mathrm{hfs}\) is the Hamiltonian for the hyperfine interaction given by \begin{equation} \mathcal{H}_\mathrm{hfs}=A_\mathrm{hfs}I\cdot J + B_\mathrm{hfs}\frac{3(I\cdot J)^2+3/2 I\cdot J -I^2J^2}{2I(2I+1)2J(2J+1)}, \end{equation} and \(\mathcal{H}_Z\) is the Hamiltonian for the Zeeman interaction given by \begin{equation} \mathcal{H}_Z=\frac{\mu_B}{\hbar}(g_J J_z+g_I I_z)B_z \end{equation}

For the \(j=1/2\) ground-state, this can be compared against the breitRabi formula [1]

:nbsphinx-math:`begin{equation}

E_{vert j=1/2,m_j,I,m_Irangle} = -frac{A_mathrm{hfs}(I+1/2)}{2(2I+1)}+g_Imu_B m B_z pm frac{A_mathrm{hfs}(I+1/2)}{2}left(1+frac{4mx}{2I+1}+x^2right)^{1/2},

end{equation}` where \(m=m_I\pm m_j\) and \begin{equation} x=\frac{g_j-g_I)\mu_B B_z}{A_\mathrm{hfs}(I+1/2)}. \end{equation}

[7]:
#Lande g-factors
gj = cs.getLandegj(0, 0.5)
gf3 = cs.getLandegf(0, 0.5, 4)
gf4 = cs.getLandegf(0, 0.5, 4)
print("Cs 6S1/2 gJ = %2.2f, F=3 gF = %2.2f, F=4 gF = %2.2f,"
      % (gj, gf3, gf4))
gj = cs.getLandegjExact(0, 0.5)
gf3 = cs.getLandegfExact(0, 0.5, 4)
gf4 = cs.getLandegfExact(0, 0.5, 4)
print("Cs 6S1/2 gJ = %f, F=3 gF = %f, F=4 gF = %f,"
      % (gj, gf3, gf4))

#Zeeman Interaction 6S1/2
#========================
#Define State of Interest
lg = 0; jg = 1 / 2;

#B-Field (G)
B = np.linspace(0, 5000, 101) * 1e-4
#Exact-Diagonalisation (Breit-Rabi)
en, f, mf = cs.breitRabi(6, lg, jg, B)

f, ax = plt.subplots(3, 1, figsize=(12, 10))
ax[0].plot(B * 1e4, en * 1e-9, 'k')
ax[0].set_xlabel(r'$B_z$ (G)')
ax[0].set_ylabel(r'$E$ (GHz)')
ax[0].set_title(r'Breit-Rabi Diagram $6S_{1/2}$')

#Plot Breit-Rabi Formula
#=======================
#   Steck Eq(25)
#                         dEhfs                dEhfs      4mx
#   E(j=1/2,mj;I,mI) = - ------- + gI*uB*m*B ± ----- (1+ ----- + x^2)^{1/2}
#                        2(2I+1)                 2        2I+1
#
#   where m=mI±mj to match sign above and x=(gj-gI)uB*B/dEhfs
uB = physical_constants["Bohr magneton in Hz/T"][0]
Ahfs, Bhfs = cs.getHFSCoefficients(6,0,0.5)
qhfs = Ahfs * (cs.I + jg)
gj = cs.getLandegjExact(lg, jg)
#Rescale to x
x=(gj - cs.gI) * uB * B / qhfs;

Ebr = np.zeros(np.shape(en))

ctr=-1
for mI in np.arange(-cs.I, 1+cs.I):
    for mj in np.arange(-jg, 1+jg):
        m=mI+mj
        ctr+=1
        if(abs(m)==(cs.I + jg)):
            Ebr[:,ctr] = qhfs * cs.I/(2 * cs.I + 1) + np.sign(mj)*1/2*(gj+2*cs.gI)*uB*B
        else:
            Ebr[:,ctr] = -qhfs/(2 * (2 * cs.I + 1)) + cs.gI * uB * m * B \
                         + np.sign(mj) * qhfs / 2 * np.sqrt(1 + 4 * m * x / (2 * cs.I + 1) + x**2)

ax[0].plot(B * 1e4, Ebr * 1e-9, 'r.');

ax[0].text(0,min(en[0,:])*1.5e-9,r"$F=3$",horizontalalignment='center',fontsize=12)
ax[0].text(0,max(en[0,:])*1.2e-9,r"$F=4$",horizontalalignment='center',fontsize=12)


#6P1/2

#Exact-Diagonalisation (Breit-Rabi)
B=np.linspace(0,2000,101)*1e-4
en,f,mf=cs.breitRabi(6,1,0.5,B)
ax[1].plot(B*1e4,en*1e-9,'k')
ax[1].set_xlabel(r'$B_z$ (G)')
ax[1].set_ylabel(r'$E$ (GHz)')
ax[1].set_title(r'Breit-Rabi Diagram $6P_{1/2}$')
ax[1].text(0, min(en[0, :]) * 1.5e-9,
           r"$F=3$",
           horizontalalignment='center', fontsize=12)
ax[1].text(0, max(en[0, :]) * 1.2e-9,
           r"$F=4$",
           horizontalalignment='center',
           fontsize=12)

#Exact-Diagonalisation (Breit-Rabi)
B = np.linspace(0, 500, 101) * 1e-4
en, f, mf = cs.breitRabi(6, 1, 1.5, B)
ax[2].plot(B * 1e4, en * 1e-9, 'k')
ax[2].set_xlabel(r'$B_z$ (G)')
ax[2].set_ylabel(r'$E$ (GHz)')
ax[2].set_title(r'Breit-Rabi Diagram $6P_{3/2}$')
ax[2].text(0, min(en[0,:]) * 2e-9,
           r"$F=2$",
           horizontalalignment='center', fontsize=12)
ax[2].text(0, max(en[0,:]) * 1.7e-9,
           r"$F=5$",
           horizontalalignment='center', fontsize=12)

plt.show()
Cs 6S1/2 gJ = 2.00, F=3 gF = 0.25, F=4 gF = 0.25,
Cs 6S1/2 gJ = 2.002319, F=3 gF = 0.249941, F=4 gF = 0.249941,
_images/ARC_3_1_additions_13_1.png
[ ]: