Simple Python tools for modelling few-level atom-light interactions

We will first review the semi-classical approach to near-resonant atom-light interactions using the density matrix to describe an atom-light system. We will derive the optical Bloch equations (OBEs) and look at how measurable quantities can be extracted to compare to experiments. In the next section we will look at two possible solution methods and the regimes in which they are valid. These methods are presented alongside relevant code snippets to give some insight into how the equations are translated into Python functions. Given lists of the system parameters, we will use Sympy to derive the OBEs from an interaction Hamiltonian and decay/dephasing operator, then expresses these as matrix equations and uses a SVD to solve them. We then look at some examples of using the code to generate plots of quantities of common interest in experiment, and use these to understand some of the phenomena observed in atom-light systems. We finish by demonstrating how the functions could be used to account for effects such as Doppler broadening in systems with a finite temperature.

In this tutorial, selected code snippets are presented to illustrate how the equations outlined are coded and solved. All functions outlined here are available in a module on GitHub (https://github.com/LucyDownes/OBE_Python_Tools) alongside Jupyter notebooks showing how the functions are created and some examples of their use. All codes used to generate the plots in this work are also available there. For most plots, parameters are in units of $2\pi\,\mathrm{MHz}$. The factor of 2π is omitted from plot axis labels for clarity. In examples, parameters were chosen such that a particular effect is clearly visible in the plot. Hence they are purely illustrative and not intended to represent any particular physical experiment. All of the code presented here was written in Python 3.8 but could very easily be made backwards compatible with Python 2.7. The code is in no way optimized for speed, it is intended to be as transparent as possible to aid insight. Variable names have been chosen to align with the notation used in this paper. Parameters (and hence variables) referring to a field coupling two levels i and j with have the subscript ij, whereas parameters that are relevant to a single atomic level n will have the single subscript n. Variable names with upper/lower-case letters denote upper/lower-case Greek symbols respectively, for example Gamma denotes Γ while gamma denotes γ.

In order to model how a single atom interacts with electromagnetic (EM) fields, we take a semi-classical approach in which the radiation field is considered as classical but the atom is considered quantum mechanically [15, 16].

The time evolution of the atomic wavefunction $|\psi\rangle$ is described by the time-dependent Schrödinger equation

$$\begin{equation} i\hbar\frac{\partial}{\partial t}|\psi\rangle = H|\psi\rangle. \end{equation} \tag{ 1 }$$

The total Hamiltonian H of the system is made up of two parts

$$\begin{equation} H = H_0 + H_{I}(t), \end{equation} \tag{ 2 }$$

where H₀ is the bare atomic Hamiltonian in the absence of any radiation and H_I describes the interaction with a time-dependent EM field. To begin with we will ignore the effects of dephasing processes and focus only on the coherent dynamics of the system.

First let us consider a two level atom with non-degenerate levels $|1\rangle$ and $|2\rangle$ having energies of $E_1 = \hbar\omega_1$ and $E_2 = \hbar\omega_2$, as shown on the left hand side of figure 1. These are known as the bare states of the system, and form a complete basis set. Any state of the atom $|\psi\rangle$ can then be written as a superposition of these bare states, which for our two-level atom can be written

$$\begin{equation} |\psi\rangle = c_1|1\rangle + c_2|2\rangle \end{equation} \tag{ 3 }$$

where the coefficients c_i represent the probability amplitude of being in state $|i\rangle$.

Zoom In Zoom Out Reset image size

In the basis of bare states we can write the atomic Hamiltonian as

$$\begin{equation} H_0 = \begin{pmatrix} \hbar\omega_1 & 0 \\ 0 & \hbar\omega_2 \end{pmatrix}. \end{equation} \tag{ 4 }$$

The eigenvalues of this Hamiltonian give the energies of the two levels. By defining the transition frequency $\omega_0 = \omega_2 - \omega_1$, we can shift the zero point energy and rewrite the Hamiltonian as

$$\begin{equation} H_0 = \begin{pmatrix} 0 & 0 \\ 0 & \hbar\omega_0 \end{pmatrix}. \end{equation} \tag{ 5 }$$

The energy difference between the states remains the same ($\hbar\omega_0$), we have just defined the energy of the lowest state to be zero.

The levels are coupled by a laser field with frequency ω. Since we are considering near-resonant excitation, we have that $\omega \approx\omega_0$, but we can quantify the difference between the laser frequency and the resonant frequency of the transition as $\Delta = \omega - \omega_0$. This is known as the detuning.

The interaction Hamiltonian H_I takes the form $-\mathbf{d}\cdot \mathbf{E}$ where E is the electric field and $\mathbf{d} = -e\mathbf{D}$ is the electric dipole operator.

Treating the EM field classically means it can be modelled as a simple oscillating field, meaning it can be written as $\mathbf{E} = E_0 \hat{\boldsymbol{e}}\cos(\omega t)$ where E₀ is the amplitude of the field, $\hat{\boldsymbol{e}}$ is a polarization vector and ω is the frequency of the field. Here we have neglected the spatial dependence of the field, an approximation known as the dipole approximation. This approximation is valid since the scale of the atomic wavefunction is much smaller than the wavelength of the applied fields, so the field can be considered uniform across the extent of the atom.

The field interacts with the atom via the dipole operator. This can be written as

$$\begin{equation} \hat{\mathbf{D}} = \sum_{i,j} \mathbf{D}_{ij}|i\rangle\langle\, j| \end{equation} \tag{ 6 }$$

where $|i\rangle$ and $|\,j\rangle$ are basis states of the bare atomic Hamiltonian. The term $\mathbf{D}_{ij} = \langle i|D|\,j\rangle$ is known as the matrix element [15]. Now we can write the interaction Hamiltonian as

$$\begin{equation} H_I = eE_0 \cos(\omega t)\sum_{i,j} \mathbf{D}_{ij}\cdot \hat{\mathbf{e}}|i\rangle\langle\, j|. \end{equation} \tag{ 7 }$$

It is convenient at this point to define a new quantity

$$\begin{equation} \Omega_{ij} = \frac{eE_0}{\hbar}\langle i|\mathbf{D}\cdot \hat{\mathbf{e}}|\,j\rangle \end{equation} \tag{ 8 }$$

which is known as the Rabi frequency. This defines the strength of the coupling between the levels i and j due to the field, and has units of frequency. It depends on both the intensity E₀ of the field and the polarization through $\hat{\mathbf{e}}$. From the definition of the Rabi frequency we can see that $\Omega_{ii} = 0$, and $\Omega_{ji} = \Omega_{ij}^*$.

Our interaction Hamiltonian now becomes

$$\begin{equation} H_I = \frac{\hbar\Omega}{2}\left(e^{i\omega t} + e^{-i\omega t}\right)\left[|1\rangle\langle 2| + |2\rangle\langle 1|\right], \end{equation} \tag{ 9 }$$

where we have dropped the subscript on the Rabi frequency since there is only one driving field. The terms in the round brackets describe the field, and can be thought of as describing the emission and absorption of a photon respectively. The operators in the square brackets describe the atom, and are related to a transition $|2\rangle\rightarrow |1\rangle$ and $|1\rangle\rightarrow |2\rangle$ (lowering and raising operators). The process whereby the atom makes the transition $|1\rangle\rightarrow |2\rangle$ after absorbing a photon, or makes the transition $|2\rangle\rightarrow |1\rangle$ by emitting a photon are more important than the other processes, so we will keep only the two terms that describe these [15]. This is an approximation known as the 'rotating wave approximation', and means the interaction Hamiltonian now has the form

$$\begin{equation} H_I = \frac{\hbar\Omega}{2}\left(e^{i\omega t}|1\rangle\langle 2| + e^{-i\omega t}|2\rangle\langle 1|\right). \end{equation} \tag{ 10 }$$

Combining this with the atomic Hamiltonian H₀, we get the total Hamiltonian

$$\begin{equation} H = \hbar \omega_0|2\rangle\langle 2| + \frac{\hbar\Omega}{2}\left(e^{i\omega t}|1\rangle\langle 2| + e^{-i\omega t}|2\rangle\langle 1|\right), \end{equation} \tag{ 11 }$$

which could also be expressed in matrix form as

$$\begin{equation} H = \frac{\hbar}{2}\begin{pmatrix} 0 & \Omega e^{i\omega t} \\ \Omega e^{-i\omega t} & 2\omega_0 \end{pmatrix}. \end{equation} \tag{ 12 }$$

In order to get rid of the explicit time dependence in the Hamiltonian, we can perform a frame transformation into a rotating frame. In this simple two-level case this can be done ¹ by introducing new variables $\tilde{c}_1(t)$ and $\tilde{c}_2(t)$ and writing the state vector as

$$\begin{equation} |\psi\rangle = \tilde{c}_1|1\rangle + \tilde{c}_2e^{-i\omega t}|2\rangle. \end{equation} \tag{ 13 }$$

By inserting this expression for the state vector and the total Hamiltonian into the time-dependent Schrödinger equation, it is possible to find equations for the time dependence of the coefficients $\tilde{c}_1(t)$ and $\tilde{c}_2(t)$. This can also be written in the form of a time-dependent Schrödinger equation with the new time-independent Hamiltonian

$$\begin{equation} \tilde{H} = \frac{\hbar}{2}\begin{pmatrix} 0 & \Omega \\ \Omega & -2\Delta \end{pmatrix}, \end{equation} \tag{ 14 }$$

where $\Delta = \omega - \omega_0$ is the detuning we defined earlier. This is the most commonly used form of the Hamiltonian, and the form we will use from here onwards.

This combined atom-light Hamiltonian allows us to explore properties of the coupled atom-light system. The effective energy levels of the atom-light system $E_{\pm}$, sometimes called eigenenergies, are found as the eigenvalues of this Hamiltonian and are given by

$$\begin{equation} E_{\pm} = -\frac{\hbar\Delta_{12}}{2}\pm\frac{\hbar}{2}\sqrt{{\Delta_{12}}^2+{\Omega}^2}. \end{equation} \tag{ 15 }$$

The eigenvectors of H are called 'eigenstates'. In the absence of the field (i.e. for $\Omega = 0$ and ω = 0) the eigenstates are the bare states of the system ($|1\rangle$ and $|2\rangle$) with eigenenergies of 0 and $\hbar\omega_0$ respectively, but when considering the coupled atom-field system the bare states are no longer eigenstates. For a near resonant field (i.e. for $\Delta\approx 0$) the eigenstates are given by

$$\begin{equation} |\pm\rangle = \frac{1}{\sqrt{2}}\left(|1\rangle \pm |2\rangle\right). \end{equation} \tag{ 16 }$$

These new states are known as dressed states, and are linear superpositions of the bare states. From the eigenenergies in equation (15) it can be seen that for increasing Ω the levels $E_{\pm}$ will be shifted further from each other; an effect known as the AC Stark effect or light shift [17]. This is illustrated on the right hand side of figure 1.

If we were to solve the Schrödinger equation for the coherent dynamics of the state vector $|\psi\rangle$ as defined earlier, we would find that the solutions are oscillatory. This oscillation of population between the two states are known as Rabi oscillations, and occur at a frequency defined by the Rabi frequency. This behavior can be seen in the upper panel of figure 2. However in real atomic systems there are incoherent processes such as decay between the levels which need to be taken into account which is not possible by considering the state vector alone. To this end we introduce the density operator defined as $\hat{\rho} = |\psi\rangle\langle\psi|$ and consider its evolution instead of that of the state vector.

Zoom In Zoom Out Reset image size

In the basis of bare states, the density operator can be written in matrix form. Using the definition of the state vector from equation (3) we have

$$\begin{equation} \hat{\rho} = \begin{pmatrix} |c_1|^2 & c_1 c_2^* \\ c_2 c_1^* & |c_2|^2 \end{pmatrix}. \end{equation} \tag{ 17 }$$

The diagonal elements of the density matrix describe the populations of the bare states $|1\rangle$ and $|2\rangle$, while the off-diagonal elements describe the coherence between states. It is worth noting that using the density matrix allows us to include more information than the state vector alone. For example it is possible to describe mixed states, states that are a statistical mixture of quantum states instead of a linear superposition, which cannot be described by a single state vector. For a pure state (one that can be described by a linear superposition of bare states) the trace of the square of the density matrix ($\mathrm{Tr}[\rho^2]$) is equal to 1, which is not the case for mixed states.

The coherent time evolution of the density matrix is given by the Liouville von–Neumann equation [15]

$$\begin{equation} i\hbar\frac{\partial\hat{\rho}}{\partial t} = \left[\hat{H}, \hat{\rho}\right], \end{equation} \tag{ 18 }$$

in which $[\hat{H}, \hat{\rho}]$ represents the commutator of the total Hamiltonian $\hat{H}$ and the density operator $\hat{\rho}$. This equation is analogous to the time-dependent Schrödinger equation for the state vector. In order to include dissipative effects we use the Master equation

$$\begin{equation} \frac{\partial\hat{\rho}}{\partial t} = -\frac{i}{\hbar}\left[\hat{H}, \hat{\rho}\right] + \hat{\mathcal{L}}(\hat{\rho}), \end{equation} \tag{ 19 }$$

where again $\hat{H}$ is the Hamiltonian of the atom-light system and the Lindblad operator $\hat{\mathcal{L}}$ describes the decay/dephasing in the system.

Zoom In Zoom Out Reset image size

The dissipative processes can be split into two separate parts; one part describing the spontaneous decay within the atom ($\hat{L}_{\mathrm{decay}}$) and another describing dephasing due to the finite linewidths of the driving fields ($\hat{L}_{\mathrm{dephasing}}$). The total operator will then be the sum of these parts

$$\begin{equation} \hat{\mathcal{L}} = \hat{L}_{\mathrm{decay}} + \hat{L}_{\mathrm{dephasing}}. \end{equation} \tag{ 20 }$$

First we focus on the operator $\hat{L}_{\mathrm{decay}}$ describing spontaneous decay between levels. The decay rate of population out of each level is given by $\Gamma_i$. For our two-level atom we assume that since the state $|1\rangle$ is the ground state, there is no decay out of this level. This amounts to saying $\Gamma_1 = 0$. We also assume that the system is closed; population remains within the two-level system and is not lost to an external reservoir. This means that any population that decays from $|2\rangle$ must end up back in $|1\rangle$. The amount of population transfer will depend on the population of the excited state. Hence the diagonal elements of the decay operator describing the change in the population of the ground and excited state will be given by $\Gamma_2\rho_{22}$ and $-\Gamma_2\rho_{22}$ respectively. Since decay can only ever cause a loss of coherence, the off-diagonal elements of the decay operator are always negative, and decay at the average rate of decay for the two levels. This means that for our two-level system the decay operator can be written in matrix form as

$$\begin{equation} \hat{L}_{\mathrm{decay}} = \begin{pmatrix} \Gamma_2\rho_{22} & -\frac{\Gamma_2}{2}\rho_{12} \\ -\frac{\Gamma_2}{2}\rho_{21} & - \Gamma_2\rho_{22} \end{pmatrix}. \end{equation} \tag{ 21 }$$

The action of decay is to damp the oscillations, meaning that eventually the system reaches an equilibrium, a so-called steady state. This is illustrated in the right-hand panel of figure 2. In continuous-wave experiments, this steady state is reached in very short timescales compared to the measurement time, so the steady state solution is often of most interest.

Secondly we establish the form of the dephasing operator $\hat{L}_{\mathrm{dephasing}}$ which accounts for the fact that the driving field can never be truly single-frequency, but will instead have a finite width. The dephasing due to the finite laser linewidths only affects the coherences (off-diagonal density matrix elements). Since we have only two levels, we only need to include effects due to a single field. In this case the dephasing operator takes the form

$$\begin{equation} \hat{L}_{\mathrm{dephasing}} = \begin{pmatrix} 0 & -\gamma_{12}\rho_{12} \\ -\gamma_{21}\rho_{21} & 0 \end{pmatrix} \end{equation} \tag{ 22 }$$

where $\gamma_{12} = \gamma_{21}$ is the linewidth of the field coupling the levels $|1\rangle$ and $|2\rangle$. Now we have expressions for the Hamiltonian and decay operator, we can substitute them into equation (19) to get an expression for the time evolution of the density matrix. In the steady-state, the solution is not time-dependent. This steady-state solution can be found by setting $\partial\hat{\rho}/\partial t = 0$ in equation (19) and solving for $\hat{\rho}$.

Once we have the values of the density matrix elements in the steady state, we would like to be able to extract useful measurable quantities to compare to experiment. The response of the medium to an applied field is characterized by the electric susceptibility χ. This macroscopic response of the system can be related to the individual atomic response as

$$\begin{equation} \chi = \frac{-2ND_{ji}^2}{\epsilon_0 \hbar\Omega_{ij}}\rho_{ji} \end{equation} \tag{ 23 }$$

where N is the number density (atoms per unit volume in the ensemble), $D_{ji} = \langle\, j|D|i\rangle$ is the dipole matrix element defined earlier, $\Omega_{ij}$ is the Rabi frequency of the field coupling the ith and jth levels and ρ_ji is the element of the density matrix describing the coherence between the ith and jth levels.

The real part of the complex susceptibility is proportional to the dispersion of the medium, while the imaginary part is proportional to the absorption coefficient α. Often the experimental quantity of interest is the amount of absorption of a laser after passing through the medium, most commonly the absorption of the probe laser which couples the lowest two levels. From the above relations we can see that the probe absorption $\alpha_{12}\propto-\Im[\rho_{21}]$, so often we will be interested in the steady state value of the density matrix element ρ₂₁.

When considering the simple case of a two-level system we did not have to pay attention to the configuration of the states, however as we add more states and fields coupling these states we need to think about how they are connected. For many-level systems there are multiple possible arrangements of energy levels, driving fields and decay pathways within an atom. For 3-level systems, three of the most commonly considered configurations are the Lambda (Λ), Vee (V) and ladder (Ξ) systems, shown in figure 4 from left to right. Here we will deal only with ladder systems, in which each level is only coupled to the levels directly above and below it in energy. This level configuration is common in experiments which excite atoms to high-lying Rydberg states through multi-photon excitation [18–21]. Obviously this simplified description means that hyperfine structure is not accounted for, but this is often not resolved in higher-lying states and so can be neglected.

Zoom In Zoom Out Reset image size

We can then extend our simple two-level model from the previous section to include an extra state $|3\rangle$ and driving field Ω₂₃ that couples the intermediate and highest energy states. A diagram of this simple three-level ladder scheme and the relevant parameters can be seen in figure 5. For this three-level system we apply the same methods and approximations as before and can write the total Hamiltonian in the interaction picture as

$$\begin{equation} H = \frac{\hbar}{2} \begin{pmatrix} 0 & \Omega_{12} & 0 \\ \Omega_{12} & -2\Delta_{12} & \Omega_{23} \\ 0 & \Omega_{23} & -2(\Delta_{12} + \Delta_{23}) \end{pmatrix}, \end{equation} \tag{ 24 }$$

where $\Omega_{12, 23}$ are the Rabi frequencies of the laser fields and $\Delta_{12, 23}$ are the corresponding detunings.

Zoom In Zoom Out Reset image size

In the ladder configuration we assume that each level decays to the level directly below it at a rate $\Gamma_n$, and that there is no decay out of the lowest energy level. In this case the diagonal elements of the decay operator $\hat{L}_{\mathrm{decay}}$ will be given by $L_{ii} = \Gamma_{i+1}\rho_{i+1, i+1} - \Gamma_{i}\rho_{ii}$ while the off-diagonal elements will be given by $L_{ij,\, i\neq j} = -\frac{\Gamma_i + \Gamma_j}{2}\rho_{ij}$ [15]. This means that for our example 3-level system we have

$$\begin{equation} \hat{L}_{\mathrm{decay}} = \begin{pmatrix} \Gamma_2 \rho_{22} & - \frac{\Gamma_2}{2}\rho_{12} & -\frac{\Gamma_3}{2}\rho_{13} \\[3pt] -\frac{\Gamma_2}{2}\rho_{21} & \Gamma_3 \rho_{33} - \Gamma_2\rho_{22} & -\frac{\Gamma_2 + \Gamma_3}{2}\rho_{23} \\[3pt] -\frac{\Gamma_3}{2}\rho_{31} & -\frac{\Gamma_3 + \Gamma_2}{2}\rho_{32} & - \Gamma_3\rho_{33} \end{pmatrix}, \end{equation} \tag{ 25 }$$

where we have already made the assumption that there is no decay out of the lowest level ($\Gamma_1 = 0$). The elements of the dephasing matrix $\hat{L}_{\mathrm{dephasing}}$ can be expressed as

$$\begin{equation} L_{ij, i\neq j} = -\left( \gamma_{i,i+1} + \dots + \gamma_{j-1, j}\right) \rho_{ij}, \end{equation} \tag{ 26 }$$

where $\gamma_{n,n+1}$ is the linewidth of the field coupling the nth and $(n+1)$th levels [22].

In Rydberg experiments it is often necessary to use three or more fields to reach the desired Rydberg state [18–21]. If another field is introduced, for example to couple to adjacent Rydberg states when performing electrometry experiments [23–25], then this again adds to the complexity of the system. These multi-level systems can often be thought of as ladder systems. Note that the codes could all be adapted to the other system configurations shown in figure 4, but the form of the Hamiltonian and decay/dephasing operators would need to be redefined.

Any linear system of equations can be expressed as a matrix of coefficients multiplied by a vector of variables. Our aim is to find a way to write the OBEs in this way, such that we can find the steady state solution using linear algebra.

3.1.1. Setting up the matrices.

First we set up a function to create the density matrix $\hat{\rho}$ itself. Given a number of levels, this function returns a matrix populated with Sympy symbolic objects which will allow us to set up the equations and manipulate the terms. For an n-level system the density matrix will have n² components.

def create_rho_matrix(levels = 3):

 rho_matrix = numpy.empty((levels, levels),

 dtype = 'object')

 for i in range(levels):

   for j in range(levels):

    globals()[

       'rho_'+str(i+1)+str(j+1)] = \

     sympy.Symbol(

      'rho_'+str(i+1)+str(j+1))

    rho_matrix[i,j] = globals()[

       'rho_'+ str(i+1)+str(j+1)]

 return numpy.matrix(rho_matrix)

Next we set up the Hamiltonian of the interaction between the atom and the light, in the form shown in equations (14) and (24).

def Hamiltonian(Omegas, Deltas):

 levels = len(Omegas)+1

 H = numpy.zeros((levels,levels))

 for i in range(levels):

    for j in range(levels):

      if numpy.logical_and(i==j, i!=0):

        H[i,j] = -2*(numpy.sum(

          Deltas[:i]))

      elif numpy.abs(i-j) == 1:

        H[i,j] = Omegas[numpy.min([i,j])]

 return numpy.matrix(H/2)

This function returns a matrix of floats. Note that we have omitted the factor of $\hbar$ when we create this Hamiltonian. This is because it will cancel with the factor of $\hbar$ in the Master equation later on. We could include it explicitly, but because it is often many orders of magnitude smaller than the values of Ω and Δ it makes the overall values within the Hamiltonian very small and we risk errors due to floating point precision.

Next we set up the matrix form of the decay operator $\hat{\mathcal{L}}$. As we did earlier, we will split this into two separate parts for simplicity. First we focus on $\hat{L}_{\mathrm{decay}}$ describing spontaneous decay. The function takes a list of values for the decay rates of the excited states $\Gamma_i$. For an n-level system, this list will have n − 1 entries as we always assume no decay from the lowest energy level ($\Gamma_1 = 0$). Again the function returns an n × n matrix populated by multiples of Sympy symbolic objects.

def L_decay(Gammas):

 levels = len(Gammas)+1

 rhos = create_rho_matrix(levels = levels)

 Gammas_all = [0] + Gammas

 decay_matrix = numpy.zeros((levels, levels),\

 dtype = 'object')

 for i in range(levels):

   for j in range(levels):

     if i != j:

      decay_matrix[i,j] = -0.5*(

        Gammas_all[i]+Gammas_all[j])\

        *rhos[i,j]

    elif i != levels - 1:

      into = Gammas_all[i+1]*\

      rhos[1+i,j+1]

      outof = Gammas_all[i]*rhos[i, j]

      decay_matrix[i, j] = into - outof

     else:

      outof = Gammas_all[i]*rhos[i, j]

      decay_matrix[i,j] = - outof

 return numpy.matrix(decay_matrix)

We can create the matrix $\hat{L}_{\mathrm{dephasing}}$ in a similar way, again from a list of values of γ_ij .

def L_dephasing(gammas):

 levels = len(gammas)+1

 rhos = create_rho_matrix(levels = levels)

 deph_matrix = numpy.zeros((levels, levels),\

 dtype = 'object')

 for i in range(levels):

   for j in range(levels):

     if i != j:

        deph_matrix[i,j] = -(numpy.sum(

          gammas[numpy.min([i,j]):\

          numpy.max([i,j])]))*rhos[i,j]

 return numpy.matrix(deph_matrix)

Now we have expressions for the Hamiltonian and Lindblad operators, we can put them into the Master equation to get an expression for the time evolution of the density matrix. As we saw in equation (19) we have that

$$\begin{equation} \frac{\partial\hat{\rho}}{\partial t} = -\frac{i}{\hbar}(\hat{H}\hat{\rho} - \hat{\rho}\hat{H}) + \hat{L}_{\mathrm{decay}} + \hat{L}_{\mathrm{dephasing}}. \end{equation} \tag{ 27 }$$

[Note that the $\hbar$ omitted earlier from the Hamiltonian would have cancelled with the $\hbar$ here, so when we define the function we omit the factor of $1/\hbar$.]

def Master_eqn(H_tot, L):

 levels = H_tot.shape[0]

 dens_mat = create_rho_matrix(levels = levels)

 return -1j*(

   H_tot*dens_mat - dens_mat*H_tot) + L

This gives us an expression for the time evolution of each of the components of the density matrix in terms of other components. To make this complex system of equations more convenient to solve we want to write them in the form

$$\begin{equation} \frac{\partial\hat{\rho}_{\mathrm{vect}}}{\partial t} = \hat{M}\hat{\rho}_{\mathrm{vect}} \end{equation} \tag{ 28 }$$

where $\hat{\rho}_{\mathrm{vect}}$ is a column vector containing all of the elements of $\hat{\rho}$, and $\hat{M}$ is a matrix of coefficients. First we set up a function to create $\hat{\rho}_{\mathrm{vect}}$, a list of the components of the density matrix $\rho_{i,j}$.

def create_rho_list(levels = 3):

 rho_list = []

 for i in range(levels):

    for j in range(levels):

      globals()[

        'rho_'+str(i+1)+str(j+1)] = \

      sympy.Symbol(\

        'rho_'+str(i+1)+str(j+1))

      rho_list.append(globals()[

        'rho_'+str(i+1)+str(j+1)])

 return rho_list

We can then create the matrix of coefficients $\hat{M}$ by using the Sympy functions expand and coeff. The expand function collects together terms containing the same symbolic object. The coeff function allows us to extract the multiplying coefficient for a specified symbolic object. This works since the only terms we have are linear with respect to the density matrix elements, e.g. we never have terms in $\rho_{ij}^2$ or higher.

def OBE_matrix(Master_matrix):

 levels = Master_matrix.shape[0]

 rho_vector = create_rho_list(levels = levels)

 coeff_matrix = numpy.zeros((levels**2, \

 levels**2), dtype = 'complex')

 count = 0

 for i in range(levels):

   for j in range(levels):

     entry = Master_matrix[i,j]

     expanded = Sympy.expand(entry)

     for n,r in enumerate(rho_vector):

       coeff_matrix[count, n] = complex(

         expanded.coeff(r))

     count += 1

 return coeff_matrix

3.1.2. Time-dependent solutions.

We now have an expression for $\partial{\hat{\rho}}_{\mathrm{vect}}/\partial t$ in terms of the matrix $\hat{M}$ multiplied by $\hat{\rho}_{\mathrm{vect}}$. The value of $\hat{\rho}_{\mathrm{vect}}$ at any time t can be found from

$$\begin{equation} \hat{\rho}_{\mathrm{vect}}(t) = \exp\left(\hat{M}t\right)\hat{\rho}_{\mathrm{vect}}(0) \end{equation} \tag{ 29 }$$

where $\hat{\rho}_{\mathrm{vect}}(0)$ is the value of $\hat{\rho}_{\mathrm{vect}}$ at time t = 0. Since we already have a way of creating the matrix $\hat{M}$, it is trivial to define functions to find the time-dependent behavior of $\hat{\rho}_{\mathrm{vect}}$.

import scipy.linalg as linalg

def time_evolve(operator, t, p_0):

 exponent = operator*t

  #scipy.linalg.expm does matrix exponentiation

 term = linalg.expm(exponent)

 return numpy.matmul(term, p_0)

3.1.3. SVD.

Often in experiments the measurement timescales are long such that transient oscillatory behavior such as Rabi oscillations are not observed, and the system reaches a so-called steady-state. In the steady state, $\partial{\hat{\rho}}_{\mathrm{vect}}/\partial t = 0$ so we need to solve the expression

$$\begin{equation} \hat{M}\hat{\rho}_{\mathrm{vect}} = 0. \end{equation} \tag{ 30 }$$

This can be solved numerically in many different ways ² , but here we will solve the system by performing a singular value decomposition (SVD) on the matrix $\hat{M}$ [26]. The full details of the theory of the SVD are beyond this discussion, it is suffice to say that the SVD transforms $\hat{M}$ into three matrices such that

$$\begin{equation} \hat{M} = U\Sigma V^{\dagger} \end{equation} \tag{ 31 }$$

where $V\,^{\dagger}$ is the conjugate transpose of V. The matrix Σ is diagonal, the elements corresponding to the singular values of $\hat{M}$. The solution to the system of equations is then the column of V corresponding to the zero singular value. Since we have the same number of equations as we have unknowns we expect only one non-trivial solution. If none of the singular values are zero ($\hat{M}$ is non-singular) then there is no non-trivial solution. Floating point precision means that it is unlikely that any of the values will be identically zero, so we set a limit and assume that values smaller than this are due to floating point errors.

Zoom In Zoom Out Reset image size

We perform the SVD using the built-in numpy.linalg.svd routine which returns the matrices U and $V\,^{\dagger}$ along with an array of singular values Σ. Before returning the values of the elements of the density matrix we normalize the sum of the populations to be 1. The function returns an array of values for $\hat{\rho}_{\mathrm{vect}}$, the order of which is the same as the elements of the output of the create_rho_list function.

def SVD(coeff_matrix):

 levels = int(numpy.sqrt(

   coeff_matrix.shape[0]))

 u,sig,v = numpy.linalg.svd(coeff_matrix)

 abs_sig = numpy.abs(sig)

 minval = numpy.min(abs_sig)

 if minval>1e-12:

    # to account for floating point precision

   print('ERROR - Matrix is non-singular')

   return numpy.zeros((levels**2))

 index = abs_sig.argmin()

 rho = numpy.conjugate(v[index,:])

  # extract populations and normalise

 pops = numpy.zeros((levels))

 for l in range(levels):

   pops[l] = numpy.real(rho[(l*levels)+l])

 t = 1/(numpy.sum(pops))

 rho_norm = rho*t.item()

  #item() gets Python scalar from array

 return rho_norm

For ease of calculation we combine all of the previous functions into one function that takes lists of the different parameters (Ω, Δ, Γ and γ) and outputs the steady state solution for the density matrix in the order that rho_list is created. Note that it outputs the solution of the entire density matrix, so we need to extract the element of interest (usually ρ₂₁).

def steady_state_soln(

 Omegas, Deltas, Gammas, gammas = []):

 L_atom = L_decay(Gammas)

  # include dephasing if gammas values given

 if len(gammas) != 0:

   L_laser = L_dephasing(gammas)

   L_tot = L_atom + L_laser

 else:

   L_tot = L_atom

  # build Hamiltonian

 H = Hamiltonian(Omegas, Deltas)

  # insert H into Master equation

 Master = Master_eqn(H, L_tot)

  # create matrix of coefficients

 rho_coeffs = OBE_matrix(Master)

  # perform SVD to get steady state solution

 soln = SVD(rho_coeffs)

 return soln

All we need to do now is pass lists of the system parameters to this function.

Obviously because the matrix method requires solving an eigenvalue problem for each set of parameters it is not fast, however it does allow for the computation of results without any further assumptions about the system parameters. If the system is in the weak-probe regime, defined as the regime where $\Omega_{12}^2/\Gamma_3(\Gamma_2 + \Gamma_3)\ll 1$, a number of assumptions can be made which allow analytic expressions for the probe coherence to be found. Since these assumptions involve treating the probe field as a weak perturbation, we keep only terms linear in Ω₁₂. In this case it can be said that no population transfer occurs, so we have that

$$\begin{equation} \rho_{ii} = \begin{cases} 1,& i = 1 \\ 0,& i>1. \end{cases} \end{equation} \tag{ 32 }$$

For a three-level system, the probe coherence in the weak-probe regime can be expressed as

$$\begin{align} \rho_{21} & = \frac{-i\Omega_{12}}{2}\left[-i\Delta_{12} +\frac{\Gamma_2}{2}+\gamma_{12}+\frac{\Omega_{23}^2}{4}\left[-i(\Delta_{12} + \Delta_{23})\right.\right.\nonumber\\ & \left.\left.\qquad\qquad\qquad +\frac{\Gamma_3}{2} + \gamma_{12} + \gamma_{23}\right]^{-1}\right]^{-1}. \end{align} \tag{ 33 }$$

which can easily be encapsulated into a function.

3.2.1. Analytic solution for an arbitrary number of levels.

It is possible to formulate an expression for the analytic solution in the weak-probe regime for a ladder scheme with any number of levels. This has to be done iteratively.

For an system with n levels (hence n − 1 driving fields) we have an expression for the probe coherence which needs to be evaluated for all fields $m = 1, 2 \dots n-2, n-1$

$$\begin{equation} \rho_{21} = \frac{2iK_1}{\Omega_{12}} \end{equation} \tag{ 34 }$$

where

$$\begin{equation} K_m = \left(\frac{\Omega_{m, m+1}}{2}\right)^2\left[Z_m + K_{m+1}\right]^{-1} \end{equation} \tag{ 35 }$$

and

$$\begin{equation} Z_m = \sum_{k = 1}^m \left(i\Delta_{k,k+1} - \gamma_{k,k+1}\right) -\frac{\Gamma_{m+1}}{2}. \end{equation} \tag{ 36 }$$

For example, if we have a three-level scheme the expression will be

$$\begin{equation} \rho_{21} = \frac{2i}{\Omega_{12}}\left(\frac{\Omega_{12}}{2}\right)^2\left[Z_1 + \left(\frac{\Omega_{23}}{2}\right)^2\left[Z_2\right]^{-1}\right]^{-1} \end{equation} \tag{ 37 }$$

where

$$\begin{equation*} Z_1 = i\Delta_{12} - \gamma_{12} - \Gamma_2/2 \end{equation*}$$

and

$$\begin{equation*} Z_2 = i(\Delta_{12} + \Delta_{23}) - (\gamma_{12} + \gamma_{23}) - \Gamma_3/2, \end{equation*}$$

which is identical to equation (33).

This has been encapsulated into the function fast_n_level which takes the same sets of parameters as the steady_state_soln matrix method and returns the value of ρ₂₁.

def term_n(n, Deltas, Gammas, gammas = []):

  # Valid for n>0

 if len(gammas) == 0:

   gammas = numpy.zeros((len(Deltas)))

 return 1j*(numpy.sum(Deltas[:n+1])) - \

 (Gammas[n]/2) - numpy.sum(gammas[:n+1])

def fast_n_level(

 Omegas, Deltas, Gammas, gammas = []):

 n_terms = len(Omegas)

 term_counter = numpy.arange(n_terms)[::-1]

 func = 0

 for n in term_counter:

   prefactor = (Omegas[n]/2)**2

   term = term_n(n, Deltas, Gammas, gammas)

   func = prefactor/(term+func)

 return (2j/Omegas[0])*func

The time taken for both the matrix and iterative solutions to find the weak-probe coherence in an n-level system is shown in figure 6. The matrix solution is around three orders of magnitude slower than the iterative solution, and has a power law dependence on the number of levels. The iterative solution scales linearly with the number of levels, as can be seen in the inset of figure 6. Obviously the exact times will vary with the hardware used, but it gives an illustration of the usefulness of the iterative method as a tool to fast computation for an arbitrary number of levels. Note that none of the code presented here has been optimized for speed. The aim is to get an insight into how the equations are constructed and solved in the simplest way possible, so there are likely much faster but possibly more complex ways to solve the equations presented here.

Below is a simple example of using the matrix method (via the steady_state_soln function) to look at the probe absorption as a function of probe detuning (Δ₁₂) in a three-level system. We create an array of values for Δ₁₂, then for each of these values we calculate the steady state solution, extract the value of ρ₂₁ and store the imaginary part in an array.

Omegas = [0.1,4] # 2pi MHz

Deltas = [0,0] # 2pi MHz

Gammas = [1,0.1] # 2pi MHz

Delta_12 = numpy.linspace(-10,10 ,200) # 2pi MHz

# create empty array to store solution

probe_abs = numpy.zeros((len(Delta_12)))

for i, p in enumerate(Delta_12):

 Deltas[0] = p # update value of Delta_12

  # Pass parameters to the function to find

  # the steady state solution for the

  # density matrix

 solution = steady_state_soln(

   Omegas, Deltas, Gammas)

  # Store the imaginary part of rho_12 in

  # an array

 probe_abs[i] = numpy.imag(solution[1])

The results of running the code for different parameters are shown in figure 7. When there is only one field ($\Omega_{23} = 0\,\mathrm{MHz}$), the system can be thought of as a two-level system. In this case the probe is maximally absorbed on resonance ($\Delta_{12} = 0\,\mathrm{MHz}$) and we see a Lorentzian absorption feature as a function of probe detuning, the width of which is dependent on the ratio between Γ₂ and Ω₁₂.

Zoom In Zoom Out Reset image size

Setting $\Omega_{12} = 0$ in equation (24) and solving for the eigenvectors yields the dressed states

$$\begin{equation} |\pm\rangle = (\Delta_{23} \pm\sqrt{{\Delta_{23}}^2 + {\Omega_{23}}^2})|2\rangle + \Omega_{23} |3\rangle, \end{equation} \tag{ 38 }$$

with eigenenergies given by

$$\begin{equation} E_{\pm} = -\frac{\hbar}{2}\left(2\Delta_{12} + \Delta_{23} \pm \sqrt{\Delta_{23}^2 + \Omega_{23}^2}\right). \end{equation} \tag{ 39 }$$

When both beams are on resonance ($\Delta_{12,23} = 0$) the expression for the dressed states becomes $|\pm\rangle = 1/\sqrt{2}\left(|3\rangle \pm |2\rangle\right)$ and we see that their eigenenergies are separated by $\hbar\Omega_{23}$. This means the probe beam that was resonant with the bare state $|2\rangle$ is no longer resonant with either of the dressed states of the coupled atom-field system, so a transparency occurs. When the probe becomes resonant with one of the dressed states we see absorption, hence in figure 7 we see two absorption peaks separated by Ω₂₃. This is known as EIT or Autler-Townes (AT) splitting, depending on the parameters used. The distinction between Autler–Townes splitting and EIT is not well defined but it is generally accepted that EIT requires the condition of weak fields (i.e. small $\Omega_{12,23}$) [27]. As Ω₂₃ increases, Autler–Townes splitting becomes the dominant effect and the transparency window widens. For $\Delta_{23} = 0$ both dressed states contain an equal amount of the bare state $|2\rangle$ so the heights of the features in the absorption spectrum will be equal, as seen in figure 7. Any detuning of the second field ($\Delta_{23} \neq 0$) will result in an unequal mixture of the bare states in the dressed states as seen in equation (38) hence an asymmetric coupling to the ground state. This leads to features of unequal heights in the absorption spectra, as seen in figure 8. The negatively (positively) detuned peaks are due to coupling with the $|-\rangle$ ($|+\rangle$) dressed state. The peaks also become further apart as Δ₂₃ is increased, since the eigenenergies are separated by $\Delta_{23} + \sqrt{\Delta_{23}^2 + \Omega_{23}^2}$.

Zoom In Zoom Out Reset image size

Using the full matrix solution, it is possible to extract all of the parameters of the density matrix. For example if we wanted to look at the populations of the states we just need to extract the relevant terms from the output. For the three-level case, the populations of the states $|1\rangle,\,|2\rangle$ and $|3\rangle$ will be the first, fifth and ninth elements respectively. An example of looking at state population as a function of probe detuning and Rabi frequency is shown in figure 9. When the probe beam is weak ($\Omega_{12}^2\ll\Gamma_3(\Gamma_2 + \Gamma_3)$) we see very little population transfer; one of the assumptions made when imposing the 'weak-probe regime'. As Ω₁₂ is increased (fainter lines in figure 9), more population is transferred between states.

Zoom In Zoom Out Reset image size

In section 3.2 we said that if the weak-probe condition was satisfied ($\Omega_{12}^2\ll\Gamma_3(\Gamma_2 + \Gamma_3)$) it was possible to find analytic solutions for the probe coherence, but that these solutions were not valid outside of the weak-probe regime. We compare the results of the full matrix solution and analytic solution in figure 10. In the first case upper panel of figure 10), the weak-probe approximation is valid (since $\Omega_{12}^2/\Gamma_3(\Gamma_2 + \Gamma_3) = 1/21 \ll 1$) and so the two methods give equal results. However the analytic solution is many orders of magnitude faster. If the weak-probe condition is not satisfied, the difference in results becomes apparent. This is seen in the lower panel of figure 10. The weak-probe condition is no longer satisfied (since $\Omega_{12}^2/\Gamma_3(\Gamma_2 + \Gamma_3) = 2.3 \gt 1$) and so the results of the analytic solution no longer match the full matrix solution.

Zoom In Zoom Out Reset image size

Now we can look at a four-level scheme with three laser fields ($\Omega_{12}, \Omega_{23}, \Omega_{34}$). For the case where $\Omega_{34} = 0\,\mathrm{MHz}$ we recover the spectrum for the three-level case as there is no coupling to the $4^{\mathrm{th}}$ level, as shown by the light purple solid line in figure 11.

Zoom In Zoom Out Reset image size

As we increase the value of Ω₃₄ a central feature appears and grows while the other features reduce in size, as shown by the dashed and dotted lines.

Similarly, for four levels we can compare the results of the full matrix solution to the analytic one, as shown in figure 12. If we go far outside the weak probe case, then the need for the full matrix solution becomes clear.

Zoom In Zoom Out Reset image size

After calculating the response of atoms moving with different velocities, their relative abundances need to be accounted for. The atoms will have a distribution of velocities as given by the Maxwell-Boltzmann distribution

$$\begin{equation} f_{v}(v_i) = \sqrt\frac{m}{2\pi k_{B}T} \exp \left(-\frac{mv_{i}^2}{2k_{B}T}\right), \end{equation} \tag{ 40 }$$

where m is the mass of an atom, k_B is the Boltzmann constant and T is the temperature of the vapour in Kelvin. We can either define this function explicitly, or recognize that this is a Gaussian distribution centered on 0 with width $\sqrt{\frac{k_{B}T}{m}}$ and use a built-in function such as scipy.stats.norm.

To compute the Doppler averaged response of the vapour we need to take the response of each velocity class and multiply its contribution to the probe absorption by the probability of an atom having that velocity (given by the Maxwell–Boltzmann distribution we defined earlier). Then, for each value of Δ₁₂ we integrate (sum) over all the weighted velocity class contributions.

5.1.1. Choice of velocity range.

In figure 14 we show an example of calculating the absorption of a weak probe beam in a thermal vapour of three-level atoms. The relative directions of the beams are taken into account with the direction variable which can take values of ±1 (same sign = co-propagating beams, opposite sign = counter-propagating beams). The actual value of the sign does not matter (as the velocity distribution will be symmetric about 0), as long as we are consistent with the signs of the wavevectors relative to each other. We also need to specify the wavelengths of the two beams, here chosen to be $\lambda_{12} = 780\,\mathrm{nm}$ and $\lambda_{23} = 480\,\mathrm{nm}$ to match commonly used rubidium two-photon EIT experiments. The beams have been set to be counter-propagating. These calculations can be quite computationally intensive, especially when using the full matrix method. Since the purpose here is to illustrate the effects of atomic motion, we will work in the weak probe regime ($\Omega_{12}^2 \ll \Gamma_3(\Gamma_2 + \Gamma_3)$) and use the iterative solution method (via the fast_n_level function) to save time.

Zoom In Zoom Out Reset image size

Omegas = [0.1,10] # 2pi MHz

Deltas = [0,0] # 2pi MHz

gammas = [0.1,0.1] # 2pi MHz

Gammas = [10,1] # 2pi MHz

lamdas = numpy.asarray([780, 480])*1e-9 # m

# unitless direction vectors

directions = numpy.asarray([1,-1])

ks = directions/lamdas # 2pi m -1

Delta_12 = numpy.linspace(-50,50, 500) # 2pi MHz

velocities = numpy.linspace(-150, 150, 501) # m/s

probe_abs = numpy.zeros((len(velocities), \

len(Delta_12)))

for i, v in enumerate(velocities):

  # loop over all atomic velocities

 for j, p in enumerate(Delta_12):

   # loop over all values of probe detuning

   # (lab frame)

  Deltas[0] = p

   # create empty array to store modified

   # detunings

  Deltas_eff = numpy.zeros(len(Deltas))

  for n in range(len(Deltas)):

   Deltas_eff[n] = Deltas[n]+(

     ks[n]*v)*1e-6 # in MHz

   # pass parameters to function

   #(including modified detuning)

  solution = fast_n_level(

    Omegas, Deltas_eff, Gammas, gammas)

  probe_abs[i,j] = -numpy.imag(solution)

The two dotted lines in the upper plot of figure 14 are the positions at which the dressed states $|+\rangle$ and $|-\rangle$ as defined in equation (39) are resonant with the probe. This will occur when

$$\begin{equation} E_{\pm} = -\Delta_{12_{\mathrm{eff}}} - \frac{\Delta_{23_{\mathrm{eff}}}}{2} \pm \frac{1}{2}\sqrt{(\Delta_{23_{\mathrm{eff}}}^2 + \Omega_{23}^2)} = 0, \end{equation} \tag{ 41 }$$

where the effective detunings are given by

$$\begin{align*} \Delta_{12_{\mathrm{eff}}} = &\, \Delta_{12} + v s_{12}k_{12} \\ \Delta_{23_{\mathrm{eff}}} = &\, \Delta_{23} + v s_{23}k_{23}. \end{align*}$$

$k_{12,23} = 1/\lambda_{12,23}$ are the magnitudes of the wavevectors (ignoring the factor of 2π) and $s_{12,23} = \pm 1$ represent the relative direction of the beams. $E_{\pm}$ are the eigenvalues of the Hamiltonian with $\Omega_{12} = 0$ and detunings replaced by the effective detunings. In the limit of large detuning where $\Delta_{23_{\mathrm{eff}}}\gg \Omega_{23}$, the equations for the lines are

$$\begin{align} v_+ = & \frac{-\Delta_{12}}{s_{12}k_{12}} \end{align} \tag{ 42 }$$

$$\begin{align} v_- = & \frac{-(\Delta_{12}+\Delta_{23})}{s_{12}k_{12} + s_{23}k_{23}}. \end{align} \tag{ 43 }$$

If fewer velocity classes were included, the integrated result would look very different and becomes meaningless. We can however use the knowledge of the equations of the asymptotic lines (found via the eigenvalues) to give an idea of the range of velocities that need to be sampled.

We can use the same method to look at the effects of atomic motion on four-level atoms. As with the previous case, we consider systems that are within the weak probe limit so we can use the faster iterative solution. Since we are only considering the 1D case we can again specify beam direction using the direction variable.

In the upper plot of figure 15 we can see three asymptotic 'branches', whose positions can again be found by finding the eigenvalues of the Hamiltonian in the case where $\Omega\ll\Delta_{\mathrm{eff}}$, and equating them to zero. Rewriting the expressions in terms of v we have

$$\begin{align} v_1 = & \frac{-\Delta_{12}}{s_{12} k_{12}} \end{align} \tag{ 44 }$$

$$\begin{align} v_2 = & \frac{-\Delta_{12}-\Delta_{23}}{s_{12} k_{12} + s_{23} k_{23}} \end{align} \tag{ 45 }$$

$$\begin{align} v_3 = & \frac{-\Delta_{12}-\Delta_{23}-\Delta_{34}}{s_{12} k_{12} + s_{23} k_{23} + s_{34} k_{34}}, \end{align} \tag{ 46 }$$

where again $k_{12,23,34} = 1/\lambda_{12,23,34}$ are the magnitudes of the wavevectors and $s_{12,23,34} = \pm 1$ represent the relative directions of the beams. Knowing the equations of these asymptotic lines can be useful when deciding how many velocity classes need to be taken into account. In the above example, there is not much extra information gained by including velocity classes $\gt\!\!100\,\mathrm{m\,s}^{-1}$ or $\lt\!\!-100\,\mathrm{m\,s}^{-1}$.

Zoom In Zoom Out Reset image size