Doppler Averaging¶

This document discusses the methods rydiqule uses to implement doppler-averaging of modelling results. The theoretical background is provided for completeness, but is fairly standard. Rydiqule’s implementation of this is less obvious in order to optimize computational efficiency by fully leveraging numpy’s vectorized operations. The primary goal of this document is to clearly outline, in a single place, the underlying conventions used by rydiqule when doppler averaging.

Assumptions¶

Rydiqule currently makes implicit assumptions when modelling atomic systems that influence the treatment of doppler averaging. First, rydiqule’s equations of motion are single atom, meaning that atom-atom interactions are ignored. Second, rydiqule only solves these equations under the optically-thin assumption, meaning that input parameters do not change. In particular, absorption of the optical fields through an extended ensemble is not considered.

Both assumptions greatly simplify doppler averaging. Specifically, these assumptions allow us to further assume that atoms in different velocity classes do not interact, meaning that doppler averaging entails a simple weighted average of the atomic response at different velocities by the velocity distribution.

The basic process under these assumptions is as follows:

Calculate the density matrix solutions to the equations of motion for a wide range of velocity classes.
Weight the density matrix solutions with the Maxwell-Distribution, assigning the relative atomic density per velocity class to each solution.
Sum the weighted solutions.

Note that rydiqule assumes a three dimensional distribution of velocities, as is the case for a vapor.

Choosing Velocity Classes to Calculate¶

The primary influence of doppler velocities on the atomic physics is via Doppler shifts of the applied optical fields: defined as \(\Delta = \vec{k}\cdot\vec{v}\), where \(\vec{k}=2\pi/\lambda\cdot\hat{k}\) is the k-vector of the the optical field and \(\vec{v}\) is the velocity vector of the atomic velocity class in question. The Doppler shift experienced by an atom is due to the sum of Doppler shifts from each component projection of the atomic velocity relative to the k-vector of the field. In practice, the optical fields are often configured such that there is a basis where some of the k-vector components are zero, reducing the number of dimensions that need to be considered. In the simplest case, all optical fields are colinear, meaning all Doppler shifts are due to velocities projected along a single axis.

Choosing which velocity classes to calculate when performing doppler averaging is a fairly complex meshing problem. Our ultimate goal is to numerically approximate a continuous integral in one to three dimensions using a discrete sum approximation. For thermal vapors, where the spread of doppler velocities is large, resulting in Doppler shifts much greater than other detunings, linewidths, or Rabi frequenices in the problem, most velocity classes only contribute minor incoherent absorption to the final result. This allows for much coarser meshes. The difficulty lies in determining which velocity classes, a-priori, can participate in generally narrow coherent processes. This difficulty scales with the number of optical fields in the problem, as each new field increases the number of possible coherent resonances.

In rydiqule, we have striven to keep the specification of Doppler classes fairly flexible, with the constraint that velocity classes along each averaged dimension are the same (ie a rectangular grid). These classes are calculated in the doppler_classes() function, which contains options for complete user specification of all classes, as well as some convenience distributions that can be specified via a few parameters.

Maxwell-Boltzmann Distribution¶

Given that rydiqule gives single-atom solutions, the total atomic response for a given set of parameters is directly proportional to the total number of atoms represented by those parameters. The Maxwell-Boltzmann distribution is used to determine the fraction of the total population that is in each velocity class.

Following the conventions used by the Wikipedia page for the Maxwell-Boltzmann Distribution, the distribution of velocities for an ensemble of three dimensions is

\[f_\vec{v}(v_x, v_y, v_z) = \left(\frac{m}{2\pi k T}\right) e^{-\frac{m}{2kT}\left(v_x^2+v_y^2+v_z^2\right)}\]

Here, \(v_i\) represents the atomic velocity component along the cartesian axes, \(k\) is Boltzmann’s constant, \(T\) is the ensemble temperature in Kelvin, and \(m\) is the atomic mass in kilograms.

This distribution has a number of general properties. To begin, it is normalized such that integrating over all velocities in 3D space will give unity. There are also a few characteristic speeds associated with this 3D distribution: the most probable speed \(v_p\), the mean speed \(\langle v\rangle\), and the rms speed \(v_{rms}\).

\[\begin{split}\begin{align} v_p &= \sqrt{\frac{2kT}{m}}\\ \langle v\rangle &= \sqrt{\frac{8kT}{\pi m}}\\ v_{rms} &= \sqrt{\frac{3kT}{m}} \end{align}\end{split}\]

Note that speed is defined as the magnitude of the velocity vector. Also note that all of these quantities are related to each other by simple numerical prefactors. Finally, observe that the distribution above can be easily separated into each cartesian component.

\[f_\vec{v} = \prod_i \frac{1}{v_p \sqrt{\pi}} e^{-\frac{-v_i^2}{v_p^2}}\]

Note that the velocity distribution of each spatial component of the velocity is independently normalized. This allows us to readily produce the appropriate weighting distribution for 1, 2 and 3 dimensional averages as needed without having to perform redundant calculations of distributions where the density matrices to be averaged do not depend on a particular \(v_i\). This function is implemented in gaussian3d().

Averaging Velocity Classes¶

Given the above assumptions, the doppler average of the density matrix solutions is given by the integral

\[\bar{\rho_{ij}} = \int \rho_{ij}(\vec{v}) f_\vec{v} d^3v\]

This integral is numerically approximated via a finite sum.

\[\bar{\rho_{ij}} \approx \sum_{klm} \rho_{ij}(v_k, v_l, v_m) f_\vec{v}(v_k, v_l, v_m) \Delta v_k \Delta v_l \Delta v_m\]

In rydiqule, the weighting function \(f_\vec{v}\) is implemented in gaussian3d(), the volume element \(\Delta v_k \Delta v_l \Delta v_m\) is calculated as the product of the gradients along each axis as calculated by numpy.gradient() on the specified velocity classes.

We again note that when all k-vectors along a particular axis are zero, \(\rho_{ij}(v_k, v_l, v_m)\) is constant along that axis and that axis of the sum can be separated and assumed to sum to unity due to normalization of the weighting distribution along each dimension.

Rydiqule’s Implementation¶

Rydiqule’s implementation of Doppler averaging is optimized to minimize duplicate calculations and fully leverage numpy’s vectorized and broadcasting operations. The general steps are are follows:

Choose the doppler velocities to use for the mesh in the average.
Generate the Equations of Motion (EOMs) for the base zero velocity class using the machinery described in Equations of Motion Generation.
Generate the part of the EOMs that are proportional to the atomic velocity components \(v_i\). This is done by generating EOMs for the system with all parameters set to zero except for the optical detunings with associated non-zero k-vector components \(k_i\), multiplied by \(v_P\) to give the most probable Doppler shifts.
Generate the complete set of EOMs for all velocity classes via a broadcasting sum of the base EOMs with the Doppler EOMs multiplied by the normalized velocity classes along each axis. Each non-zero spatial axis that is to be summed over is pre-pended as an axis to the EOM tensor, as described in Stacking Conventions.
Solve the entire stack of EOMs.
Weight the EOMs according to their velocity classes via the Maxwell-Boltzmann distribution and the discrete velocity volume element, as described above.
Sum the solutions along the velocity axes.

Internally, rydiqule defines the necessary components for Doppler averaging via three quantities:

the normalized velocity classes \(d\), provided by doppler_classes()
the most probable speed \(v_P\) (in m/s), provided by the user as a class attribute
the optical k-vector \(\vec{k} = 2\pi /\lambda\cdot\hat{k}\) in (Mrad/m), provided for each coupling that has Doppler shifts to be averaged over

This construction has the benefit of allowing for meshes (ie velocity classes) to be defined in a general way relative to the distribution width \(v_P\), making them easily re-usable for any velocity distribution that obeys the Maxwell-Boltzmann distribution.

Migrating Doppler averaging from v1 to v2¶

With the release of v2 of rydiqule, how the user provides the above quantities has changed for both Sensor and Cell.

In v1, the 'kvec' parameter of the coupling was defined as the most probable Doppler shift vector (ie \(\vec{k}*v_P\)). This has been changed in v2 such that 'kvec' is now defined as the optical k-vector only (in units of Mrad/m), and \(v_P\) is provided separately at Sensor instantiation or by manually updating the vP. Put simply, moving Sensor simulations from v1 to v2 means no longer multiplying the k-vector by \(v_P\), and providing the vP attribute.

For Cell, v1 code followed the same old convention. Now that Cell has improved ARC integration, couplings in Cell take the 'kunit' argument which defines the unit propagation axis only. The \(v_P\) and \(2\pi/\lambda\) factors are calculated automatically and applied to any coupling with 'kunit' defined.