Implementing and constraining higher fidelity kinetics for DPAL models
Hal J. Cambier Timothy J. Madden







 











[]hal.cambier.ctr@us.af.mil






Kirtland AFRL









   Ionization, hydrocarbon breakdown, and other exotic processes can harm diode-pumped alkali laser (DPAL) performance and components.  We develop a physical picture of these processes, including those that drive a non-Maxwell-Boltzmann distribution of electrons, and describe an efficient approach to solve these kinetics while resolving trace species, and enforcing conservation laws.  Comparing the model to time-dependent experiments suggests that recombination and supporting processes are weaker than naïvely expected under relevant conditions, while methane seems to improve performance in the lab more than it does in the model.  Overall, this work highlights the importance of tracking the true electron energy distribution, and how incisive experiments with time-dependent driving are.  We also use the model to emphasize how ionization may pose more immediate heat loading problems in devices.




Implementing and constraining higher fidelity kinetics for DPAL models
    T. J. Madden
    December 30, 2023
======================================================================














§ INTRODUCTION



   
   
   For decades, the kinetics of alkali vapors have garnered interest given the role they play in atomic physics experiments, atomic line filters, thermionic generators, etc.  More recently, mid-infrared diodes have become spectrally narrow enough to excite individual fine-structure components of the first resonance level in alkali like potassium, rubidium, and cesium.  This allowed for a fine-structure-specific three-level lasing scheme in buffered alkali vapors, and introduced a new application: the diode-pumped alkali laser (DPAL).

   Past research into optically-pumped alkali vapors outside the DPAL context has shown it is relatively easy for the alkali atoms to reach higher excitation levels, ionize, or undergo chemical reactions, causing concerns for DPALs as well.  In earlier work, Wu<cit.> presented a model of DPAL ionization including laser rate equations, but assumed that the electrons were in a Maxwell-Boltzmann distribution with the same temperature as the buffer gas.  Oliker et al.<cit.> studied the effects of similar kinetics in a model that included computational fluid dynamics, and ray-tracing.  Wallerestein et al.<cit.> considered an even broader set of processes involving the bound electrons, and emphasized building intuition for the hierarchy of processes based on their timescales, and likely availability of reactants.

   
   However, Zatsarinny et al.<cit.>, and Markosyan et al.<cit.>, did consider a non-Maxwell-Boltzmann (non-MB) electron energy distribution through the use of finely-grained lookup tables<cit.>.  We also relax the Maxwell-Boltzmann assumption, but model, and evolve the electron distribution a bit differently, and focus on experimental comparisons.   Because processes like electron impact transitions (EIT), electron impact ionization (EII), and recombination are sensitive to electron energy, relaxing assumptions about the energy distribution avoids underestimating ionization, and its effects.  This generalization also poses numerical issues like stiffness, and maintaining positivity, which have been tackled for decades as well.  Some relevant techniques see little use in this context though, so we discuss the implementation of one approach that proved instrumental.  For the readers' benefit, we also mention issues with simple implementations of other approaches.

   The paper is organized as follows.  We give an overview of major processes, including those affecting the electron energy distribution in Section <ref>, while relegating details on inputs to Section <ref>-<ref>.  We describe the numerical implementation in Section <ref>.  In Section <ref>, we compare the model to experiments by Zhdanov et al.<cit.>.  Here, we also demonstrate, and explain the tendency for non-MB electrons to increase heat loading.  We present our conclusions in Section <ref>.





§ DEVELOPING THE PHYSICAL PICTURE



   
   A typical DPAL medium includes around 0.9 atmospheres of helium, 0.1 atmospheres of some hydrocarbon buffer, and 1-10 parts per million of alkali vapor (e.g. 10^13-10^14/cm^3), at temperatures around 400-500K.  Diodes pump alkali atoms to the upper, more degenerate, resonance sub-level, “D2”, (4 ^2P_3/2 for K).  Fine-structure mixing collisions with the buffer drive transitions to the lower resonance sub-level, “D1” (4 ^2P_1/2), which is the state used for lasing (fig. <ref>).  Higher pressure designs may achieve enough mixing with less or no molecular buffer, thus avoiding thermal-lensing, and buffer chemistry issues.  For example, Zhdanov et al. <cit.> showed that just two atmospheres of helium sufficed to mix D1 and D2 in potassium.

   
   
   

   
   Deviation from ideal operation begins when energy-pooling collisions (EP) between resonant-level alkali populate “pioneer” Rydberg levels with rates<cit.> spanning 10^-11 to 10^-9cm^3/s.  Collisional ionization (CI), (single) photo-ionization (PhI), are the fastest deleterious pathways out of these levels.  Decays (unhindered by radiation trapping) to the 5 ^2S_J and 3 ^2D_J multiplets, and energy-quenching collisions (EQ) will counter these processes.  Buffer collisions dominate EQ, and Earl and Herm <cit.> find that methane quenches the 5 ^2P_J level with a roughly 60Å^2 cross section.  For our default model, we assume this 5 ^2P_J quenching populates 3 ^2D_J and 5 ^2S_J, and that methane quenches 4 ^2D_J and 6 ^2S_J to 5 ^2P_J with a similar cross section given the similarity of all energy gaps involved (see Section <ref>).  Helium quenching becomes more significant at much higher levels, which we do not resolve in this work.  The following subsections discuss some processes, or aspects of them, that tend to receive less attention, while the appendices give more specific implementation details.

   

 §.§ Collisional and photon ionization

   
   Explicit data for potassium is scarce, but data on other alkali, and the intermolecular potential curves<cit.> indicate that Penning ionization involving a resonance-level, and a pioneer-level atom – producing just an atomic ion – forms the primary collisional ionization (CI) channel.  Looking at rubidium, Barbier et al.<cit.> found that associative ionization between one resonance-level, and one higher energy alkali atom becomes orders of magnitude weaker than Penning ionization, and they were able to explain this theoretically by invoking electronic exchange<cit.>.  This leaves collisions between one D1/D2, and one 5 ^2S_J/3 ^2D_J atom as the most plausible route for direct K_2^+ production, but processes described below can influence the dimer cation fraction far more effectively.

   
   
   Predictions by Aymar et al. <cit.>, and Zatsarinny et al. <cit.> for photo-ionization (PhI) cross-sections near D1 and D2 wavelengths imply that photo-ionization exceeds any CI channel above about 10kW/cm^2 for 10^14/cm^3 alkali densities.  Note that both the dominant CI, and PhI channels give free electrons an initial kinetic energy of 0.3-0.6eV.

   
   

   

 §.§ Elastic and inelastic energy transfer (free electrons)


   A new host of processes enter the picture once ions appear.  Fig. <ref> sketches gain and loss timescales for electrons in a typical potassium DPAL.  For the usual extent of population inversion, electrons experience net energy gain from collisions with potassium.  Inversion between D1 and D2 even gives a noticeable contribution, a fact previously recognized by Markosyan and Kushner<cit.>, and Zatsarinny et al.<cit.>.  Even if all the alkali were ionized, energy transfer via elastic collisions with helium still dominate  electron-electron energy transfer.  Energy transfer through inelastic collisions with a hydrocarbon buffer can easily surpass the energy transfer to helium as fig.<ref> illustrates (cf. Section <ref>).

   

 §.§ Hydrocarbon breakdown mechanisms


   Hydrocarbons may also present a target for electron impact breakdown, or excited-alkali reactions.  For methane specifically, CH_4 +  e^- → CH_3 +  H^- dominates electron-based breakdown<cit.>.  Data on subsequent pathways is sparse for heavy alkali, but Penning detachment predictions for other excited alkali-H^- collisions<cit.>, and measurements for oxygen anions striking excited oxygen molecules<cit.> suggest that KH+ e^- formation using the abundant K^* will have cross sections of order 10^-14cm^2.  Penning detachment using the ground state is much weaker, and has a high threshold<cit.>.

   Regarding  excited-state chemical reactions, Azyazov et al.<cit.> excited Rubidium to its second resonance level, and attempted to measure both reactive, and non-reactive quenching with various molecules.  For methane in particular, they obtained a reactive branching ratio of 0.04±0.03, and detected no RbH – the most likely product expected based on their own theoretical calculations – thus concluding  the reaction rate was negligible.   Given this limited, negative data for reactive pathways for hydrocarbon breakdown, we dot not include such processes for now.

   

 §.§ Recombination, and related dimer association rates


   Regarding recombination, the introduction references have already noted a hierarchy where dissociative recombination (DR) dominates neutral- and electron- mediated three-body recombination (NMR, EMR), and radiative recombination trails far behind.  The more energetic, non-MB nature of the electrons renders recombination – especially three-body channels – more difficult than previously expected.  Unfortunately, the dominant DR channel is still riddled with uncertainties in the energy-dependent cross section, product states, and auxiliary processes affecting the dimer ion population like association.  

   As far as product states for DR, the rules of thumb relating potential crossings to favored product states <cit.>, and analogies to experiments that do resolve products (e.g. Le Padellec et al.<cit.> for CN^+) indicate that of DR will likely be: one ground state atom, and one excited to the pooling levels, or levels just above.
   

   Arimondo et al.<cit.> examined DR for argon-buffered cesium vapor at relevant temperatures and densities, and also reported a best fit rate for association, Cs^+ + Cs + Ar→Cs_2^+ + Ar, of 2.3×10^-23cm^6/s, but noted the error might be an order of magnitude.  To our knowledge, the nearest similar measurements are for Cs^+ + Cs + Cs→Cs_2^+ + Cs at similar temperature, but lower buffer pressure, by Bergman and Chanin<cit.>, and Morgulis and Korchevoi<cit.>, who reported rates on the order of 10^-30cm^3/s, and 10^-26cm^3/s, respectively.

   

   The critical, and uncertain nature of the association rate warrants some brief technical discussion, especially to motivate the lower values that we later find necessary.  In a two-step termolecular reaction model<cit.>, the K and K^+ form a metastable ( K_2^+)^* with some rate constant, k_f.  This metastable state dissociates in a characteristic lifetime, τ, in the absence of stabilizing third-body collisions, which have characteristic rate constant k_s = n_Bσ vstab where n_B is the buffer density.  For an intermediate population in equilibrium, the effective association rate constant is:
   
   
    k_eff =  k_f k_s  1/τ + k_s .

   
   The association rate reduces to k_eff≈ k_f when stabilization proceeds much faster than metastable disintegration, or k_eff≈ k_f k_s τ in the opposite limit.  For an ion-induced dipole interaction[high ionization conceivably screens this interaction, but our estimates for the Debye length even when every alkali is ionized are on the order of 1000Å], the appropriate Langevin cross section just depends on (dipole) polarizability, α, and collision energy, E, as σ_L = π(2α/E)^1/2 (in atomic units).  Taking polarizabilities from Mitroy et al. <cit.>, leads to k_f ≈ 6× 10^-9.  Depending on whether the stabilizing atom “sees” the ion, and has a Langevin cross section to stabilize, or just a typical momentum transfer cross section, the stabilization rate will be k_s ≈ 10^-12 or k_s ≈ 10^-8.

   In order for the three-body rate constant in Arimondo et al. to be as large as 10^-23, yet have the process scale with buffer density, the characteristic τ must be 10^-3s or 10^-7s depending on which limit of k_s applies, but each value is roughly a thousand times greater than the respective 1/k_s, contradicting the condition τ≪ 1/k_s to have buffer density dependence.  For an atmosphere or so of buffer, the k_f ≈ 6× 10^-9 estimated above coincides instead with a “three-body” rate around ≈10^-26.

   
   

   

 §.§ Miscellaneous processes


   Dissociative collisions by excited alkali, and electron impacts – called dissociative excitation (DE) – can counteract association.  We model the former process based on measurements for sodium by Tapalian and Smith<cit.> (Section <ref>).  Peak cross sections, and energy dependence for DE are fairly uniform among a variety of cation dimers<cit.>), so we also included a toy model of DE (Section <ref>), but found that it played a weak role.

   
   Neutral dimer association can divert population from the lasing cycle, counteracted by dissociation due to excited alkali<cit.> (Section <ref>).  The neutral and ion dimer populations can interact directly through charge exchange, but at the moment we do not model this.








§ IMPLEMENTING THE MODEL




   Now we must simulate all the physics discussed above.  Processes involving just the thermalized, massive particles follow standard rate equations.  For the basic three-level cycle, laser intensity, and pump term, we adopt the model of Hager and Perram<cit.>, which we summarize below.  To represent the non-MB electrons, we bin the spectrum to obtain ordinary differential equations, but this introduces issues with energy conservation which we also address here, following which we explain how to resolve these numerical issues.

   We model the core laser kinetics by source contributions (Δ s_x) to densities of the first three levels, n_i (starting at i=1 for the ground state), and the two-way laser photon density, ψ_L:
   
   
    Δ s_1    =  σ_31(n_3-2n_1)ω + σ_21(n_2-n_1)ψ_L + n_2 Γ_21 +n_3 Γ_31, 
    Δ s_2    = -σ_21(n_2-n_1)ψ_L - n_2Γ_21 +γ_mix(n_3-2exp[-θ]n_2), 
    Δ s_3    = -σ_31(n_3-2n_1)ω - n_3 Γ_31 -γ_mix(n_3-2exp[-θ]n_2), 
    Δ s_L    = (rt^4exp[2σ_21(n_2-n_1)l_g]-1)ψ_L/τ_RT + {f× n_1 A_10/τ_1}

   
   where
   
   
    ω = I_P,in/hν_Pσ_31(n_3-2n_1)l_g(exp[σ_31(n_3-2n_1)l_g]-1) t_P (1+t_P^2 r_P exp[σ_31(n_3-2n_1)l_g])

   
   is the absorption rate based on input pump intensity I_P,in.  The g_i are level degeneracies, A_ij; radiative transition rates, Γ_ij; total decay rates, γ_fs; the fine-structure mixing rate, τ_RT; the cavity round-trip time, l_g; the gain length, r_x the output coupler reflection coefficients, t_x; the intracavity transmission coefficients (set to 1 throughout), θ=(E_3-E_2)/(T), where T and  are temperature and the Boltzmann constant.  The last term in the laser equation is a rough estimate of the spontaneous emission seed for lasing, which we retain since our method is still based on evolving the equations in time, whether we apply to finding steady-state solutions or not.  The laser term contains a potentially very large exponential factor times the inverse cavity timescale.  We can address this extreme stiffness by keeping the overall effective timescale above 0.01τ_RT via
   
   
    dψ_Ldt→ a^-1tanh[a(r_L exp[2 σΔn l_g] - 1 )]

   
   with a=0.01, for example.  Testing with different a, and the raw equations showed that this still preserved behavior on timescales of interest while avoiding situations where the code had to switch to excessively small timesteps.  We also track contributions to the general buffer specific energy density (the thermal bath) due to mixing of the D1 and D2 energy levels (E_2, E_3):
   
   
    Δ s_eb = γ_mix (n_3-2 n_2 exp[-θ]) (E_3-E_2)

   

   
   Note that energy exchange with the thermal bath also gets tallied for the other processes (like pooling, and quenching).

   Ensuring global energy conservation with a binned electron spectrum, and irregular grid, is not trivial.  Solutions include: modifying the target-bin branching ratios to conserve energy at the cost of straying from the discretized differential cross section, and adopting a better, higher-order “sub-grid” model than a flat-top distribution (see e.g. Le and Cambier<cit.>).  The sub-grid approach avoids unsavory tweaks to the cross sections, and better handles transitions with energy less than a bin width.  For now though, re-weighting is simpler to implement, and test various methods with, but the higher order discretization warrants further consideration.  This accounts for ionization, recombination, and inelastic collisions for free electrons, which just leaves energy transfer with the buffer.  We base our formulation on the Lorentz model with a correction for the finite temperature of the buffer particles, (see e.g. Loureiro & Amorim <cit.> where it is given in terms of velocity magnitude)
   
   
    ∂_t n_E = ∂_E[ E-T/2τ_exch(E) n_E + E Tτ_exch(E)∂_E n_E ],

   
   where E is electron energy, n_E the electron differential energy density, and τ_exch the characteristic energy exchange timescale.  We implement the discretized version as a series of upwinded advection fluxes between neighboring bins, plus a diffusion term, so the resulting source terms are more automatically conservative.  See Section <ref> for details.  
   



   

 §.§ Conservation-enforcing projection ala Sandu (2001)

   


   Resolving trace populations while avoiding negative densities and conservation error is a familiar challenge, motivating various approaches over the years.  For example, Preussner and Brand<cit.>, and Bertolazzi<cit.> both focus on preserving non-negativity with semi-implicit and implicit methods respectively.  Instead of focusing on the integration method, Sandu<cit.> remarks that finding a point obeying conservation rules, closest to some arbitrary method's guess, just defines a linearly-constrained quadratic optimization problem.  This approach is highly general, avoids longterm drift in error, and forms the backbone of our code.

   The projection relies on expressing conservation laws as a linear constraint:
   
   
    A^⊺ n' - A^⊺ n_0 - x = 0

   
   where n_0 is some previous, trusted population vector, and x=0 absent any external source/sinks.  Sandu reviews how to find A given any stoichiometry, but alkali number and charge conservation rows can be written by inspection for our problem:
   
   
      A_N^⊺    = (   1,,1,            2,         1,         2,         1, 1,     0,   ),       A_N^⊺ n = n_alk(t_0) 

      A_Q^⊺    = (   0,,0,            0,         1,         1,        -1,-1,     0,   ),       A_Q^⊺ n = 0 

                                                                       
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   

   In principle, neutral buffers will contribute their own rows, and internal energy conservation can be included in a similar fashion.  In practice, buffer conservation holds well either way, while adding internal energy means adding a large external source term (the pump).  For large timesteps, the energy conservation correction tends to interfere with the others, so some numerical error from the (time) integration method can still affect energy for now.  This is not an issue for steady-state-seeking simulations, as any equilibrium is conservative by construction of the source terms.

   This A then features in the optimization problem
   
   
    n^n+1 = argmin(12 n^n+1 (G  n^n+1) - (G ñ^n+1)^⊺ n^n+1) : A^⊺ n^n+1=y_0,  n^n+1≥.

   
   where ñ^n+1 is the solver's uncorrected guess, and G is the error metric for the numerical solver.  In Sandu's example G is a diagonal matrix
   
   
    G_ii = 1/(N_species(tol_abs+tol_rel|n_i|)^2).

   
   The weighting assumes extra importance in a problem like ours with many levels of the same atom (i.e. a long row of zeros in A^⊺).  Without any for example (G_ij≡δ_ij), the sea of trace populations will ebb and flow when the projection corrects changes in large, dynamic populations.

   This tool frees us to take steps Δ n using a simple integration method like linearized implicit Euler:
   
   
    (I-Δt J)Δ n =  sΔt

   
   where s is the total source vector, and J the Jacobian.



   

 §.§ Alternative approaches

   


   We consider it worth mentioning here some alternative approaches, and immediate issues we faced applying them.

   Computational singular perturbation (CSP)<cit.> takes the opposite approach: instead of addressing accuracy issues for implicit methods, it addresses stiffness issues for explicit ones.  It achieves this by breaking up the source term into modes with characteristic timescales, typically found by straight-forward eigen-decomposition, or block diagonalization of the Jacobian, with a possible higher order correction for evolution of the bases.  `Fast' modes can reach a quasi-equilibrium with respect to the slow ones (or `exhaust').  For example, the laser levels in our problem often reach a quasi-equilibrium that evolves slowly with respect to the slow changes caused by pooling, etc.  CSP works well on the `conventional' and three-level kinetics of our problem, but transfer between the electron bins alters modes on a fast (Courant-Friedrichs-Lewy) timescale, generating fast modes which do not exhaust for a long time without the (expensive) higher order corrections.  Simple attempts at splitting off processes that only shift electrons in energy helped little as the CSP modes still changed drastically between steps.  Efficient application of CSP would likely require representing the electron spectrum differently, e.g. via moments, or some other coarser expansion.

   Applying the semi-implicit algorithm in Preussner and Brand is straightforward, and should complement the projection step nicely: it avoids negative densities that increase the cost of the projection, while the projection should reduce the need for very small steps to ensure conservation over long times.  However, this approach still required small timesteps to avoid oscillations associated with the radiative component, and to a lesser extent, electron drag and excitation/de-excitation.  Given the lower cost per step, some variation that splits off the radiative component may still offer a practical path forward for larger, more expensive engineering simulations




§ APPLYING THE MODEL




   Knowing the physical ingredients and main uncertainties, numerical techniques and issues, we can now compare predictions to experiments in a time-dependent, ionization-prone regime, as well as draw some general lessons from simulating general, steady-state situations.


   

 §.§ Testing against time-dependent experiments

   


   To avoid complications from thermal build-up and other slow processes, Zhdanov et al. <cit.> ran a set of experiments with time-varying, sub-millisecond pump pulses for a potassium vapor at 190^∘C with varying buffer composition.  Their main experiments used 500 Torr of buffer with a varying percentage of methane, and a peak pump power of 160W, translating to an intensity around 30kW/cm^2 for their reported beam profile.  They reported laser output, as well as 5 ^2P_J-4 ^2S_J fluorescence, but not in absolute units.  They also measured 5 ^2P_J fluorescence for trials with 200 Torr of pure helium, and pure argon, and saw that neither composition lased.

   We first ran the model with just the basic three-level model turned on to make sure it reproduced the lasing threshold, using the reported pump linewidth, laser beam profile, gain and cavity lengths, and pressure-broadening of the D1, and D2 transitions based on Pitz et al.<cit.> (c.f. Section <ref>).     
   Initial runs with the default kinetics, and slight variations, established some basic tenets of an `alkali-depletion' paradigm under these conditions:
   
   
      
   1   past ∼1-10% (alkali) ionization, impact ionization and DR dominate free electron population gain and loss
      
   2   this leads to sigmoidal growth of ion fraction, where drag, and pump intensity can also skew its shape, saturation timescale, and final ion fraction
      
   3   at late times, the barrier to more ionization is re-energizing low energy, post-impact electrons back above ionization thresholds before they recombine
      
   4   increased buffer drag delays and diminishes the peak in 5 ^2P_J population matching the trend of their 200 Torr experiments, but the early growth is never so linear in time as seen in the experiments

   

   
   
   

   The simulations in fig. <ref> illustrate the basic depletion mechanism, where the association rate for K_2^+ production was set to the equivalent of a 10^-26cm^6/s three body rate, and atomic-level quenching was reduced to 20% of its default value.  Since we assumed it most significantly quenches the 5 ^2P_J, 6 ^2S_J, and 4 ^2D_J levels, further reduction of the DR cross section, or faster impact ionization out of these levels would generate similar behavior.

   
   
   

   Fig. <ref> shows the evolution of the electron energy distribution in the pure helium run, a spectrum for one of the methane-added runs, and a buffer-temperature Maxwell-Boltzmann distribution with the same final electron count as the pure helium case.  The non-MB curve is roughly one-tenth the MB curve below 0.1eV, so already the total recombination rate is significantly changed.  The curve for the methane case highlights methane's ability to suppress not only number, but mean electron energy.  Figures <ref>-<ref> break down various contributions to the source terms for the ground state, D1, and D2, which demonstrates the role of impact ionization here, and explains a strange phenomenon described below.

   
   
   

   
   
   

   
   
   

   For pure helium, the rise in laser output slightly overshoots the other mixtures before cresting in fig.<ref>, and the source term figures reveal how moderate ionization briefly boosts laser gain to cause this.  From lasing onset – seen as a spike around 20-30μs – to 50μs or so into the simulation, the sum of three-level, impact transfer processes, and eventually association into K_2^+ all drive significant ground state depletion.  Meanwhile, D1 experiences net growth, thus raising laser gain until severe ionization depletes both populations.  Net electron impact mixing of the D1, and D2 levels also switches sign around this time, as D2 atoms sufficiently outnumber D1 atoms.

   Artificially raising absorption cross sections by a factor of 5 dramatizes this effect in fig.<ref>.  This also leads to delayed, saturated curves for the methane mixtures resembling their output curves in the experiment, but if this mechanism does play a role in explaining the responses with methane, it can not be the sole cause.  For now, this overshoot serves as a hint for where we should examine the physics and methods further, e.g. perhaps the modeled electrons are too hot at early times leaving excitation too strong versus de-excitation, or perhaps the K_2^+ association rate is still too high.

   
   
   


   

 §.§ Ionization can drive significant heat loading

   


   Predictions of vigorous heat generation by free electrons in a pumped, alkali vapor date at least as far back as a paper by Measures<cit.>.  Additional, unexpected heating can hurt beam quality, among other effects, but this aspect of ionization receives relatively less attention, so we review some simple arguments behind strong heating.

   From fig. <ref>, one sees that a typical electron born around 0.3-0.5 eV, will frequently gain 1.6eV by de-exciting a K(4P) atom on a sub- timescale.  Often after one such jump, energy loss to the buffer is suddenly much faster.  As ionization and recombination are orders of magnitude slower, this means a typical electron can convert a pump photon worth of energy to the buffer hundreds of times while free.  

   Meanwhile, bound electrons converting the ≈ 0.01eV D1-D2 gap to heat at a 1-10 GHz rate via fine-structure mixing present the greatest heat source in ionization-free models, and the logical point of comparison:
   
   
    f_ioni n_A 
          n_A .

   
   This shows the drag heat load approaches the inherent mixing one around 10% ionization.

   Electron energy loss to methane will largely go into local heat as well: Menard-Bourcin et al. <cit.> give cross sections for vibrational-translational energy transfer for methane in helium, and the resulting rates exceed radiative decays based on Yurchenko et al. <cit.> by several orders of magnitude.  The code accounts for this.

   Fig. <ref> shows heat load, and fraction of it from mixing, helium energy transfer, methane energy transfer, atomic level quenching, and other sources, for a small grid of alkali density, and methane buffer fraction at a fixed value of angle-averaged and frequency-integrated pump intensity, J̅ of 100kW/cm^2.  At moderate alkali density, adding a small amount of methane drops ionization and associated heat load substantially.  For the same total buffer density at higher alkali densities, methane does little besides replace helium as the pathway for heat generation.  Total buffer density must be raised for such conditions.
   
   
   





§ CONCLUSIONS



   Active DPAL conditions fundamentally lead to a non-Maxwell-Boltzmann distribution for free electrons, and ignoring this effect underestimates the propensity for alkali to ionize, as well as the extra heat load associated with ionization.  By comparing the model to experiments performed by Zhdanov et al., whose variations in time and composition probed the transition from good to poor performance, we confirmed this basic point, but also showed important discrepancies.  The model predicted a slight, temporary boost to laser gain at early times for the pure helium case, in contrast to the experiment.  The model also over-predicted the efficacy of methane.  The latter is less surprising in the face of quenching uncertainties, while the former constrains more central aspects of the model, like impact ionization, drag, alkali dimer ion association. 

   To reduce modeling uncertainties, we have started implementing the higher-order discretization of the distribution, obtaining more accurate impact ionization estimates, and tying the higher fidelity kinetics to a more detailed spatial model of the laser.  We also intend to investigate three-body recombination with non-MB electrons more rigorously.

   The comparison illustrates the utility of time-dependent pumping on kinetics-relevant timescales for setting stronger constraints on models.  Further variations on this idea, as well as more quantitative diagnostics for the basic three levels and excited states, are warranted.  Direct measurements, or higher-quality estimates of the most uncertain processes would obviously be most helpful too.









   We thank Dr's. Habib Najm, John Shadid, and their colleagues for helpful discussions, especially regarding CSP.  HJC thanks Ben Oliker for fielding basic questions at the start of the project.  We also thank an anonymous referee for detailed feedback, which significantly improved the paper.










   

§ NEUTRAL AND HEAVY ION KINETICS

   
   

   We included all individual fine-structure levels from 4 ^2S_1/2 through 6 ^2S_1/2, and used radiative transition rates from calculations by Nandy et al.<cit.> except for D1 and D2 transitions to ground.  Processes involving free electrons do not have rates, but energy-dependent cross sections.  Details, or references appear in subsequent appendices where not already given.  

   For the basic three levels, the mixing rate γ_mix was based on
   
   
    γ_mix = (n_He Q_He + n_CH_4 Q_CH_4) ⟨ v ⟩_th

   
   with Q_He = 18.7Å, Q_CH_4 = 58.9Å, and ⟨ v ⟩_th the mean thermal velocity.  Broadening rates for D1, and D2 absorption were taken from Pitz et al.<cit.>:
   
   
    Γ_D1/ = 19.84 P_He + 13.08 P_CH_4

   
   
    Γ_D2/ = 19.84 P_He + 27.78 P_CH_4

   
   We used these to calculate line-center cross sections, and overlap with the pump spectral profile, which Zhdanov et al. said was a Gaussian with 12 GHz full-width half-maximum.


   As noted, our pooling rates our informed by the results of Namiotka et al.<cit.>.  We did not make any assumptions about hybrid-J (e.g. 4 ^2P_3/2+4 ^2P_1/2) pooling rates, but know from other alkali<cit.> that they are likely comparable, so we expect the model to underestimate pooling.
   
   


   The Earl and Herm<cit.> quenching cross section for 5 ^2P_J by methane was not fine-structure specific, nor was the destination state resolved.  For now, we just took a similar multiplet-to-multiplet cross section for the neighboring states with similar energy separation, and divided it by the number of levels involved, and scaled the rates up or down all together.  For the simulations shown, we used the cross sections in table <ref>.

   

   As discussed above, the forward dimer ion association rate was still modeled as an effective three-body rate
   
   
    K^+ + 4 ^2 S_1/2 + He k→ K_2^+ + He .

   
   We took the neutral association rate to be
   
   
    2× 4 ^2 S_1/2 + He k→ K_2 + He .

   


   Collisional dissociation by excited atoms, again based on Tapalian and Smith<cit.>, and Ban et al. <cit.>, was modeled as
   
   
    4 ^2 P_J +  K_2  σ v→ 3 × 4 ^2 S_1/2

   
   
    4 ^2 P_J +  K_2^+ σ v→ 2 × 4 ^2 S_1/2 +  K^+

   

   Reverse rates are calculated by detailed balance, where Tango et al.<cit.>, and Magnier et al.<cit.> provide details necessary to compute the molecular partition functions, as well as heat taken from, or put into, the thermal bath for the participants on each side to be fully thermalized.



   

§ ELECTRON/NEUTRAL-MEDIATED RECOMBINATION (EMR/NMR)

   

   As discussed in <cit.>, three-body recombination involves the probability an electron is close enough to an ion to be captured while colliding with a third body to remove energy and momentum:
   
   
    s_3BR∼(n_e 4π3 r_c^3) (n_+ n_3 ⟨ (π r_t^2) v_rel⟩)

   
   where r_c = e/⟨ E_e ⟩ or e/kT for a thermal population (in cgs units), and π r_t^2 is the energy-momentum transfer cross section.  The r_t equals r_c if the third body is a thermal electron.  This simple formalism applies most literally in the case of NMR at low density where third body impacts almost immediately re-ionizing the atom are unlikely.  Otherwise, the same impact ionizations counteracting electron-mediated recombinations should counteract neutral-mediated ones.

   


   

§ ELECTRON IMPACT ENERGY TRANSFER (EIT)

   


   Phelps et al. <cit.> measured cross sections for many impact transitions from the ground state, and the usual Klein-Rosseland relation gives the corresponding de-excitation cross sections in terms of lower/upper level degeneracies g_l/g_u, and transition energy;
   
   
    σ_dex(E) = g_lg_uE+E_tEσ_exc(E+E_t).

   
   Note that the data were not fine-structure specific, so we assumed the branching ratios were purely based on degeneracies.  Other intermultiplet transitions were estimated based on the model in Vriens & Smeets <cit.>.  While one-electron, intramultiplet impact transitions are dipole forbidden, two-electron processes - with a threshold - are not, and Moores & Sheorey <cit.> calculated such a cross-section for K(4 ^2P_J) mixing.

   Regarding methane, Itikawa et al. <cit.> provide mode-specific cross sections for the first vibrational modes, which are significantly larger than cross sections to excite rotational modes<cit.>.  We do not track individual levels of methane, as they will remain Boltzmann-distributed far more easily than potassium, so we use this assumption, plus the mode information (energies, degeneracies) in Itikawa et al. to include de-excitation from thermally excited methane.




   

§ ELECTRON IMPACT IONIZATION (EII)

   

   Runs discussed here used a total cross section with the incident energy scaling derived by Gryzinski <cit.>:
   
   
    σ^G_EII = π⟨ r ⟩^2 1x(x-1x+1)^3/2(1+23(1-12x)log[2.7+√(x-1)])

   
   where x is the ratio of incident kinetic to ionization energy, E_inc/W_i, and ⟨ r ⟩^2 the squared mean radius of the outermost electron(s).  

   To simplify matters, a similar probability distribution for the secondary electron's energy based on <cit.> was used for whatever total cross section was chosen above - as opposed to a certain differential cross section -  
   
   
    𝒫(E_s) = 1/(E̅arctan((E_p-W_i)/2))1+(E_s/E̅)^2

   
   where E_s, and E_p are the secondary's and primary's (initial) electron energies, and E̅ is a parameter in their model set to 0.8W_i for results presented here.  Since a source bin can map to many target bins in EII, the adjustment for global energy conservation in this case is done by scaling and skewing the entire product distribution slightly, leaving a linear system with two constraints (number and energy conservation), and two unknown parameters to solve.



   

§ DISSOCIATIVE RECOMBINATION AND EXCITATION (DR, DE)

   

   DR cross sections tend to share a power-law “envelope” inversely proportional to energy, and values up to 10^-15cm^2 at 0.1 eV, but they can also fluctuate wildly in strength, and product branching ratio (see e.g. Little et al.<cit.>).  Some of the reduction that better reproduces data, may really reflect steep, sporadic drops instead.  Once more urgent issues are addressed, examining the effect of toy notches versus overall re-scaling may be worthwhile.

   The explicit formula for “default” DR was σ(E) = 10^-15cm^2 (E/eV)^-1, while the reduced-recombination runs used 0.05 times this, and the product branching ratio was a constant 50/50 split between 6 ^2S_J and the 4 ^2D_J multiplet.

   The toy model we used for DE was another power-law with sudden activation at the dissociation energy of K_2^+, so σ_DE = 3×10^15cm^2 (E/eV)^-1 : E>0.827eV.

   



   

§ DRAG

   
   

   For clarity, we spell out how we discretize drag.  At every wall indexed j, between cells indexed j-1 and j, the source contributions to bin densities, and the buffer thermal energy pool are:
   
   
    Δ s_j = -f[E_j-T/2τ_jn_*Δ E_*] - 2[E_j Tτ_jn_j/ΔE_j-n_j-1/ΔE_j-1E_j+1-E_j-1]

   
   
    Δ s_j-1 = f [E_j-T/2τ_jn_*Δ E_*] + 2[E_j Tτ_jn_j/ΔE_j-n_j-1/ΔE_j-1E_j+1-E_j-1]

   
   
    Δ s_eb = 12( E_j+1-E_j-1) Δs_j-1

   
   where f=1, and n_*=n_j for E_j>T/2, and f=-1, n_*=n_j-1 for E_j<T/2.  The local energy exchange rate with the buffer, τ_j^-1 = (m_e/m_He) n_Heσ_mt(E_j) v_e(E_j), where σ_mt is the momentum transfer cross section, n_He the buffer density, m_He the buffer particle mass, m_e the electron mass, and v_e the electron velocity.







