Hybrid-PIC modeling of laser-plasma interactions and hot electron generation in gold hohlraum walls

Indirect-drive thermonuclear-ignition experiments¹ are presently being conducted on the national ignition facility (NIF).² These experiments use the x-ray heating caused by the deposition of laser energy into a blow-off plasma on the walls of a hohlraum to drive the implosion of a DT-fuel capsule.³ While progress has been made in understanding these experiments since NIF became operational, substantial disagreement remains between simulated and measured neutron yields.⁴ The interaction of the NIF lasers with the inner surface of the hohlraum can generate a hot non-thermal electron population, as well as thermal x-rays. The hot electrons can be produced from multiple laser-plasma interactions (LPI) including parametric instabilities, such as stimulated Raman scattering (SRS) and two plasmon decay (TPD).⁵ Hot electrons generated in this way can penetrate the capsule and deposit unwanted preheat.⁶ Presently, simulations of NIF shots are conducted using radiation-hydrodynamic fluid codes. Fluid codes generally model the collisional absorption of laser energy in the hohlraum by a ray-tracing algorithm⁷ based on the eikonal approximation (⁠ $λ ≪$ relevant lengthscales). However, parametric instabilities such as SRS and TPD cannot be modeled by this approach, as they require a full electro-magnetic (EM) treatment of the fields. Moreover, the non-thermal effects on the electron energy distribution (EEDF) cannot be captured by a fluid treatment.

In this paper, we describe our initial attempts to use a hybrid particle-in-cell (PIC) code to study the problem of parametric instabilities in the blow-off plasma in the hohlraum walls. We have used the code “Chicago” for the work presented.⁸ Chicago is a fully relativistic, fully electromagnetic, massively parallel, particle-in-cell code that is capable of conducting 1D, 2D, and 3D simulations. Chicago can conduct fully kinetic calculations as necessary and fluid-particle simulations when such calculations are accurate. The hybrid implicit PIC algorithms in Chicago make it possible to conduct, with modern computational resources, physically realistic simulations of hohlraum blow-off plasmas. Because typical PIC constraints such as the need to resolve the plasma frequency and Debye length are relaxed in the implicit algorithm, simulations of such plasmas are possible even for spatial scales of several millimeters and time scales of tens of nanoseconds.

The NIF lasers, with a wavelength $λ ∼ 1 / 3$ μm, interact with the high-Z (usually gold) hohlraum walls for timescales on the order of a ns, during which time the blowoff plasma at the wall can hydrodynamically expand into the interior over distances on the order of 100 μm. On such time and lengthscales, although suitable for a hybrid approach, it is not feasible to resolve the laser frequency and wavelength on a simulation grid, as is required to observe the growth of processes such as SRS and TPD. For this reason, we employ a two-tiered approach to hohlraum modeling. An initially thin blowoff plasma is assumed to be present at the walls. Large-scale simulations (⁠ $Δ x > λ$⁠) of the hohlraum wall can be conducted using a fast quasi-neutral PIC-based magneto-hydrodynamic (MHD) algorithm⁹ with “fluid” species, in which Lagrangian ion macroparticles are advanced by momentum equations based on fluid equations, in conjunction with a generalized Ohm's law. In a recent publication, the fluid algorithms in Chicago have been benchmarked¹⁰ against the fluid code Hydra¹¹ for several test problems with good agreement found between the pure fluid code and the hybrid-PIC approach. Moreover, the MHD algorithm in Chicago was found to agree well with experimental data in simulations of a dense plasma focus.¹² The deposition of laser energy into the wall is modeled by tracking macroparticle photons which obey ray-tracing equations of motion. From such simulations, we can observe the time evolution of the hohlraum walls and characterize the density gradients and other plasma conditions in locations in which LPI instabilities are expected to occur. From these results, we can transition to smaller-scale highly resolved (⁠ $Δ x ≪ λ$⁠) simulations using traditional kinetic PIC methods, from which we can fully model all of the non-linear laser-plasma interactions, as well as assess the details of the EEDF. Chicago has also been used extensively to study LPI, such as beatwave generation^13,14 in underdense plasmas, with good agreement found between simulation and theory. For notational purposes, simulations of the hohlraum walls using the hybrid PIC method will be referred to as “hohlraum” simulations and the highly resolved kinetic PIC simulations as the “LPI” simulations.

From the hohlraum simulations, we can observe the time evolution of the walls. The TPD and SRS instabilities are expected to occur at or near the quarter-critical surface in the plasma, the surface where the local plasma density is one quarter of the critical density associated with the laser wavelength. Approximate theoretical predictions of the linear growth rates of these instabilities, with effects due to both collisions and plasma inhomogeneity included, which have been compiled by Wen et al.,¹⁵ can be used to assess whether the walls in the hohlraum simulations would be expected to exhibit instability growth. According to such estimates, laser intensities $≳ 10 15$ W/cm² may be sufficiently large to introduce the TPD instability in gas-filled hohlraums due to the detailed nature of the time evolution of the quarter-critical surface. In vacuum hohlraums, by contrast, the stabilizing effect of collisions should be sufficient to suppress instability for $I ≤ 10 15$ W/cm².

In 1D and 2D LPI simulations, which use initial conditions based on the hohlraum results, we can clearly observe the SRS and TPD instabilities, as well as collisional absorption. We find that SRS and TPD intensity thresholds calculated by simulation are in good agreement with theoretical estimates for the intensity thresholds given in Eqs. (2) and (3) of Sec. II, which justifies the use of the theoretical thresholds to estimate susceptibility to instability in the hybrid simulations. We have observed hot electron generation due to these instabilities. From this, we have concluded that the gold hohlraum walls in gas-filled hohlraums can be TPD unstable at relevant intensity levels, and could be a source of hot electrons. SRS, if unstable, can also generate hot electrons. But our analysis indicates that SRS threshold intensities are generally higher than those for TPD, and that SRS is likely stable for NIF relevant intensities in the parameter regimes we consider. The possible effects of hot electrons generated by such means on an inertial confinement fusion target are outside of the scope of present paper.

Significant modifications to the Chicago code were required to carry out this work. Since the NIF hohlraum often contains a background gas fill,¹ the single-ion MHD model discussed in Ref. 9 was extended to allow for multiple ion species. Similarly, the ray-tracing capability required to model collisional absorption of laser energy was also introduced into the code. These algorithms are discussed in detail in  Appendixes A and  B, respectively. Moreover, the theory of the SRS and TPD instabilities is reviewed in  Appendix C. Hohlraum simulation results are discussed in Sec. II, with 1D and 2D LPI simulation results discussed in Secs. III and IV, respectively. Section V presents conclusions and possible future directions.

Results from 1D gold hohlraum simulations are discussed in this section. The hohlraum wall is modeled as a gold plasma with the initial density profiles of gold ions and electrons shown in the top left of Fig. 1. Gold ions in the wall (⁠ $x ≤ 0.1$ mm) have a density of $5.9 × 10 22$ cm⁻³ which corresponds to solid density. In the initial simulations, we assume a vacuum hohlraum interior. Simulations with a background He gas are discussed later. A 5 μm density ramp is applied to the gold profile at the vacuum interface. Although the initial ramp length is short compared to its final length after being exposed to a several ns pulse, it is still relatively long compared to the wavelength of the laser, $λ = 1 / 3$ μm. So the use of the eikonal approximation to model collisional absorption should be well justified. The gold plasma is given initial temperature values of $T e = T i = 0.15$ eV. A laser with intensity $I = 10 15$ W/cm² is injected from the right hand side of the simulation space. At a laser wavelength of 1/3 μm, the critical density is $n c = 10 22$ cm⁻³ (see  Appendix B for details). A 0.1 ns linear temporal ramp is applied to the electric field amplitude of the laser pulse.

For the vacuum hohlraum simulations, the MHD algorithm is used in concert with the Eulerian fluid algorithm.¹⁶ An MHD formulation suitable for use in PIC codes has been described in Ref. 9. This algorithm has now been extended to allow for multiple ion species. Some details on this extension are given in  Appendix A. The Eulerian fluid algorithm allows simulations with ∼1 particle/cell/species. This feature, combined with the increased timestep allowed by the MHD approximation, allows for relatively fast modeling of dense plasmas over hydrodynamic timescales. Equation of state (EOS) data are obtained from tables pre-calculated by the Propaceos code.^17,18 The electron-ion collision frequency is obtained from the Lee-More model with Desjarlais corrections (LMD),^19,20 and we use a single-group (gray) radiation solve^16,21 to calculate the radiation transport in the hohlraum. We use a variable cell size with a minimum value of $Δ x ∼ 0.5$ μm in the region of the density ramp. Laser propagation and energy deposition are modeled by a ray-tracing algorithm⁷ which has been implemented in a form suitable for use in PIC codes, in which macroparticle “ray” photons are propagated through the simulation space. With this algorithm, which is described in  Appendix B, collisional absorption of laser light in the hohlraum walls can be modeled without requiring spatial and temporal resolution of the laser wavelength. The timestep is limited by the Courant condition, as ray photons cannot cross a cell in a single timestep. Only a few ray photons need to be injected per timestep to get convergent results for the heating of the plasma due to collisional absorption.

We define the critical point, x_c, by the equality $n e ( x c ) = n c$⁠, and the quarter-critical point, x_q, where $n e ( x q ) = n c / 4$⁠. The density inhomogeneity lengthscales are also defined (in 1D) as

where i can be either c or q. Snapshots of the gold ion and electron density, shown in Fig. 1, demonstrate that x_c and x_q move into the vacuum on the order of tens of microns for a pulse of 2 ns, along with shock compression at the original plasma-vacuum interface. The radiation temperature quickly reaches a steady-state value of $∼ 200$ eV in the optically thin “quasi-vacuum” region, which we define roughly here as the region, $x ≳ x c$⁠.

To quantitatively describe the blow-off plasma properties in order to assess the possibility of issues with LPI instabilities, we have tabulated the instantaneous positions of x_c and x_q, along with the density inhomogeneity gradients, L_c and L_q, which can be calculated by Eq. (1) using the local density slopes at the critical points. The results for the critical surfaces, x_i and L_i, along with some relevant local plasma parameters are plotted as a function of time in Fig. 2. During the 2 ns simulation time, the critical surfaces both increase roughly linearly with time and move outward at a velocity on the order of 100 μm/ns. The critical lengthscale, L_c, increases roughly linearly with time after about 1 ns, while the growth of L_q rises rapidly early (t < 0.5 ns) but increases much more slowly at later times. The local value of T_e and $Z ¯$ at x_c and x_q are also plotted as a function of time in the figure. The behavior at x_q is qualitatively consistent with theoretical models of laser-driven ablation in plasmas in which absorption is dominated by inverse bremsstrahlung, which predict a self-similar isothermal expansion of the ablation region in which L_c increases linearly in time with a slope given by the local ion sound speed.^22,23 Such a quasi-steady-state situation, $L c ∝ t , T e ( x c )$ and $Z ¯ ( x c )$ constant, is seen to be established in Fig. 2 for t > 0.5 ns.

Using the data compiled in Fig. 2, we can now estimate whether the simulated hohlraum is susceptible to the SRS or TPD instabilities.⁵ In SRS, the incident EM wave from the laser decays into a scattered EM wave and an electron plasma wave (or plasmon), while in TPD the EM wave decays into two plasmons. Based on frequency and wave number matching conditions, we find that SRS can occur for $n e ≤ n c / 4$⁠, and TPD when $n e ∼ n c / 4$ (see  Appendix C for details). In the LPI simulations discussed in Sec. III, we find a strong back-scattered EM signal and a localized plasmon signal near the quarter-critical surface, $n e ∼ n c / 4$⁠. In this case, the angular frequency of the laser is approximately twice the local plasma frequency, ω_pe, and the scattered wave is at the plasma frequency. The back-scattered wave is consistent with the SRS instability, and the localized plasmon signal is consistent with both SRS and TPD. The maximum SRS growth rate at the quarter-critical density, including the effects of both plasma inhomogeneity (finite L_q) and electron-ion collisions, can be written^15,24 as

where I₁₄ is the local laser intensity (at $n c / 4$⁠) in units of 10¹⁴ W/cm², $λ μ m$ is the (vacuum) laser wavelength in microns, $L μ m$ is L_q in microns, $a o = 0.00199$⁠, and b_SRS= 221. The electron-ion collision frequency, ν_e, is calculated from the LMD model as a function of n_e, T_e, and $Z ¯$ . The effect of Landau damping should be added to the collisional damping term²⁵ in Eq. (2). But collisional damping of the plasmon is found to dominate over Landau damping in the parameter range considered here, so it is not included. For TPD, the corresponding growth rate is given by^15,26

where T_kev is the electron temperature in keV, and $b T P D = 13.1$⁠. We note here that Eqs. (2) and (3) describe only the growth rate of the linear phase of the instability. There are many kinetic and nonlinear effects which affect the evolution of the instabilities, including electron trapping,²⁷ and rescattering.²⁸ In this paper, we are primarily interested in determining the stability of the hohlraum wall, rather than the details of the nonlinear evolution of LPI, and for this purpose we may use the thresholds calculated from linear theory.

Threshold intensities for both SRS and TPD, $I T / S R S$ and $I T / T P D$⁠, may be obtained by setting $γ ′$ to zero in Eqs. (2) and (3), respectively, and solving for intensity. In the limit of a collisionless plasma (⁠ $ν e → 0$⁠), we obtain values for the inhomogeneous thresholds, $I T I / S R S$ and $I T I / T P D$⁠. In the limit of a homogeneous plasma (⁠ $L q → ∞$⁠), both SRS and TPD have the same collisional threshold, $I T C$⁠. The SRS and TPD thresholds, $I T / S R S$ and $I T / T P D$⁠, as calculated by Eqs. (2) and (3) are shown in plot (a) of Fig. 3 for a 1/3–μm laser in a 500-eV gold plasma at the quarter-critical density. The full SRS and TPD thresholds are shown as solid red and blue lines, respectively. The corresponding inhomogeneous thresholds are shown as dashed lines, and the collisional threshold as a dotted black line. The purely inhomogeneous growth rates predict quite low intensity thresholds for relatively large values of L_q, but the stabilizing effect of the electron-ion collisions in the gold results in a full threshold $> 10 15$ W/cm² by more than a factor of two at a temperature of 500 eV. We see also that at smaller values of L_Q that $I T / S R S ≫ I T / T P D$⁠. At $n e ≃ n c / 4$ in a gold plasma the LMD model predicts a collision frequency which drops monotonically with temperature, which results in a lower collision threshold. This effect can be seen in plot (b) of Fig. 3, which shows $I T / T P D$ at 500, 1000, and 1500 eV. The TPD threshold drops below 10¹⁵ W/cm² only for relatively large values of L_q and temperatures approaching a keV. From Fig. 2, we can see that for the hohlraum simulation, L_q does not exceed 50 μm, and T_e reaches a peak temperature of less than 700 eV., at a time near the end of the temporal ramp of 0.1 ns. Thus, Eqs. (2) and (3) predict stability against both SRS and TPD at an incident laser intensity of 10¹⁵ W/cm².

This can also be seen clearly in Fig. 4 where we have plotted the local laser intensity at x_q as a function of time. In the plot, we show the incident laser intensity, I_inc, as a black line, and the local laser intensity, $I ( x q )$⁠, as black dots. The local laser intensity can be calculated from the ray photon number density, n_ph, tabulated on the grid, by the formula $I = n p h ℏ ω c$⁠. Since photons at x_q have either propagated from the vacuum to x_q, or have continued to x_c and been reflected back to the quarter-critical surface, we find that $I ( x q ) ≤ I i n c$⁠, as these photons have deposited much of their incident energy into the plasma via collisional absorption. For these simulation parameters, we find that the photons lose almost all of their energy prior to reaching x_q, with very little energy propagating all the way to x_c and then being reflected back. So the reflected photon intensity at x_q is very small. The fact that $I ( x q )$ decreases with time after the ramp on the incident pulse is due to the fact that as the plasma expands outward and heats, the photons deposit more of their energy prior to reaching the quarter-critical surface. Explicitly $I − 1 d I / d t ∼ ν e n e / n c$ [see Eq. (B13) in  Appendix B]. Additionally, we have plotted the local intensity threshold values for the SRS and TPD instabilities. The full SRS and TPD thresholds are shown as blue and red circles, respectively. The SRS and TPD inhomogeneous thresholds are shown as blue and red diamonds, and the collisional threshold is shown as a green diamond. Each black circle on the plot represents an intensity value calculated at the instantaneous quarter-critical surface. According to the threshold estimates, only times for which the local intensity is less than the threshold values will be stable against TPD and/or SRS. That is, stability requires the black circles to fall below both the red and blue circles. So the hohlraum simulation at 10¹⁵ W/cm² is stable at all times against both SRS and TPD. As the plasma wall evolves in time, the inhomogeneous thresholds drop with increasing L_q, but the effect of the collisions stabilizes the hohlraum even at large values of L_q. As stated above, the plasma temperature does not increase enough for the collision threshold to be exceeded at the quarter-critical surface.

We next consider a series of 1D hohlraum simulations with varying incident laser intensity. In Fig. 5, we have shown the gold ion density profiles for simulations with intensities of 10¹⁴ (black), $3 × 10 14$ (red), and 10¹⁵ W/cm² (blue) at t = 2 ns. As would be expected, the plasma expansion proceeds more rapidly at higher intensities. Other quantities such as x_q, L_q, $T e ( x q ) , Z ¯ ( x q )$⁠, and the vacuum radiation temperature also exhibit reasonable scaling behavior with intensity, with each of these quantities increasing with laser intensity. The results found here cannot be quantitatively compared with the scaling laws compiled in the theory Ref. 22, as this treatment assumes a Spitzer collision frequency, $ν e ∝ T e − 3 / 2$⁠, which is not valid for gold plasma with $Z ¯ ≫ 1$ for $n e ∼ n c$⁠.

The quarter critical laser intensities, $I ( x q )$⁠, and threshold intensities, $I T / T P D , I T I / T P D$⁠, and $I T / C$⁠, are shown as a function of time in Fig. 6 for the three different intensity values. The quarter-critical intensity values are shown as solid circles, and the full threshold values as open circles. The purely inhomogeneous thresholds are shown as open diamonds. For all three intensity values, the local laser intensity is well below the full TPD threshold (⁠ $≤$ SRS threshold), suggesting stability for incident intensities $≤ 10 15$ W/cm². It is interesting to note that for all three values of incident intensity, the resulting inhomogeneous threshold intensity at x_q is roughly the same. Since the temperature and density lengthscale both scale similarly with I_inc, and $I T I / T P D$ is a function of $T e / L q$ [Eq. (3)], the inhomogeneous threshold value appears to be independent of I_inc.

Since the NIF hohlraum generally contains a background helium gas,¹ we now perform 1D hohlraum simulations in the presence of a background gas. These simulations are performed in the same manner as the vacuum-gold simulations. We again utilize the MHD algorithm, but now make use of the multi-ion capability recently added to Chicago, and described in  Appendix A. Proper behavior of the gold-helium interface in the simulation requires the extra terms in the ion momentum equation which are introduced in the multi-ion formulation. The initial simulation geometry is the same as in the top left plot of Fig. 1, with the inclusion of the uniform He background gas in the region x > 0.1 mm. The 1/3–μm wavelength ray-photons, with $I i n c = 10 15$ W/cm², are injected from the right simulation boundary and propagate through 2 mm of He gas before striking the gold. In Fig. 7, number density profiles for gold, helium, and the single electron species are shown for simulations with a He mass density of (top row) $ρ H e = 0$⁠, (middle) 0.03 mg/cm³, and (bottom) 0.96 mg/cm³. Densities are shown t = 0.2, 0.5, and 1 ns (left to right). The full 2 mm extent of the He gas is not shown in the figure. The introduction of the background He gas is seen to have some retarding effect on the hydrodynamic expansion of hohlraum wall compared to the vacuum case. For the low density case (⁠ $ρ H e = 0.03$ mg/cm³), there is a forward propagating shock wave seen at the location of Au-He interface, with a shock front seen in both gold and helium. Note however that the shock occurs in the gold at a density (⁠ $∼ 10 20$ cm⁻³) which is far below the quarter-critical value. For densities $≳ 10 21$ cm⁻³, the gold and electron density profiles for the vacuum and low density fill cases are virtually the same.

At the high helium density value, $ρ H e = 0.96$ mg/cm³, shocks emerging from the material interface propagate both forward and backward, with the result that there are perturbations in the gold density which occur even in the region around the quarter-critical density. This can be seen more explicitly in Fig. 8 which shows a closeup of the quarter-critical region for the high density helium fill run at $t = 0.57$ ns and 0.61 ns. We see that at the earlier time there is no longer one, but three quarter-critical surfaces (top plot), each with a corresponding gradient length-scale, calculated from the local density slope (the slope of the thick dashed lines in the plot). The temperature at each of the three critical surfaces can vary significantly as well (bottom plot). At the slightly later time, the electron density profile has evolved in such a way that there is again only one quarter-critical surface. For both the low density and high density gas fill cases, we have again plotted x_q, L_q, $T e ( x q )$⁠, and $Z ¯ ( x q )$ in Fig. 9. At the low helium density, the data are nearly indistinguishable from the vacuum results shown in Fig. 2, so the vacuum data are not replotted here. For the high helium density simulation, we have plotted data for each quarter-critical surface at each time value. There is always either one value of x_q or three at each time. Other multiplicities were not observed. Note that there is an interval from $∼ 0.2$ to 0.6 ns, where there is an elevated value of $T e ( x q )$ in the high density case. The corresponding L_q values in this time window can range widely from a few tenths to several hundred microns.

The resulting values for $I ( x q )$ and $I T / T P D$ are plotted for the high helium density case in Fig. 10. The corresponding data for the low helium density case are so close to the vacuum data shown in Fig. 4 that they are not shown. Since there can be either one or three quarter-critical surfaces (see Fig. 8) in the high helium density case, we have plotted the resulting data in three separate subplots to aid in visualizing in the data. At each time, the data from the multiple quarter-critical surfaces are plotted from left to right in the descending order of x_q, that is, the left-most plot shows the data corresponding to the largest value of x_q. If only one value of x_q is present, there is no corresponding data point in the middle and right plots.

We see in plot (a) that there is a time interval from about 0.4–0.6 ns during which the intensity approaches the TPD threshold, and actually becomes unstable at $t ∼ 0.43$ ns. In this time interval [see Figs. 9(b) and 9(c)], there exist quarter-critical surfaces which have both relatively large values of L_q (several tens of microns) and temperatures over 1 keV. A gold plasma with $I ∼ 10 15$ W/cm² and $λ = 1 / 3 μ$m was shown to be TPD instable in this regime in plot (b) of Fig. 3.

Although the presence of the He gas might be expected to increase the stability of the hohlraum, since the growth of L_q should be retarded by the fill. The hybrid-PIC simulations suggest that there are circumstances in which the plasma can become more susceptible to TPD due to the detailed time evolution of the quarter-critical surface(s) of the hohlraum wall. The decrease in threshold with L_q accounts for the fact that there is only a finite spatial region for which $n e ∼ n c / 4$ in the presence of a density gradient. A glance at the black curve on the top of Fig. 8 shows that it is possible to have a local maximum in n_e near the critical surface, which case L_q becomes very large. This does not necessarily imply an extremely long interaction region. In such a case, the density profile in the interaction region can be more carefully treated.²⁹ We note however that for the unstable region seen in Fig. 10 for $t ∼ 0.43$ ns, the quarter-critical density surface is nowhere near a density maximum, and a linear approximation for n_e near x_q is appropriate.

In the above analysis, we have used the instability threshold estimate of Eqs. (2) and (3) in concert with 1D vacuum hohlraum simulations, and have determined that the hohlraum walls should be SRS and TPD stable for $I ≤ 10 15$ W/cm² for pulse lengths around 1 ns. Including a He gas (⁠ $ρ H e ∼ 1$ mg/cm³) into the hohlraum introduces shock structures into the gold blowoff plasma, which can somewhat unpredictably affect the TPD stability. Notably, we may have several LPI resonant regions due to the blowoff plasma structure. Depending on the details of the density and temperature profiles, which vary in time, the plasma may be either more or less sensitive to TPD, and may become TPD unstable at quarter-critical surfaces with large L_q and large temperature. In Secs. III and IV, we present the results of fully kinetic LPI simulations to assess the reliability of the threshold estimates used above, and observe the generation and interplay of the SRS and TPD instabilities in parameter regimes suggested by the hybrid hohlraum simulation results.

In this section, we describe the results of 1D fully kinetic LPI simulations. Since these simulations are 1D, there is no possibility of observing TPD (the TPD growth rate vanishes when the plasmon wave vectors are parallel to the laser wave vector⁵) However, at a sufficiently high intensity, SRS growth is clearly visible. As in the hohlraum simulations, the laser wavelength is 1/3 μm (⁠ $n c = 10 22$ cm⁻³). Initial density profiles are shown in Fig. 11. The figure shows how the L_q value can be varied in parametric studies which follow later in the section, with x_q always located at x = 0. For now, however, we consider the case for L_q = 50 μm. The maximum plasma density is chosen to be $< n c$ so that the EM energy may propagate through the plasma without reflection (in the absence of SRS). This choice was made to mimic the situation in the hohlraum simulations in which most of the laser energy is absorbed prior to reaching the critical surface, and thus, there is no significant reflected laser energy incident on the quarter-critical surface. Initial values of temperature, $T e = T i = 500$ eV, $L q = 50$ μm, and $Z ¯ = 50$ were chosen to be in rough agreement with the data in Fig. 2 for an incident intensity value of 10¹⁵ W/cm². An electromagnetic wave is again launched from the right simulation boundary.

The simulations are performed using fully explicit kinetic particles for both electrons and ions. The ray photons of the hohlraum simulations are replaced by actual injected EM fields. It is not possible to retain the EOS modeling in kinetic simulations which are run at short timescales. Radiation can also not be easily included. Such effects in a kinetic PIC framework would require some kind of collisional-radiative modeling capability which is available only in a very rudimentary form in Chicago at present. For this reason, the ion-charge state in the following simulations is fixed. The timestep is chosen so that $τ / Δ t ∼ 50$ where $τ = 2 π / ω o$⁠. The incident laser angular frequency is denoted as ω_o in accordance with the nomenclature introduced above and in  Appendix C, as well as to distinguish it from the general angular frequency variable used in fast Fourier transforms (FFTs) of simulation data presented later in this section. The cell size is $Δ x ≃ λ / 30$⁠, and the plasma is represented by $∼ 500$ electron macroparticles per cell, with fewer gold ion particles required. Electron-ion collisions are again calculated by the LMD model.

These highly resolved simulations are run for about ten thousand laser periods, which is much less than the ns length pulse durations considered in the hohlraum simulations of Sec. II. On such short timescales, it is not possible to keep the electron and ion temperatures equilibrated. As a result, in the LPI simulations, the electrons are still heated by collisional absorption and other mechanisms, but the ions essentially remain at their initial temperature. Although we have carefully chosen the initial conditions in the highly resolved LPI simulations to agree with the conditions found in the hohlraum simulations, in this sense also, the two simulation types are not entirely self-consistent, which must be considered when comparing results from the two sets of simulations.

Some initial results are shown in Fig. 12 for which the incident intensity is chosen to be 10¹⁵ W/cm². In the short timescale LPI simulations, the temporal ramp on the laser intensity is chosen to be 3 laser periods. On the left plot of Fig. 12, we have plotted the measured intensity at $x = 9$ μm, and x = −11 μm. We note here that the “intensity” in the context of EM waves means the measured Poynting flux, and thus contains (in principle) contributions from EM energy propagating in both forward and reverse directions, although the incident EM energy is injected only from the right. For illustrative purposes, we show the results using both ray-tracing photons and actual EM fields. For the pure EM simulation (red trace), we have applied a low-pass filter to the data to obtain effectively “cycle-averaged” results. For both methods, we see that for this incident intensity, at x = 9 μm, there is essentially no reflection of laser energy, which would be evidenced by a change in the measured intensity at that position. From the data at x = −11 μm, we see that about 60% of the head of the pulse is transmitted through the plasma, with the remaining 40% of the energy deposited into the plasma by collisional absorption. As the simulation continues, the transmission coefficient gradually increases with time. This is due to the bulk heating of the plasma which can be observed in the right plot of Fig. 12, which shows T_e as a function of x at $ω o t = 0$⁠, 2400, and 4800. As T_e increases from its initial value during the simulation, the collision frequency drops which reduces the absorption of energy by the plasma. At intensities significantly lower than 10¹⁵ W/cm², there is no significant electron heating on such short timescales and the transmission coefficient remains essentially constant in time. The results of the ray photon and EM wave treatments are seen to be in agreement, which is to be expected as long the density gradient scalelengths in the problem are long compared to the laser wavelength, and the intensity is low enough to avoid LPI instabilities such as SRS, which cannot be modeled by the ray-tracing approach.

The SRS instability, visibly absent at 10¹⁵ W/cm², can be clearly seen by increasing the laser intensity by a factor of ten. Some results at this higher intensity are shown in Fig. 13. The measured intensity at x = 9 μm begins to drop noticeably for $ω o t ∼ 2000$⁠, which indicates the presence of reflected wave energy. The precipitous drop in transmission at x = −11 μm at about the same time also indicates enhanced absorption in the plasma. An FFT of the transverse field, E_z, at the same two positions shows the presence of SRS lines. At x = 9 μm, the strong line $ω ∼ ω o / 2$ is consistent with a directly back-scattered wave, where $ω s ∼ ω o − ω p e ∼ ω p e$⁠, where ω_pe is calculated at $n c / 4$⁠, and ω_s is the angular frequency of the scattered EM wave.

The presence of higher harmonics at $ω / ω o ∼ 1.5$⁠, 2, and 2.5, which are not present if the FFT is taken over a time range restricted to the linear phase of SRS growth (⁠ $ω o t < 2000$⁠), indicates that nonlinear effects are also significant. The FFT of E_z at x = −11 μm also shows a prominent forward-scattered SRS peak at $ω o + ω p e$⁠. But this line is also only present after the linear phase of the instability has concluded. For values of $L q < 100$ μm, there appears to be no forward-scattered SRS peak in the linear regime. At significantly larger values of L_q, there are both strong forward and backward SRS peaks at $ω / ω o ∼ 0.5$ in the linear regime, which is consistent with the theory of Ref. 25.

For the same simulation, snapshots of electron density and electric field components, E_z and E_x, are shown in Fig. 14 in a small spatial region centered around x_q = 0. Since the incident EM wave vector is in the −x-direction and is polarized in the virtual z component, any electric field component in the x-direction is indicative of a plasmon, for which the electric field and wave number are parallel. The figure shows results at $ω o t = 1200$⁠, 1600, and 2000. At the earliest time, the first trace of the plasmon field component is visible. At the intermediate time, a perturbation in the electron density can be seen to be emerge above the particle noise. The beating of the incident and scattered EM waves can also now clearly be seen. At the latest time, the plasmon signal has become comparable in amplitude to the incident field. At later times, the behavior becomes increasingly non-linear in character.

The growth rate of the 1D LPI simulations is calculated by integrating the square of the instantaneous plasmon field over the interaction region and fitting the results to obtain $2 γ ′$ during the linear phase of the instability growth. We have performed a parametric study in which both the laser intensity and L_q are varied. The density profiles are modified to achieve various values of L_q as illustrated in Fig. 11, while the intensity values range from 10¹⁵ to 10¹⁶ W/cm². All of the simulations retain the 1/3 μm wavelength, use a fixed $Z ¯$ of 50, and use initial temperatures of $T e = T i = 500$ eV. The results of this study are presented in Sec. IV together with results from 2D simulations.

In this section, we extend the LPI simulations of Sec. III to two dimensions. A contour plot of the initial gold density profile for a 2D Cartesian LPI simulation is shown in Fig. 15. The axial density profile is the same as shown in Fig. 11 for L_q = 50 μm. The initial density profile is uniform in the transverse (z) direction. The simulation extends over 6 μm in the z-direction with $Δ z ∼ Δ x$⁠. The laser propagates in the −x-direction and is polarized in the z-direction. The laser intensity is also uniform in the z-direction, which is a reasonable approximation to make for beam with a waist much larger than a few μms. All other physical and numerical parameters are unchanged from the 1D simulations.

The presence of a plasmon signal, due to either SRS, or TPD in 2D, is still inferred by the presence of an E_x component in the field. In Fig. 16, we have plotted the contours E_x in the interaction region near x_q for a 2D simulation with $I = 10 16$ W/cm² and L_q = 50 μm. The results are shown for $ω o t = 400$⁠, 800, and 1200 (from left to right). For the same parameters in 1D, we saw a strong SRS signal. This can be seen in the 2D results as well at $ω o t = 400$⁠, as there is clearly a plasmon signal with $k →$ parallel to $x ̂$⁠. Given the lack of structure in the z direction, and since the growth rate for TPD vanishes in this propagation direction, this should be the analog of the SRS instability exhibited in 1D. By $ω o t = 800$⁠, some 2D structure begins to emerge. This structure grows in amplitude with time and is seen to be localized around the quarter-critical surface at $ω o t = 1200$⁠. Careful inspection of the E_x signal shows field structures consistent with localized plasmons propagating at ±45° to the x-axis, which is consistent with the fact that strong TPD growth is predicted in this case. By $ω o t = 1200$⁠, this TPD signal is much larger in amplitude than the SRS (note that the three plots in Fig. 16 all have different ranges on the color bar in order to emphasize the plasmon signal at each time).

In Fig. 17, we have shown the electron density, and electric field components at $ω o t = 2000$ for both 1D (bottom row) and 2D (top) simulations. In two dimensions, we can see both an SRS signal (structure independent of z) at x > 1 μm which corresponds roughly to the 1D signal and a higher amplitude TPD signal more localized around the quarter-critical surface. From the 2D results, the amplitude E_x is already comparable to the unperturbed amplitude of E_z, suggesting that the linear phase of the instability growth has ended. At later times, the simulation results become increasingly nonlinear in character. The 1D and 2D results for E_z both have a reflected EM wave for x > 1 μm, which is consistent with back-scattered SRS.

The effect of the LPI instabilities, in concert with collisional absorption, on the EEDF, f(E), can be seen in Fig. 18. The EEDFs are constructed at each time from a histogram in energy created from all macroelectrons in the range from −2.5 to 2.5 μm in x, with the result normalized so that

The initial EEDF corresponds to a maxwellian distribution at 500 eV. At later times, EEDFs take on a distinct non-maxwellian character, with a high temperature tail becoming visible for $ω o t$ values greater than a few hundred. The high temperature tail increases rapidly until about $ω o t ∼ 2400$⁠, after which time the EEDF seems to converge, and the saturation of the instabilities is visible in the time-dependence of the plasmon signal.

We now consider a parametric study of both 1D and 2D LPI simulations in which both the laser intensity and L_q are varied. The growth rates for the instabilities during the linear phase can be calculated from the simulation data as described in Sec. III. The growth rates calculated from the PIC simulations are then compared to the theoretical maximum growth rate for SRS (in 1D) and TPD (in 2D) given by Eqs. (2) and (3). These equations are evaluated by using the simulation values of T_e, I, and ν_e at the quarter-critical surface (x_q = 0 for all LPI simulations). When electron collisions are included in the LPI simulations, however, both the electron temperature and the intensity at x_q vary as a function of time, as was shown in Fig. 12. Note that when collisions are artificially neglected, both I and T_e at x_q are roughly constant before the onset of the instabilities. In Fig. 19, we show the time variation of the quarter-critical electron temperature and plasmon field component (top), as well as the upstream and quarter-critical intensity for a 1D simulation with L_q = 200 μm Which intensity and temperature should be used in the growth rate formulas? It is seen that the time of peak quarter-critical intensity, t_max, corresponds roughly to the emergence of the plasmon signal above the noise floor, this allows us to define an effective intensity and temperature, $I ( x q , t m a x )$ and $T e ( x q , t m a x )$ which can be used in Eqs. (2) and (3). In the entire parameter study, the effective intensity never falls below $0.9 × I i n c$ and the effective temperature never exceeds 1 keV. The growth rates for both 1D and 2D simulations are shown in Fig. 20. The growth rates calculated directly from the simulations are shown in the figure as closed circles, while the open squares represent the growth rates from Eqs. (2) and (3), in which we have used $I ( x q , t m a x )$ and $0.8 × T e ( x q , t m a x )$ for the intensity and electron temperature, respectively. The factor of 0.8 in the temperature calculation was chosen to give the best agreement with the theory and the simulations in the homogeneous limit. It is plausible to assume that the linear phase of the instability growth begins somewhat prior to t_max [see Fig. 19], which justifies assuming a temperature value somewhat less than $T e ( x q , t m a x )$⁠, since the local intensity has a fairly broad maximum at t_max we use the full value of $I ( x q , t m a x )$⁠.

In the 1D simulations we have, of course, compared the simulation results to the theoretical SRS growth rate. In 2D, TPD and SRS can coexist but for finite values of L_q, and temperatures on the order of 1 keV, TPD growth should dominate SRS. For that reason, we have plotted the TPD growth rate from Eq. (3) as squares in the lower plot of Fig. 20. In 1D, there is good agreement between theory and simulation for large values of L_q, while there is somewhat poorer agreement as L_q decreases. One explanation for this is that Eq. (2), as derived in Ref. 24, is strictly valid in the regime where the second term in the parentheses is much smaller than unity, so the equation is expected to break down as L_q decreases. The effect was also noted in the model for SRS described in Ref. 15, in which numerically calculated growth rates exceeded Eq. (2) as L_q was decreased. Another issue worthy of mention is that, as L_q decreases and the growth rates decreases, it becomes increasingly more difficult to measure the growth rate, as the simulations must be run for longer times and the signal-to-noise ratio degrades. So the effective error bars on the simulation data become greater. In 2D, the simulation growth rates are in relatively good agreement with Eq. (3). In both 1D and 2D we found no noticeable instability growth at any value of L_q for incident intensities at or below 10¹⁵ W/cm² for runtimes out to $10 4 ω o − 1$⁠. This is also consistent with the theory results (see Fig. 3) for temperatures below a keV.

It is possible, by methods described in Refs. 15 and 29, to calculate more accurate growth rates than the approximate expressions given by Eqs. (2) and (3). However, for this purpose of this paper, these expressions, when compared to PIC results, have been shown to be accurate enough justify their use in Sec. II to estimate SRS and TPD thresholds in the hybrid simulations of hohlraum wall expansion.

We have performed hybrid PIC simulations of gold hohlraum walls based on NIF parameters and investigated the possibility of hot electron production by TPD and SRS. The simulations of gas-filled hohlraums on hydrodynamic timescales required the development of a multi-ion PIC-based MHD algorithm in order to properly model the details at the interface of the fill gas and the blow-off gold plasma. The collisional absorption of laser energy was modeled by a PIC-based ray-tracing algorithm. The hohlraum simulations using the algorithm discussed here were 1D, but the algorithms are fully generalized for multiple-dimensions and can be used to model multiple laser beams of arbitrary cross-section and pulse shape. The relatively simple ray-tracing treatment used here to model collisional absorption can be extended to allow for more complicated phenomena such as magnetic fields and coupling between plasma wave modes.³⁰ But it is ultimately limited to linear phenomena.

From the results of the hohlraum simulations, theoretical intensity thresholds for TPD and SRS suggest that the gold walls in a vacuum hohlraum should be stable for $I i n c ≤ 10 15$ W/cm². The presence of a He gas with mass density on the order of 1 mg/cm³ may however make the walls less stable to TPD, since shock waves emerging from the Au-He interface can distort the plasma density and temperature profile near the quarter-critical surface in such a way that the TPD threshold intensity is exceeded.

Results from the 1D hohlraum simulations were used to initialize highly resolved LPI simulations performed on much shorter timescales, in which both the SRS and TPD instabilities were observed for $I ≳ 10 15$ W/cm² for values of L_q consistent with the hohlraum simulations. A non-thermal hot electron tail was also observed in the EEDF. In 1D and 2D parametric studies, we have found the simulation growth rates in good general agreement with the approximate theory values given by Eqs. (2) and (3) for $L q >$ a few tens of microns. The LPI simulations were performed on timescales shorter than the electron-ion equilibration time, which leads to unequal electron and ion temperatures. Moreover, the LPI simulations were performed with a fixed charge state and an ideal equation of state model. In this sense, the “handoff” between the hybrid hohlraum simulations and the kinetic LPI simulations was not completely self-consistent. Taking these caveats into account, we nonetheless believe that the general conclusions stated above should be valid.

The authors would like to thank Dr. A. S. Richardson of NRL for a useful discussion regarding some aspects of the work. This work was supported by Sandia National Laboratory. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the National Nuclear Security Administration under DE-AC4-94AL85000.

In a previous publication,⁹ we have outlined a PIC-based MHD algorithm which was valid for a single ion species. In order to allow for hohlraum simulations in the presence of background gases, we have extended this algorithm to allow for multiple MHD ion species. In this section, we describe the formulation of the multi-ion MHD equations, and then briefly describe how they are implemented in the hybrid PIC code Chicago.

In the following, we briefly outline the procedure used to derive the multi-ion MHD momentum equations. The derivation given here is similar to the one presented by Glocer et al.³¹ We allow for an arbitrary number, N_i, of MHD ion species and a single electron species, with quasi-neutrality assumed

where n_e and n_k are the electron and ion densities, respectively, and $Z ¯ k$ is the charge state of the k-th ion. The charge density can then be written as

where

and $ρ k = ( Z ¯ n k / n e )$⁠.

The initial non-MHD fluid momentum equations are given for the k-th ion species

and the single electron species

where $v →$⁠, T, and P, with the appropriate subscripts denoting the species, are the fluid velocity, temperature, and pressure, respectively. The general binary collision frequency, $ν α β$⁠, for particles of species α and β can be given by the Spitzer value (Chicago uses a general binary form given in Ref. 32). Alternately, the electron-ion collision frequency can be calculated by the Lee-More model with Desjarlais corrections (LMD).^19,20 Conservation of momentum implies that $n α m α ν α β = n β m β ν β α$⁠. Equations (A4) and (A5) are an extension of the usual plasma two-fluid model to allow for multiple ion species.³² We have also allowed for an arbitrary equation-of-state (EOS) model from which $Z ¯$⁠,P, and the internal energy may be obtained as a function of the local densities and temperatures.¹⁶ We have also included electro-thermal force term $∼ − β ∇ T e$ where β is a dimensionless transport coefficient.³³ In a magnetized plasma, the transport coefficients (ν and β in this case) are tensor quantities with different values in directions parallel and perpendicular to the magnetic field. This can be included in the formulation but we have neglected it here, since we are considering primarily unmagnetized plasmas.

We have not explicitly included viscosity in the momentum equations but it may be included above by simply adding a viscosity term³³ to the right hand sides of Eqs. (A4) and (A5). It is sometimes useful to include a small amount of artificial viscosity to simulations to avoid numerical problems at shock fronts. In our PIC implementation, this is roughly analogous to what is done in other fluid codes.

Other than the relatively straight-forward extension to allow for multiple ions and the inclusion of the electro-thermal and viscous forces, the fluid energy equations are essentially unchanged from those presented in Ref. 16, which include flux-limited thermal conductivity for both ions and electrons and a radiation term in the electron energy equation to account for absorption and emission of radiation.

The derivation of the MHD equations from the multi-fluid equations is similar to the standard single-ion formulation.³⁴ A generalized Ohm's law is obtained by first multiplying Eq. (A4) by $e Z ¯ k / m k$ and summing over ion species, k. Equation (A5) is then multiplied by $− e / m e$⁠. The resulting two equations are then added together to obtain

where

is the plasma conductivity and

In deriving Eq. (A6), we have dropped terms of order $( m e / m k )$ and made use of the low-frequency approximation $ω ≪ ω p e$⁠, where $ω − 1$ is a characteristic time scale and ω_pe is the electron plasma frequency. We have also neglected the Hall term³⁵ which is of order of $( ω c e / ν e )$⁠, where ω_ce is the electron cyclotron frequency. This is appropriate for unmagnetized and/or weakly magnetized plasmas, and is consistent with the approximation of using scalar transport coefficients.

The MHD ion momentum equations are now obtained by using Eq. (A6) to replace the electric field in Eq. (A4), with the result

The second line on the right hand side of Eq. (A9) contains the effect of ion-ion collisions, while the third and fourth lines contain multi-ion MHD correction terms which include the effects of ion-ion interactions, including collisional terms mediated by the electron species (as evidenced by the presence of ν_ek). These terms are required for a proper MHD treatment of multiple-ion systems, such as the gas-filled hohlraum simulations discussed in Sec. II. Note that in the limit of a single ion species (N_i = 1 and $ρ k = 1$⁠) only the first line on the right hand side of Eq. (A9) is non-zero and the ion equation of motion reduces to the standard MHD result.⁹ 

The multi-ion MHD formulation described above has been implemented into the hybrid-PIC code Chicago and is carried out similarly to the single-ion implementation described in Ref. 9. Fluid ion macroparticles are advanced in time with Eq. (A9) in conjunction with Eq. (A6). Electric and magnetic fields are advanced using the full Maxwell equations with the current density in Ampère's law being supplied by Eq. (A6). Fields, densities, conductivity, and mechanical forces (pressure gradients, etc.) are gathered on the grid nodes and are interpolated to the ion position during the particle advance. The electrons comprise a virtual species with attributes (⁠ $v → e$⁠, T_e) carried by ion particles. The virtual electron weight is determined by the appropriate charge state of the particular ion it is attached to. Both electron and ion energy equations are advanced on the grid and used to update macroparticle temperatures.

The electron energy equation includes a radiation source/sink term which accounts for both emission and absorption of radiation. And the radiation energy density (⁠ $∝ T rad 4$⁠, where $T rad$ is the radiation temperature) is calculated by the multi-group diffusion approximation,²¹ with a flux limiter included in the diffusion operator to approximately handle the transition from optically thick regions (i.e., hohlraum walls) to vacuum.

The EOS data needed for both the plasma and the opacity data needed for the radiation field solve in the diffusion approximation are pre-calculated by the Propaceos code.^17,18 The MHD algorithm may also employ the Eulerian fluid algorithm described in Ref. 16. In this case, following the Lagrangian particle advance, ion MHD macroparticles are remapped back onto the Eulerian grid.

In this appendix, we describe the implementation of a ray-tracing algorithm into Chicago. A very general treatment of ray-tracing methods in plasma physics has been given by Tracy et al.³⁰ In the following discussion, we follow this approach to derive the eikonal equations for an EM wave propagating in an underdense unmagnetized plasma. The same equations have previously been implemented in fluid codes by Kaiser.⁷ A PIC implementation is briefly described below.

In order to introduce some nomenclature used in the following, we present here a brief review of the theory of EM wave propagation in under-dense unmagnetized plasma (see, for example, Ref. 34). The wave equation for an EM wave with an amplitude, E, propagating in a uniform unmagnetized cold plasma is given by

where

The usual dispersion relation is obtained by setting D to zero. Surfaces of constant phase propagate at

where

where

and

Wave packets propagate at the group velocity which is given by

The effect of collision absorption⁵ can be included by adding a sink to the equation for the EM energy density.

where

In the treatment of Tracy et al.³⁰ the eikonal equations of motion can be written in the form of Hamilton's equation in a 6D $x →$- $k →$ phase space.

with a dispersion Hamiltonian given by

Using the dispersion function from Eq. (B2) with a spatially dependent electron plasma frequency, $ω p e ( x → ) ∝ n e 1 / 2$⁠, we obtain the ray equations⁷ 

Integrating these equations in time sweeps out the path of the ray trajectory in space. When collisions are included the Lagrangian ray energy density is given by

We note that this formulation is based on the eikonal (geometric optics) approximation: $λ ≪$ all other relevant lengthscales. Moreover, the absorption model assumes that $ν e / ω ≪ 1$⁠.

An arbitrary wave packet (i.e., laser pulse) is broken into discrete rays which are treated as PIC photon macroparticles which have the following attributes.

• a particle weight which is proportional to the number, N, of real photons represented by the macro-particle.

• Instantaneous photon position.

• Instantaneous photon group velocity.

• photon energy (⁠ $ℏ ω$⁠).

The photon position and group velocity are advanced by Eq. (B12), and the time evolution of the photon weight is given by Eq. (B13) with $u R = N ℏ ω$⁠. The energy lost by the macrophoton is deposited locally and interpolated onto the grid with the lost photon energy added to the internal energy of the local electron population.

The implementation given here is appropriate for unmagnetized plasmas, but the approach can be generalized by using a different dispersion hamiltonian.³⁰ The ray-tracing model described here behaves as if laser light is s-polarized (⁠ $E → · ∇ n e = 0$⁠), which is adequate for describing collisional absorption, but for p-polarized light (⁠ $E → · ∇ n e ≠ 0$⁠), the component of the electric field parallel to $∇ n e$ can couple to a Langmuir wave. This is a separate mechanism of depositing laser energy into the plasma known as resonance absorption⁵ which cannot be captured by the present formulation. Of course, this effect and other parametric instabilities are all captured fully by kinetic PIC simulations in which the laser wavelength and frequency are resolved by the grid cell size and timestep.

In this appendix, we briefly review the parametric instabilities of stimulated Raman scattering (SRS) and two plasmon decay (TPD). There are, of course, many other possible instabilities observed in laser plasma interactions, but these two particular processes are of interest here because they are found in the small-scale highly resolved kinetic simulations described in Secs. III and IV. We note that both SRS and TPD are fluid instabilities, but electrons which are accelerated in the daughter waves can generate non-maxwellian features in the EEDF. There can also exist purely kinetic instabilities, which can be captured only in fully kinetic simulations. Much of the material below has been previously summarized by Kruer.⁵ 

In stimulated Raman scattering, an EM wave (⁠ $ω o , k → o$⁠) is scattered by a electron plasma wave (plasmon), with frequency and wave number matching conditions given by⁵ 

where ω_s and $k → s$ are the angular frequency and wave number of the scattered EM wave in the plasma, and

and $v e = T e / m e$⁠. At low enough temperatures,

In this case, a resonant interaction can occur with $ω o ≳ 2 ω p e$ or when $n e ≲ n c / 4$⁠. For directly backscattered radiation at $n ∼ n c / 4$ in which $ω o ∼ 2 ω p e$ and $ω s ∼ ω p e$⁠, the maximum SRS growth rate is given by^5,25

where $v o s = ( e E / m e ω o )$ is the local electron quiver velocity and E is the peak electric field amplitude. The quiver velocity can, of course, be related to the intensity, $I = 1 / 2 ϵ ε o c E 2$⁠, where $ϵ = 1 − n e / n c$ which is $∼ 3 / 4$ for $n e ∼ n c / 4$⁠.

Equation (C4) is strictly valid only in a homogeneous and collisionless plasma. The presence of both electron-ion collisions and density inhomogeneities decreases the actual growth and imposes a threshold intensity on the instability.

In a collisionless inhomogeneous plasma characterized by an inhomogeneity lengthscale,

the maximum effective SRS growth is modified to read^15,24

where

Instability only occurs when N > 3∕8, or^15,24,36

The collisional threshold in an inhomogeneous plasma is obtained by requiring that the instability growth rate meet or exceed the plasmon (energy) damping rate^5,25

In Eq. (C9), we have assumed that plasmon damping is dominated by electron-ion collisions rather than Landau damping or other mechanisms.

Combining Eqs. (C4) and (C9), we may write the effective SRS growth rate in a collisionless plasma as²⁴ 

This equation is reprinted in practical units in Eq. (2) of Sec. II.

An EM wave (⁠ $ω o , k → o$⁠) can decay into two electron plasma waves. The process requires the frequency and wave nunmber matching conditions

Again, at low temperatures,

which implies that there is a resonant interaction with $ω 0 ∼ 2 ω p e$ or when $n e ∼ n c / 4$⁠. The maximum growth rate occurs when both plasmon wave numbers are at an angle of 45° to both $k o →$ and $v → o$⁠, in which case⁵ 

where γ is again given by Eq. (C4). Equation (C9) also holds for TPD in a collisional plasma.³⁷ So the growth rates for SRS and TPD in a homogeneous plasma should be the same. In the absence of collisions, the TPD inhomogeneous growth rate is modified to read^15,26

where

where α is a very weak function of the parameter $β = 36 ϵ v e 4 / v o s 2 c 2$⁠, and has a value on the order of one over many orders of variation in the magnitude in β.²⁶ 

Equation (C14) leads to an inhomogeneous threshold of

In Refs. 37 and 5,α = 1, while Ref. 26 has $α = 0.92$ for $β ≪ 1$⁠, and 1.38 for $β ≫ 1$⁠. In the presence of collisions, the TPD growth rate becomes

Choosing the small-β value for α (as in Ref. 15), this equation has also been written in practical units in Eq. (3) in Sec. II.