The walls of the hohlraum used in experiments at the national ignition facility are heated by laser beams with intensities $ \u223c 10 15 $ W/cm^{2}, a wavelength of $ \u223c 1 / 3 $ μm, and pulse lengths on the order of a ns, with collisional absorption believed to be the primary heating mechanism. Xrays generated by the hot ablated plasma at the gold walls are then used to implode a target in the hohlraum interior. In addition to the collisional absorption of laser energy at the walls, nonlinear laserplasma interactions (LPI), such as stimulated Raman scattering and two plasmon decay, are believed to generate a population of suprathermal electrons which, if present in the hohlraum, can have a deleterious effect on target implosion. We describe results of hohlraum modeling using a hybrid particleincell code. To enable this work, new particlebased algorithms for a multipleion magnetohydrodynamic (MHD) treatment, and a particlebased raytracing model were developed. The use of such hybrid methods relaxes the requirement to resolve the laser wavelength, and allows for relatively largescale hohlraum simulations with a reasonable number of cells. But the nonlinear effects which are believed to be the cause of hot electron generation can only be captured by fully kinetic simulations with good resolution of the laser wavelength. For this reason, we employ a twotiered approach to hohlraum modeling. Largescale simulations of the collisional absorption process can be conducted using the fast quasineutral MHD algorithm with fluid particle species. From these simulations, we can observe the time evolution of the hohlraum walls and characterize the density and temperature profiles. From these results, we can transition to smallerscale highly resolved simulations using traditional kinetic particleincell methods, from which we can fully model all of the nonlinear laserplasma interactions, as well as assess the details of the electron distribution function. We find that vacuum hohlraums should be stable to both two plasmon decay and stimulated Raman scattering instabilities for intensities $ \u2264 10 15 $ W/cm^{2}. In gasfilled hohlraums, shocks may be induced in the blowoff gold plasma, which leads to more complex density and temperatures profiles. The resulting effect on LPI stability depends strongly on the details of the profile, and it is possible for the gasfilled hohlraum to become unstable to two plasmon decay at 10^{15} W/cm^{2} if the quartercritical surface reaches temperatures exceeding 1 keV.
I. INTRODUCTION
Indirectdrive thermonuclearignition experiments^{1} are presently being conducted on the national ignition facility (NIF).^{2} These experiments use the xray heating caused by the deposition of laser energy into a blowoff plasma on the walls of a hohlraum to drive the implosion of a DTfuel capsule.^{3} While progress has been made in understanding these experiments since NIF became operational, substantial disagreement remains between simulated and measured neutron yields.^{4} The interaction of the NIF lasers with the inner surface of the hohlraum can generate a hot nonthermal electron population, as well as thermal xrays. The hot electrons can be produced from multiple laserplasma interactions (LPI) including parametric instabilities, such as stimulated Raman scattering (SRS) and two plasmon decay (TPD).^{5} Hot electrons generated in this way can penetrate the capsule and deposit unwanted preheat.^{6} Presently, simulations of NIF shots are conducted using radiationhydrodynamic fluid codes. Fluid codes generally model the collisional absorption of laser energy in the hohlraum by a raytracing algorithm^{7} based on the eikonal approximation ( $ \lambda \u226a $ relevant lengthscales). However, parametric instabilities such as SRS and TPD cannot be modeled by this approach, as they require a full electromagnetic (EM) treatment of the fields. Moreover, the nonthermal 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 particleincell (PIC) code to study the problem of parametric instabilities in the blowoff plasma in the hohlraum walls. We have used the code “Chicago” for the work presented.^{8} Chicago is a fully relativistic, fully electromagnetic, massively parallel, particleincell code that is capable of conducting 1D, 2D, and 3D simulations. Chicago can conduct fully kinetic calculations as necessary and fluidparticle 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 blowoff 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 $ \lambda \u223c 1 / 3 $ μm, interact with the highZ (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 twotiered approach to hohlraum modeling. An initially thin blowoff plasma is assumed to be present at the walls. Largescale simulations ( $ \Delta x > \lambda $) of the hohlraum wall can be conducted using a fast quasineutral PICbased magnetohydrodynamic (MHD) algorithm^{9} 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^{10} against the fluid code Hydra^{11} for several test problems with good agreement found between the pure fluid code and the hybridPIC approach. Moreover, the MHD algorithm in Chicago was found to agree well with experimental data in simulations of a dense plasma focus.^{12} The deposition of laser energy into the wall is modeled by tracking macroparticle photons which obey raytracing 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 smallerscale highly resolved ( $ \Delta x \u226a \lambda $) simulations using traditional kinetic PIC methods, from which we can fully model all of the nonlinear laserplasma 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 quartercritical 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.,^{15} 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 $ \u2273 10 15 $ W/cm^{2} may be sufficiently large to introduce the TPD instability in gasfilled hohlraums due to the detailed nature of the time evolution of the quartercritical surface. In vacuum hohlraums, by contrast, the stabilizing effect of collisions should be sufficient to suppress instability for $ I \u2264 10 15 $ W/cm^{2}.
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 gasfilled 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,^{1} the singleion MHD model discussed in Ref. 9 was extended to allow for multiple ion species. Similarly, the raytracing 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.
II. 1D HYBRID GOLD HOHLRAUM SIMULATIONS
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 \u2264 0.1 $ mm) have a density of $ 5.9 \xd7 10 22 $ cm^{−3} 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, $ \lambda = 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^{2} 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^{−3} (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.^{16} 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 precalculated by the Propaceos code.^{17,18} The electronion collision frequency is obtained from the LeeMore model with Desjarlais corrections (LMD),^{19,20} and we use a singlegroup (gray) radiation solve^{16,21} to calculate the radiation transport in the hohlraum. We use a variable cell size with a minimum value of $ \Delta x \u223c 0.5 $ μm in the region of the density ramp. Laser propagation and energy deposition are modeled by a raytracing algorithm^{7} 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 quartercritical 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 plasmavacuum interface. The radiation temperature quickly reaches a steadystate value of $ \u223c 200 $ eV in the optically thin “quasivacuum” region, which we define roughly here as the region, $ x \u2273 x c $.
To quantitatively describe the blowoff 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 \xaf $ 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 laserdriven ablation in plasmas in which absorption is dominated by inverse bremsstrahlung, which predict a selfsimilar 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 quasisteadystate situation, $ L c \u221d t , \u2009 T e ( x c ) $ and $ Z \xaf ( 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.^{5} 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 \u2264 n c / 4 $, and TPD when $ n e \u223c n c / 4 $ (see Appendix C for details). In the LPI simulations discussed in Sec. III, we find a strong backscattered EM signal and a localized plasmon signal near the quartercritical surface, $ n e \u223c 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 backscattered 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 quartercritical density, including the effects of both plasma inhomogeneity (finite L_{q}) and electronion collisions, can be written^{15,24} as
where I_{14} is the local laser intensity (at $ n c / 4 $) in units of 10^{14} W/cm^{2}, $ \lambda \mu m $ is the (vacuum) laser wavelength in microns, $ L \mu m $ is L_{q} in microns, $ a o = 0.00199 $, and b_{SRS}= 221. The electronion collision frequency, ν_{e}, is calculated from the LMD model as a function of n_{e}, T_{e}, and $ Z \xaf $ . The effect of Landau damping should be added to the collisional damping term^{25} 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,^{27} and rescattering.^{28} 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 $ \gamma \u2032 $ to zero in Eqs. (2) and (3), respectively, and solving for intensity. In the limit of a collisionless plasma ( $ \nu e \u2192 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 \u2192 \u221e $), 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 500eV gold plasma at the quartercritical 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 electronion collisions in the gold results in a full threshold $ > 10 15 $ W/cm^{2} 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 \u226b I T / T P D $. At $ n e \u2243 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^{15} W/cm^{2} 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^{15} W/cm^{2}.
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 \u210f \omega 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 quartercritical surface, we find that $ I ( x q ) \u2264 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 quartercritical surface. Explicitly $ I \u2212 1 d I / d t \u223c \nu 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 quartercritical 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^{15} W/cm^{2} 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 quartercritical 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^{14} (black), $ 3 \xd7 10 14 $ (red), and 10^{15} W/cm^{2} (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 ) , \u2009 Z \xaf ( 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, $ \nu e \u221d T e \u2212 3 / 2 $, which is not valid for gold plasma with $ Z \xaf \u226b 1 $ for $ n e \u223c n c $.
The quarter critical laser intensities, $ I ( x q ) $, and threshold intensities, $ I T / T P D , \u2009 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 quartercritical 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 ( $\u2264$ SRS threshold), suggesting stability for incident intensities $ \u2264 10 15 $ W/cm^{2}. 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,^{1} we now perform 1D hohlraum simulations in the presence of a background gas. These simulations are performed in the same manner as the vacuumgold simulations. We again utilize the MHD algorithm, but now make use of the multiion capability recently added to Chicago, and described in Appendix A. Proper behavior of the goldhelium interface in the simulation requires the extra terms in the ion momentum equation which are introduced in the multiion 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 rayphotons, with $ I i n c = 10 15 $ W/cm^{2}, 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) $ \rho H e = 0 $, (middle) 0.03 mg/cm^{3}, and (bottom) 0.96 mg/cm^{3}. 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 ( $ \rho H e = 0.03 $ mg/cm^{3}), there is a forward propagating shock wave seen at the location of AuHe interface, with a shock front seen in both gold and helium. Note however that the shock occurs in the gold at a density ( $ \u223c 10 20 $ cm^{−3}) which is far below the quartercritical value. For densities $ \u2273 10 21 $ cm^{−3}, the gold and electron density profiles for the vacuum and low density fill cases are virtually the same.
At the high helium density value, $ \rho H e = 0.96 $ mg/cm^{3}, 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 quartercritical density. This can be seen more explicitly in Fig. 8 which shows a closeup of the quartercritical 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 quartercritical surfaces (top plot), each with a corresponding gradient lengthscale, 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 quartercritical 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 \xaf ( 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 quartercritical 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 $ \u223c 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 quartercritical 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 quartercritical surfaces are plotted from left to right in the descending order of x_{q}, that is, the leftmost 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 \u223c 0.43 $ ns. In this time interval [see Figs. 9(b) and 9(c)], there exist quartercritical surfaces which have both relatively large values of L_{q} (several tens of microns) and temperatures over 1 keV. A gold plasma with $ I \u223c 10 15 $ W/cm^{2} and $ \lambda = 1 / 3 \mu $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 hybridPIC simulations suggest that there are circumstances in which the plasma can become more susceptible to TPD due to the detailed time evolution of the quartercritical 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 \u223c 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.^{29} We note however that for the unstable region seen in Fig. 10 for $ t \u223c 0.43 $ ns, the quartercritical 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 \u2264 10 15 $ W/cm^{2} for pulse lengths around 1 ns. Including a He gas ( $ \rho H e \u223c 1 $ mg/cm^{3}) 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 quartercritical 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.
III. 1D LPI SIMULATIONS
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^{5}) 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^{−3}). 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 quartercritical surface. Initial values of temperature, $ T e = T i = 500 $ eV, $ L q = 50 $ μm, and $ Z \xaf = 50 $ were chosen to be in rough agreement with the data in Fig. 2 for an incident intensity value of 10^{15} W/cm^{2}. 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 collisionalradiative modeling capability which is available only in a very rudimentary form in Chicago at present. For this reason, the ioncharge state in the following simulations is fixed. The timestep is chosen so that $ \tau / \Delta t \u223c 50 $ where $ \tau = 2 \pi / \omega 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 $ \Delta x \u2243 \lambda / 30 $, and the plasma is represented by $ \u223c 500 $ electron macroparticles per cell, with fewer gold ion particles required. Electronion 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 selfconsistent, 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^{15} W/cm^{2}. 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 raytracing photons and actual EM fields. For the pure EM simulation (red trace), we have applied a lowpass filter to the data to obtain effectively “cycleaveraged” 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 $ \omega 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^{15} W/cm^{2}, 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 raytracing approach.
The SRS instability, visibly absent at 10^{15} W/cm^{2}, 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 $ \omega o t \u223c 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 $ \omega \u223c \omega o / 2 $ is consistent with a directly backscattered wave, where $ \omega s \u223c \omega o \u2212 \omega p e \u223c \omega 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 $ \omega / \omega o \u223c 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 ( $ \omega o t < 2000 $), indicates that nonlinear effects are also significant. The FFT of E_{z} at x = −11 μm also shows a prominent forwardscattered SRS peak at $ \omega o + \omega 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 forwardscattered SRS peak in the linear regime. At significantly larger values of L_{q}, there are both strong forward and backward SRS peaks at $ \omega / \omega o \u223c 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 −xdirection and is polarized in the virtual z component, any electric field component in the xdirection is indicative of a plasmon, for which the electric field and wave number are parallel. The figure shows results at $ \omega 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 nonlinear 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 \gamma \u2032 $ 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^{15} to 10^{16} W/cm^{2}. All of the simulations retain the 1/3 μm wavelength, use a fixed $ Z \xaf $ 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.
IV. 2D LPI 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 zdirection with $ \Delta z \u223c \Delta x $. The laser propagates in the −xdirection and is polarized in the zdirection. The laser intensity is also uniform in the zdirection, 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^{2} and L_{q} = 50 μm. The results are shown for $ \omega 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 $ \omega o t = 400 $, as there is clearly a plasmon signal with $ k \u2192 $ parallel to $ x \u0302 $. 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 $ \omega o t = 800 $, some 2D structure begins to emerge. This structure grows in amplitude with time and is seen to be localized around the quartercritical surface at $ \omega o t = 1200 $. Careful inspection of the E_{x} signal shows field structures consistent with localized plasmons propagating at ±45° to the xaxis, which is consistent with the fact that strong TPD growth is predicted in this case. By $ \omega 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 $ \omega 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 quartercritical 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 backscattered 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 nonmaxwellian character, with a high temperature tail becoming visible for $ \omega o t $ values greater than a few hundred. The high temperature tail increases rapidly until about $ \omega o t \u223c 2400 $, after which time the EEDF seems to converge, and the saturation of the instabilities is visible in the timedependence 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 quartercritical 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 quartercritical electron temperature and plasmon field component (top), as well as the upstream and quartercritical 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 quartercritical 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 \xd7 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 \xd7 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 signaltonoise 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^{15} W/cm^{2} for runtimes out to $ 10 4 \omega o \u2212 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.
V. SUMMARY AND DISCUSSION
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 gasfilled hohlraums on hydrodynamic timescales required the development of a multiion PICbased MHD algorithm in order to properly model the details at the interface of the fill gas and the blowoff gold plasma. The collisional absorption of laser energy was modeled by a PICbased raytracing algorithm. The hohlraum simulations using the algorithm discussed here were 1D, but the algorithms are fully generalized for multipledimensions and can be used to model multiple laser beams of arbitrary crosssection and pulse shape. The relatively simple raytracing 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.^{30} 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 \u2264 10 15 $ W/cm^{2}. The presence of a He gas with mass density on the order of 1 mg/cm^{3} may however make the walls less stable to TPD, since shock waves emerging from the AuHe interface can distort the plasma density and temperature profile near the quartercritical 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 \u2273 10 15 $ W/cm^{2} for values of L_{q} consistent with the hohlraum simulations. A nonthermal 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 electronion 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 selfconsistent. Taking these caveats into account, we nonetheless believe that the general conclusions stated above should be valid.
ACKNOWLEDGMENTS
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 DEAC494AL85000.
APPENDIX A: MULTIION MHD FORMULATION FOR PIC CODES
In a previous publication,^{9} we have outlined a PICbased 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 multiion MHD equations, and then briefly describe how they are implemented in the hybrid PIC code Chicago.
1. Formulation of multiple ion MHD
In the following, we briefly outline the procedure used to derive the multiion MHD momentum equations. The derivation given here is similar to the one presented by Glocer et al.^{31} We allow for an arbitrary number, N_{i}, of MHD ion species and a single electron species, with quasineutrality assumed
where n_{e} and n_{k} are the electron and ion densities, respectively, and $ Z \xaf k $ is the charge state of the kth ion. The charge density can then be written as
where
and $ \rho k = ( Z \xaf n k / n e ) $.
The initial nonMHD fluid momentum equations are given for the kth ion species
and the single electron species
where $ v \u2192 $, T, and P, with the appropriate subscripts denoting the species, are the fluid velocity, temperature, and pressure, respectively. The general binary collision frequency, $ \nu \alpha \beta $, for particles of species α and β can be given by the Spitzer value (Chicago uses a general binary form given in Ref. 32). Alternately, the electronion collision frequency can be calculated by the LeeMore model with Desjarlais corrections (LMD).^{19,20} Conservation of momentum implies that $ n \alpha m \alpha \nu \alpha \beta = n \beta m \beta \nu \beta \alpha $. Equations (A4) and (A5) are an extension of the usual plasma twofluid model to allow for multiple ion species.^{32} We have also allowed for an arbitrary equationofstate (EOS) model from which $ Z \xaf $,P, and the internal energy may be obtained as a function of the local densities and temperatures.^{16} We have also included electrothermal force term $ \u223c \u2212 \beta \u2207 T e $ where β is a dimensionless transport coefficient.^{33} 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^{33} 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 straightforward extension to allow for multiple ions and the inclusion of the electrothermal and viscous forces, the fluid energy equations are essentially unchanged from those presented in Ref. 16, which include fluxlimited 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 multifluid equations is similar to the standard singleion formulation.^{34} A generalized Ohm's law is obtained by first multiplying Eq. (A4) by $ e Z \xaf k / m k $ and summing over ion species, k. Equation (A5) is then multiplied by $ \u2212 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 lowfrequency approximation $ \omega \u226a \omega p e $, where $ \omega \u2212 1 $ is a characteristic time scale and ω_{pe} is the electron plasma frequency. We have also neglected the Hall term^{35} which is of order of $ ( \omega c e / \nu 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 ionion collisions, while the third and fourth lines contain multiion MHD correction terms which include the effects of ionion 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 multipleion systems, such as the gasfilled hohlraum simulations discussed in Sec. II. Note that in the limit of a single ion species (N_{i} = 1 and $ \rho k = 1 $) only the first line on the right hand side of Eq. (A9) is nonzero and the ion equation of motion reduces to the standard MHD result.^{9}
2. PIC implementation of multiple ion MHD equations
The multiion MHD formulation described above has been implemented into the hybridPIC code Chicago and is carried out similarly to the singleion 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 \u2192 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 ( $ \u221d T rad 4 $, where $ T rad $ is the radiation temperature) is calculated by the multigroup diffusion approximation,^{21} 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 precalculated 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.
APPENDIX B: RAY TRACING (EIKONAL) EQUATIONS FOR EM WAVE PROPAGATION IN UNDERDENSE UNMAGNETIZED PLASMAS AND PIC IMPLEMENTATION
In this appendix, we describe the implementation of a raytracing algorithm into Chicago. A very general treatment of raytracing methods in plasma physics has been given by Tracy et al.^{30} 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.^{7} A PIC implementation is briefly described below.
1. Electromagnetic wave propagation and dispersion relation for an underdense plasma
In order to introduce some nomenclature used in the following, we present here a brief review of the theory of EM wave propagation in underdense 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^{5} can be included by adding a sink to the equation for the EM energy density.
where
2. Ray tracing (eikonal) equations for EM wave propagating in unmagnetized plasmas and pic implementation
In the treatment of Tracy et al.^{30} the eikonal equations of motion can be written in the form of Hamilton's equation in a 6D $ x \u2192 $ $ k \u2192 $ phase space.
with a dispersion Hamiltonian given by
Using the dispersion function from Eq. (B2) with a spatially dependent electron plasma frequency, $ \omega p e ( x \u2192 ) \u221d n e 1 / 2 $, we obtain the ray equations^{7}
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: $ \lambda \u226a $ all other relevant lengthscales. Moreover, the absorption model assumes that $ \nu e / \omega \u226a 1 $.
3. PIC implementation of ray tracing
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 macroparticle.

Instantaneous photon position.

Instantaneous photon group velocity.

photon energy ( $ \u210f \omega $).
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 \u210f \omega $. 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.^{30} The raytracing model described here behaves as if laser light is spolarized ( $ E \u2192 \xb7 \u2207 n e = 0 $), which is adequate for describing collisional absorption, but for ppolarized light ( $ E \u2192 \xb7 \u2207 n e \u2260 0 $), the component of the electric field parallel to $ \u2207 n e $ can couple to a Langmuir wave. This is a separate mechanism of depositing laser energy into the plasma known as resonance absorption^{5} 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.
APPENDIX C: THEORY OF STIMULATED RAMAN SCATTERING AND TWO PLASMON DECAY
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 smallscale 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 nonmaxwellian 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.^{5}
1. Stimulated Raman scattering
In stimulated Raman scattering, an EM wave ( $ \omega o , k \u2192 o $) is scattered by a electron plasma wave (plasmon), with frequency and wave number matching conditions given by^{5}
where ω_{s} and $ k \u2192 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 $ \omega o \u2273 2 \omega p e $ or when $ n e \u2272 n c / 4 $. For directly backscattered radiation at $ n \u223c n c / 4 $ in which $ \omega o \u223c 2 \omega p e $ and $ \omega s \u223c \omega p e $, the maximum SRS growth rate is given by^{5,25}
where $ v o s = ( e E / m e \omega 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 \u03f5 \epsilon o c E 2 $, where $ \u03f5 = 1 \u2212 n e / n c $ which is $ \u223c 3 / 4 $ for $ n e \u223c n c / 4 $.
Equation (C4) is strictly valid only in a homogeneous and collisionless plasma. The presence of both electronion 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 electronion 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^{24}
2. Two plasmon decay
An EM wave ( $ \omega o , k \u2192 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 $ \omega 0 \u223c 2 \omega p e $ or when $ n e \u223c n c / 4 $. The maximum growth rate occurs when both plasmon wave numbers are at an angle of 45° to both $ k o \u2192 $ and $ v \u2192 o $, in which case^{5}
where γ is again given by Eq. (C4). Equation (C9) also holds for TPD in a collisional plasma.^{37} 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 $ \beta = 36 \u03f5 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 β.^{26}
Equation (C14) leads to an inhomogeneous threshold of
In Refs. 37 and 5,α = 1, while Ref. 26 has $ \alpha = 0.92 $ for $ \beta \u226a 1 $, and 1.38 for $ \beta \u226b 1 $. In the presence of collisions, the TPD growth rate becomes