A fully quantum calculation of broadening and shifting coefficients of the D1 and D2 spectral lines of alkali-metal atoms colliding with noble-gas atoms

We use the Baranger model to compute collisional broadening and shift rates for the D1 and D2 spectral lines of M + Ng, where M = K, Rb, Cs and Ng = He, Ne, Ar. Scattering matrix elements are calculated using the channel packet method, and non-adiabatic wavepacket dynamics are determined using the split-operator method together with a unitary transformation between adiabatic and diabatic representations. Scattering phase shift differences are weighted thermally and are integrated over temperatures ranging from 100 K to 800 K. We find that predicted broadening rates compare well with experiment, but shift rates are predicted poorly by this model because they are extremely sensitive to the near-asymptotic behavior of the potential energy surfaces.


Introduction
This research uses the Baranger model to simulate collisional line broadening of relevant alkali vapor-noble gas mixtures under varied conditions (e.g., varying temperature and pressure). The particular mixtures of interest are those in typical use in optically pumped alkali laser (OPAL) systems. The Baranger model builds directly from the work of Jabloński [1] and is a fully quantum-mechanical model. Like Anderson-Talman [2][3][4], Baranger assumes the reference frame of the emitter/absorber atom. Baranger uses the impact approximation, which assumes that the duration of a 1 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. collision is much shorter than the time between collisions. The impact limit forces one to focus more on the core features of the collision-broadened spectral line than on the wings or satellite features. The only predictions found in the literature have calculated broadening and shift under adiabatic potentials for lighter alkali (Li, Na, K) perturbed by He [5,6], primarily for astrophysical application. These calculations are limited by the semi-classical treatment of collisions [7,8] and the neglect of fine structure transitions [6]. Recent interest in the behavior of the non-adiabatic fine structure transitions of atomic alkali as they collide with noble gases has been generated by applications in astrophysics and the development of OPALs [9][10][11][12][13][14][15].
This research develops a model for line broadening in which the time evolution of the alkali vapor-noble gas system is handled through wavepacket propagation. The quantummechanical time-evolution operator for the system is governed by the Hamiltonian, and we will use a fast Fourier transform (FFT) and its inverse to transform the wavefunction of the system between the momentum and position representations, as appropriate, in order to operate with the momentumdependent and position-dependent portions, respectively, of the time-evolution operator. The normal method of examining atomic collisions is to approximate the colliding system of two atoms as a diatomic molecular system. This allows one to describe the system using appropriate Hund's states [16,17]. It is through this approximation to molecular dynamics that we will utilize difference potentials in the context of this research.
This research exhibits several new features which set it apart from the current state of the field. First, the full ab initio potential energy surfaces are used; these potential energy surfaces have been calculated through many-body calculations by Blank [18,19]. Second, collisions are treated quantummechanically and non-adiabatically and include spin-orbit and Coriolis coupling. Third, calculations are made with the only approximations beyond those of the impact limit being those imposed by the Boltzmann (thermal) distribution of energies [20]. Finally, this research uses as its foundation the same potential energy surfaces as the work of Blank [19,21], which allows us to compare results from the Anderson-Talman and Baranger models for the same set of potential energy surfaces. We can use this approach to deconvolve differences in the models from differences in computational methods which might yield dissonant results.

S-matrix elements
The Baranger model requires that we know the S-matrix (scattering matrix) elements in order either to integrate directly using or to perform the calculation of phase differences which are then integrated. Here, we begin by generating initial states which are then used to generate Moller reactant states. The Moller reactant states are propagated through an interaction under the action of the Hamiltonian and back out to calculate Moller product states and then correlation functions. We then use the Fourier transform of the correlation function to calculate the S-matrix elements which we can use in the Baranger model.

Generation of Moller states
The scattering operator S identifies how reactants |Ψ in in the infinite past map to products |Ψ out in the infinite future, The scattering operator can be defined in terms of the channel Moller operators [22,23], where, in atomic units, We can use completeness to write the incoming reactant (or outgoing product) state in the form where the |k γ γ are a separable set of reactant and product states and γ represents the full set of internal quantum states of the reactants and products [24,25]. The channel Moller operators are then used to compute reactant and product Moller states: The method we use is to begin with a Gaussian wavepacket at t = 0. We propagate the wavepacket backward as if it were a free particle for a long enough time that it does not overlap significantly with the centrifugal effective potential. We then propagate this 'intermediate Moller state' forward in time under the full Hamiltonian until t = 0. This effectively generates an intermediate state (at infinity) that would have evolved into a pure Gaussian wavepacket under no potential but that instead evolves into the relevant Moller reactant state under the full Hamiltonian of the system. Since we calculate the Moller reactant states in the asymptotic limit of the potential energy surfaces, they do not depend on the molecular state of the system but only on J and the reduced mass, μ, of the system.

Calculation of the correlation function
Having calculated the Moller reactant state, we propagate the wavepacket through the collision process to determine the Moller product state. The correlation function is a measure of the time-dependent overlap between the Moller product state and the Moller reactant state; that is, the projection of the Moller product (time-evolved) state onto the Moller reactant (initial, or t = 0) state or, in our collision process, the projection of the outbound state (the state after the collision) onto the inbound state (the state before the collision). The timedependent correlation function, C(t), has the form (in atomic units) We begin our propagation at an interatomic separation of 100 Bohr, and we consider anything farther out than 20 Bohr to be 'asymptotic' with regard to the interaction potential. However, the centrifugal effective potential reaches farther out for relevant values of the total angular momentum J, so even if we place the initial wavepacket at around 100 Bohr we still see a significant difference with J. We therefore need to generate the relevant Moller reactant states, one for each value of J, which we can use in the channel packet method [24]. Given infinite amounts of time and computational resources, the obvious method of generating a Moller reactant state would be to generate a Gaussian wavepacket starting an infinite amount of time before the collision (t = −∞) and then propagate that wavepacket until t = 0 to form the initial state. Since time and computational resources are finite, however, we must choose a suitably large time for 't = −∞' such that the Moller reactant states can be calculated in a reasonable amount of time but that the wavepacket at the time we call 't = −∞' does not overlap so much with the centrifugal effective potentials for relevant values of J that it misbehaves significantly at low kinetic energies.
In order to calculate the scattering matrix, or S-matrix, elements, we first calculate the correlation function. We propagate the Moller reactant states through the collision process to determine the Moller product states and then calculate the time-dependent correlation functions using equation (1). The wavepacket is propagated using the split operator method, in which the time evolution of wavepackets is given in atomic units by [26] We use a unitary transformation between adiabatic and diabatic representations to ensure that the potential and kinetic energy terms operate correctly.

The Hamiltonian
For the system of A 2 Π 1/2 , A 2 Π 3/2 , and B 2 Σ 1/2 states, the fullycoupled Hamiltonian is a 6 × 6 matrix [24]: The omission of the J+1 2μR 2 terms in the (2, 5), (3, 6), (5, 2), and (6, 3) elements of the third matrix has a small effect on our calculations, but it permits that 6 × 6 matrix to be approximated in 3 × 3 block-diagonal form. The remainder of our calculations will use the upper-left 3 × 3 block-diagonal portion of the Hamiltonian with an associated reduction in computational effort.

Coupled (3 × 3 matrix).
We begin from equation (3), which we call the fully-coupled 6 × 6 potential energy matrix in the diabatic representation. The first approximation we make is the omission of the J+1 2μR 2 terms as discussed above, which transforms equation (3) into block-diagonal form with two identical 3 × 3 blocks; we use the top-left 3 × 3 block with the understanding that each state is two-fold degenerate in spin. We then have for our 3 × 3 coupled potential energy matrix where Π and Σ are diabatic potentials. We use this effective potential matrix to generate coupled S-matrix elements for alkali-metal atoms colliding with noble-gas atoms. 1 × 1 matrices). We generate uncoupled S-matrix elements by making the further approximation that the off-diagonal Coriolis terms (the (1, 2) and (2, 1) elements of the 3 × 3 matrix) are zero. This allows us to diagonalize the potential matrix in terms of the adiabatic potentials of the three (now uncoupled) excited states:

Uncoupled (three
2.3.3. The ground state (1 × 1 matrix). So far we have only been concerned with excited states, but we also require calculations for the ground state of the system in order to calculate scattering phase shift differences, which in turn we need for integration within the Baranger model. The ground state does not couple to any of the excited states through collision (only through radiation), so we can express the ground state potential as a 1 × 1 matrix that can be placed in block-diagonal arrangement with the 3 × 3 matrix to form a 4 × 4 matrix. Such a matrix treatment is not really necessary, however, and we use it only for the sake of computational efficiency. The 1 × 1 ground state effective potential matrix is:

Calculation of the S-matrix elements
Once we have the correlation functions from equation (1), we can calculate the scattering matrix elements, S, by calculating the FFT of the correlation function and dividing by the channel packet expansion coefficients: This yields a scattering matrix element as a function of energy.
Since the total Hamiltonian of the system depends on J, so does the scattering matrix element. We thus have S-matrix elements as functions of J and E.

Calculation of scattering phase shifts
While one use the squared magnitudes of the S-matrix elements to calculate associated scattering cross sections [24], we use their real and imaginary parts to calculate a phase shift corresponding to each S-matrix element: Since the arctangent function is periodic, we have to check for the start of a new cycle in phase, and then add 2π to allow the total phase to accumulate. We then calculate the scattering phase shift difference between a given excited state and the ground state: Once we have the scattering phase shift difference for the entire range (in J and E) over which the collision can be said to occur, we can subtract an overall constant phase from the entire data set; here we determine an arbitrary zero of phase, in the same way that the zero of potential energy is an arbitrary choice.

Integration of the Baranger equation
Our work with the Baranger model begins with the line profile: where n is the number density of perturbers (that is, noble gas atoms) and the function g(τ ) is an integration of the singleperturber time-evolution operator over the collision phase space and the distribution of velocities [16,20,[27][28][29][30]: With an ideal gas, a Boltzmann distribution, and the Langer modification [31] for the relation between impact parameter and angular momentum , we find the broadening and shift coefficients, respectively: The left-hand side of each of equations (12) and (13) gives the broadening and shift, respectively, per unit pressure at a given temperature. Here we have reinserted as appropriate to end up with a calculation of broadening and shift in units of MHz/torr. These results give half-widths so as to be consistent with Blank [21].
We can recast these in terms of rates per concentration rather than rates per pressure. We do this by multiplying equations (12) and (13) β wpc In order to calculate the broadening and shift rates, we must perform the sums in equations (14) and (15) over kinetic energy (E) and total angular momentum (J). J is constrained to be of half-integer quantity, while E is quantized by the same energy resolution as in the calculation of S-matrix elements [24]. In this case, our energy resolution is ΔE = 2 −13 × 0.01 Hartree.
We use the ab initio potentials developed by L Blank [18,19] to calculate scattering matrix (S-matrix) elements by the channel packet method. We use those S-matrix elements to calculate the phase shift of a given state during a collision. We use the phase shifts (and, more importantly, the phase difference between a given excited state and the ground state) rather than a more direct calculation using S −1 f f S ii because calculations of the phase differences provides an intermediate check of the viability of the calculation (that is, whether the phase difference vanishes at high values of E and J). We then use the calculated phase difference θ J (E) to numerically integrate equations (14) and (15) to find the broadening and shift rates of the given spectral line. Figure 1 is a contour plot of [1 − cos θ J (E)] times the Boltzmann distribution for the temperature of interest, which gives us most of what is summed over energies and angular momenta in equation (12). Figure 2 shows a similar plot of [sin θ J (E)] times the Boltzmann distribution for the temperature of interest, which gives us most of what is summed over energies and angular momenta in equation (13). We can see immediately the effective (E, J) 'collision phase space', or region of interest in which the collision can be expected to broaden or shift the spectral line at that temperature. Here the Boltzmann distribution acts as an envelope for the function θ J (E), which does not itself depend on temperature. We can therefore use the same phase difference data for calculations in a wide range of temperatures, limited primarily by the cutoff energy defined by the limits of our phase difference data (which provides an upper limit for perturber kinetic energy and thus an upper limit for temperature).
We then sum the results shown in figures 1 and 2, as prescribed by equations (12) and (13), respectively, to find the broadening and shift coefficients at a given temperature. Now, we will consider addition schemes to account for collisional coupling between the excited states.

No addition (no collisional coupling)
Here, we assume no collisional coupling. For example, if the system is in the A 2 Π 1/2 state before the collision, it will be in the A 2 Π 1/2 state after. This is particularly effective for the A 2 Π 1/2 state, which couples negligibly weakly to the other  states during collision, but we will deal with the coupling of the A 2 Π 3/2 , and B 2 Σ 1/2 states in later sections. This allows us to compare our computational implementation of the Baranger model with the results of the Anderson-Talman model [2-4, 18, 21].
So far, we have used a Boltzmann distribution of kinetic energies, whereas the usual implementation of the Anderson-Talman model assumes all collisions occur at the thermal average velocity. If we change our implementation of the Baranger model to use a thermal average kinetic energy, the broadening and shift rates become    Figures 3 and 4 show broadening and shift rates versus temperature for each of the uncoupled states in Rb + He, and we have calculated the broadening and shift rates for the A 2 Π 1/2 , A 2 Π 3/2 , and B 2 Σ 1/2 excited states using the Baranger model and plotted those results along with equivalent calculations using the Anderson-Talman model by Blank [21].

Allard addition
Allard and Kielkopf point out that such a calculation is 'nontrivial in all but two-level atoms (atoms with only one potential difference curve, or adiabatic processes) because of fine-structure transitions between excited states that occur during the collision' [16]. That is, fine-structure mixing  produces a set of coupled equations which must be solved numerically; further, such calculations were prohibitively computationally-intensive at that time. However, this is a critical problem in an OPAL because a two-level system generally will not perform as a laser. Such effects are part of the physical processes involved in spectral line broadening. A perturber can, for instance, propagate inward (toward the emitter atom) along one potential surface, go through a transition, and then propagate outward (away from the emitter atom) along a different potential surface [16,17,31,32].
In order to account for coupling in the 2 P 3/2 manifold during collision, we perform what we call 'Allard addition' (or 'Allard coupling' [16,20]), which couples the A 2 Π 3/2 and B 2 Σ 1/2 states by weighting their phase differences equally. As a result, equations (14) and (15) are changed (only for the D 2 line; the D 1 line is calculated as in the uncoupled case) in the following way:    θ J f 1 and θ J f 2 are the uncoupled scattering phase shift differences corresponding to the two states on the 2 P 3/2 manifold (that is, the A 2 Π 3/2 and B 2 Σ 1/2 states). Figures 5 and 6 show the broadening and shift rates as functions of temperature for the D 1 lines for all nine M + Ng pairs, and figures 7 and 8 show the broadening and shift rates as functions of temperature for the D 2 lines for all nine M + Ng pairs.

Baranger addition
In the Allard addition case, the states in the 2 P 3/2 manifold are coupled in a 50/50 split, as shown by the factors of 1 2 in equations (16) and (17). We can modify the addition somewhat, to account for a variable coupling during collision; to do this, we replace the factor of 1 2 with the probability for being in each corresponding state after the collision, which corresponds to the square of each state's corresponding Smatrix element [16,20,30]. The broadening and shift coefficients that result are similar to those resulting from Allard addition, except the factors of 1 2 become variable coupling coefficients that are based on the S-matrix elements, and the phase shift differences are calculated using the 3 × 3 coupled Hamiltonian from equation (4). The broadening and shift rates using Baranger addition are: where The D 1 line broadening and shift coefficients are calculated using the coupled (3 × 3) scattering phase shift differences corresponding to the 2 P 1/2 manifold, and the D 2 line results depend on the two states on the 2 P 3/2 manifold. Here we draw a distinction between two different things both called 'coupling'; we have a coupling of the potential energy surfaces through action under the Hamiltonian (which are used to generate S-matrix elements and phase shifts), and we have coupling of phase shifts as a way of calculating what we measure as D 1 and D 2 lines instead of transitions from individual surfaces (this coupling we categorize as either 'Allard' or 'Baranger' coupling or addition). Figures 9 and 10 show broadening and shift rates as functions of temperature for the D 1 lines of all nine M + Ng pairs in the Baranger addition case, and figures 11 and 12 show broadening and shift rates as functions of temperature for the D 2 lines of all nine M + Ng pairs in the Baranger addition case. Because this form of coupling requires data about the behavior of the scattering matrix elements and not just phase differences, these results cannot be extended in the energy regime in the same way as the other cases. Thus, the broadening and shift coefficients must be limited to lower temperatures in order to prevent error due to truncation of the Boltzmann distribution at the maximum energy. In this work, we limit calculations to a maximum temperature of 500 K for the Baranger addition, whereas we extend to 800 K for the Allard addition.

Discussion
The shift coefficients are extremely sensitive to the initial Moller reactant states. This sensitivity is caused by the sine term in equation (13); for small phase shift differences, sin θ J (E) ≈ θ J (E) but cos θ J (E) ≈ 1, so small but nonzero phase shift differences cause the integrand in equation (12) to vanish but the integrand in equation (13) to remain nonzero. This nonvanishing term then multiplies the Boltzmann distribution and causes a J-independent ridge. Such a ridge appears for any nonzero offset phase as well, but a small but nonzero phase shift difference appears if the Moller reactant state generation has not propagated far enough into the distant past to escape the centrifugal effective potential. We also believe that the broadening is dominated by the close-in behavior of the potential energy surfaces while the shift is extremely sensitive to the near-asymptotic behavior of the potential energy surfaces.
To calculate the phase shift differences, phase shifts for the excited and ground states were extended linearly from the energy limits of our calculations (E = 0.0075 Hartree) to a larger energy (E = 0.012 Hartree) in order to accommodate calculations at higher temperatures. However, incorporating higher energy collisions to go to higher temperatures also requires us to include higher values of J. For example, the entire collision phase space for Rb + He at 100 K can be handled with maximum energy of 0.002 Hartree and maximum J of 65.5. Increasing the temperature to 394 K requires us to consider energies up to 0.007 Hartree and a maximum J of 110.5 to capture the entire collision phase space. Increasing the temperature to 800 K requires a maximum E of 0.012 Hartree and maximum J of 130.5 to capture the entire collision phase space. In other words, calculating at higher temperatures requires larger energies and larger values of J. We can extend phase shifts linearly in energy, but we cannot extend in J without losing critical information about that part of the collision phase space. Calculations at significantly higher temperatures will require calculations at higher values of J (and E) to capture the full collision process. Such work will be necessary to perform broadening and shift calculations at higher temperatures than about 800 K. There is good work ongoing that compares with this research [33,34].
Ultimately, agreement among broadening rates is not sufficiently good to identify conclusively which model is 'correct' for a given set of ab initio potential energy surfaces, at least at the temperatures at which experimental data have been measured. In most cases, the predictions of the Baranger and Anderson-Talman models diverge at low temperatures, so low-temperature experiments may provide a needed discriminator between the models.
A potential source of error in this implementation of the Baranger model is that the generation of the reactant Moller states does not include the off-diagonal Coriolis terms for the Hamiltonian, but include only the diagonal terms (that is, the centrifugal effective potential). As before, we can examine the reactant Moller states (at t = 0) and both the diagonal and off-diagonal Coriolis terms. We find, however, that the offdiagonal Coriolis terms do not contribute significantly to the Hamiltonian at separations as far as where we start the propagation (100 Bohr), and thus ignoring the off-diagonal Coriolis terms should not introduce a significant source of error.
There is still a great deal of theoretical work to be done in this area, from the calculation of potential energy surfaces to refinement of our scattering model and the Baranger model. Any error in the ab initio potential energy surfaces is reflected in the final results. In particular, we suspect errors in the surfaces for M + Ar because both the Baranger model and the Anderson-Talman model give results that vary significantly from experiment for these pairs. It is not clear to what degree this implementation of the wavepacket propagation technique and the Baranger model are sensitive to differences in the potential energy surfaces. It is a theoretically straightforward, but computationally intensive, process to replace the potential energy surfaces with new inputs. One could use different classes of potential, such as the Lennard-Jones (6-12) potential [35,36] that is commonly used in the Anderson-Talman model. Hager has achieved some success with the Anderson-Talman model using a 6-8 potential [37], and such a potential could be tested in the context of the Baranger model. Testing different sorts of potentials with more localized and controllable characteristics might give more information and confirmation about what parts of the potential energy surfaces give rise to which characteristics in the broadening and shift rates and intermediate calculations such as the scattering phase shift differences. There is work ongoing to generate better potential energy surfaces [38] and to utilize them to simulate reactions involving the higher manifolds [39].
The primary challenge to extending this model to predict line shape shift and broadening for other atom-atom collisions of interest is obtaining the appropriate scattering phase shifts for the collision. This requires the atom-atom potential energy curves associated with the spectral feature of interest along with the governing Hamiltonian, here given for alkali-metal noble-gas collisions in equation (3).
Another approach to generating Moller states might simply start with a Gaussian wavepacket at a very large separation distance, which could ameliorate the problem with generating reactant Moller states; in essence, the Gaussian wavepacket becomes our reactant state for which we can generate an analytic form. However, the improvement of a single problem is counteracted by the introduction of two additional problems. First, the reactant state has to be propagated through the collision process and back out to where it started; this counteracts any computational savings one might have gleaned from the lack of Moller state propagation. Second, propagating from a larger separation requires a larger computational grid in order to accommodate the space containing the wavepacket and the origin, which in turn requires FFT code capable of accommodating such a large space. This second problem might be lessened by adopting a moving reference frame that is just large enough to accommodate the wavepacket as it spreads, but we have not attempted this and we are unsure to what degree new error might be introduced through the new propagation algorithm.
Finally, we see only the distant past (or what we call the 'infinite' past) and distant future before and after the collision [24]. Because we can only look at the distant past and future, we are stuck with the impact limit of Baranger, which assumes that the duration of a collision is short compared with the time between collisions. Any work to take us out of the impact limit will necessarily involve being able to view events that occur during a collision, rather than just the distant past and future, and will require a complete reworking of the computational algorithm.