Kinetic-ion, quasineutral, fluid-electron particle-in-cell simulations of interpenetrating carbon–carbon plasma flows in 2D RZ cylindrical geometry are presented. The simulations are initialized with solid density targets that are subsequently irradiated by 1014 W/cm2 intensity lasers using a raytracing package. The ablation, interpenetration, heating, slowing, entrainment, and stagnation of the plasma flows evolve self-consistently within the code. The particle density, velocity phase space, and fits to the velocity distribution functions are used, along with analytical collisional stopping rates, to interpret the dynamics of the flow evolution. Comparisons to multifluid simulations are described and used to highlight ion-kinetic effects in the setup. Synthetic Thomson scattering diagnostic signals are generated using detailed knowledge of the plasma distribution functions. The large scale of the system, 1 × 1 mm for 2 ns, and the detailed dynamics extracted demonstrate that such hybrid codes are powerful tools for the design and evaluation of laboratory-scale high-energy-density plasma physics experiments.
I. INTRODUCTION
Kinetic-ion, quasineutral, fluid-electron hybrid codes occupy a unique region of simulation capability. On one hand, they model ions kinetically: the ability of the ion distributions to be non-Maxwellian allows the codes to intrinsically capture many effects that are difficult, if not impossible, to capture in a fluid-ion formulation. Examples of phenomena that may require a non-Maxwellian treatment of ions include: flow interpenetration,1–13 ion-viscosity, nonisotropic pressure,14 acceleration of low density ion beams, finite ion gyroradii effects, diffusion and kinetic mix at shocked interfaces,15–18 and nonthermal fusion production.19–22 On the other hand, such hybrid codes are not restricted to resolving electron spatial scales such as the Debye length, the electron skin depth or the electron gyroradius, and their associated time scales. Thus, it is possible to use these hybrid codes to model much larger spatial/temporal scales and use higher dimensionality than is practical using a fully kinetic method. However, we note that when not including kinetic electrons we cannot model inertial electron dynamics or other electron kinetic effects. An additional advantage of hybrid over fully kinetic methods is that they have a well-defined electron temperature, and thus it is easier to interface the system with physics modules such as equation-of-state (EOS), ionization, and radiation transport.
Different flavors of kinetic-ion hybrid codes have been used for many years to capture the details of physics related to ion heating/shock accleration,23,24 magnetized pinches,25–27 space weather,28–31 magnetic fusion energy,32–34 and high-energy-density physics (HEDP).7,35–37 There are many vibrant areas of research in HEDP science that relate directly to kinetic-ion physics. For instance, laser-irradiated, shock-driven capsule implosions show unexpected neutron discrepancies in hydrodynamically scaled experiments38 that may have a kinetic-ion origin.39–41 Due to the importance of high-velocity ions in nuclear fusion, even relatively small deviations from a Maxwell–Boltzmann (MB) distribution can cause dramatic changes in the fusion yield and neutron spectra.19–22 Studying the ion kinetics across a high-mach number shock front is another active area of experimental research42,43 that intrinsically requires a kinetic-ion treatment due to the long ion-collisional mean-free-paths relative to the shock front.
There are many situations specific to inertial confinement fusion (ICF) related to kinetic-ion physics.44 For instance, on the National Ignition Facility (NIF), hohlraums with low gas fill (i.e., 0.03–0.6 mg/cc of He) are found to have advantages when compared to previous designs using >1 mg/cc fill, such as improved coupling from the laser to the hohlraum wall, reduced suprathermal electron generation, and minimized nonlinear laser-plasma interaction effects.45–47 However, at the lowest densities, the standard, single-fluid ICF-design codes48 are found to incorrectly predict the laser-beam propagation, resulting in an inaccurate prediction of the capsule symmetry.49,50 Given the plasma parameters from simulated flows in the hohlraum, a collisional mean-free-path of 0.5 mm was calculated;50 this is non-negligible considering the ∼2 mm separation between the hohlraum wall and capsule. This large mean-free-path suggests the presence of interpenetrating flows. This interpenetration likely requires a kinetic-ion approach to accurately capture the plasma dynamics and, thus, to accurately model laser beam propagation. On the other hand, the large scale of the system, 3 × 5 mm over 10 ns, makes fully kinetic simulations prohibitively expensive. It is precisely in this regime that kinetic-ion, fluid-electron hybrid codes have the opportunity to make a large impact on experimental design as they are able to model full-scale laboratory systems with a more complete physical model than is available in the single-fluid approximation.
Interpenetration relevant to ICF has been studied previously in experiments2,8,11,51 and simulations1,3–7,9,13 using laser-ablated opposing foils. The simulation work has included both multifluid treatments1,3,4,6,9,13 and hybrid kinetic-ion codes, including Vlasov–Fokker–Planck continuum5 and particle-in-cell (PIC)7 approaches. Such work has explored the expansion of the two foils and how the plasma-flow transitions from weakly to strongly collisional. The general consensus is that multifluid simulations do a relatively good job of capturing the general hydrodynamics of the system. However, the multifluid simulations do not capture increased penetration of the fastest ions7 or the temperature anisotropies of the flows.5,7
To investigate the physics of hohlraum-relevant flow-interpenetration on a scale convenient for direct diagnostic observation, a recent experiment was executed on the OMEGA laser by Le Pape et al.52 The experimental setup was azimuthally symmetric, extended around 1 × 1 mm in space, and occurred over 2.0 ns. Like the NIF setup mentioned previously, this experiment is an ideal situation for a multifluid or kinetic-ion hybrid code; to perform such a calculation using a fully kinetic code would be extremely computationally challenging. Conversely, a single-fluid code would not capture the relevant physics of interpenetration.
In this article, we use a kinetic-ion, quasineutral, fluid-electron hybrid particle-in-cell (PIC) code to model this experiment on the OMEGA laser system at full-scale geometry. This model is initialized with a stationary, solid-density, cold, 5 eV, target. A laser irradiates the target via a raytracing package, and the system evolves self-consistently over time and space. As the plasma ablates outward from the opposing surfaces, the plasmas interpenetrate, slow and heat each other, and eventually stagnate. We investigate the evolution of the ion densities and velocity phase space, understanding their progress using collisional stopping and temperature equilibration rates. Finally, we create a synthetic Thomson scattering diagnostic to understand the role of kinetic ions in the experimental observables.
II. SIMULATION PARAMETERS
A. Numerical details
The simulation was performed in a 2D cylindrical RZ geometry using the PIC code Chicago.36 The cell-width was 2.5 μm with a simulation box of 1040 × 4000 μm in R × Z. A time step of 5.3 fs was used. The simulation was initialized with 1600 macroparticles per cell in the solid targets. A macroparticle-splitting algorithm was used to attain adequate macroparticle numbers; every 0.053 ns the macroparticles were split to retain a minimum of 1000 macroparticles per cell. No macroparticle merging was performed. While both the outer ring (OR) and inner rod (IR) are carbon ions, they are described as two distinct species in the model. The laser is modeled via “ray” particles using the geometric optics approximation to include refraction and uses inverse bremsstrahlung energy deposition;53 see Appendix B of Ref. 36 for a detailed description. Symmetric, reflecting boundary conditions were used on the edges of the simulation box. We note that the bulk of the plasma flow does not reach the Z boundary over the duration of the simulation. In the R direction, a shock is launched into the solid material. This shock reaches the R boundary, but the reflection of this shock against the boundary does not reach the critical surface and is thus not expected to change the ablation dynamics substantially.
The simulation method uses a kinetic-ion, hydrodynamic-PIC formulation36,54 that follows the motion of kinetic ions, while modeling electrons as fluids and calculating the electric field via a generalized Ohm's law. The derivation of this method, for the magnetized fluid-ion case, is found in Appendix A of Thoma et al.,36 and a summary of the method for unmagnetized kinetic-ions is found in Sec. III A of Ref. 22. The ion velocity equation uses an unmagnetized, kinetic-ion approach, which neglects thermoelectric effects [i.e., Eq. (15) of Ref. 22]. The ion–electron scattering uses the average velocities of a species and is described using the electron–ion, νej, and ion–electron, νke, collision frequencies via a fluid scattering model.55 Ion–ion scattering is described via a binary collision method56,57 similar to Nanbu and Yonemura's method,58 but with modifications to more accurately handle arbitrary weighted particles.57 Electrons are modeled as a virtual species with attributes (i.e., temperature, density) carried by the ion particles. The electron temperature is advanced via fluid energy equations59 with an electron thermal flux limiter60 of 0.15. Equation-of-state (EOS) tables from Propaceos61,62 are used to model ionization and calculate opacities. EOS properties are calculated for a given ion species independently of the other ion species in the problem. This may not be valid at cold temperature and high density when highly coupled effects are important; in these cases, ion densities and temperatures of all species may be required. However, in this simulation, the ion species mix together only in the hotter, lower-density plasma blow off, where such effects are negligible. The ionization state of the plasma depends only on the electron temperature and is not influenced by mix. We include radiation power loss using a quasi-optically thin approximation. While not used here, we note that the code contains a binary nuclear fusion algorithm63 that may be useful in some applications. The hybrid code can also include a magnetic field via a generalized Ohm's law, but this was not included in the present simulations. From previous work,64,65 we expect that the setups modeled here may attain strong, order MG, magnetic fields via the Biermann battery effect (i.e., ). However, this work also shows that the fields are confined to locations near the target surface and, thus, should not play a major role in the interpenetration. The magnetic-field-driven filamentary Weibel instability has also been observed in interpenetration experiments in the past.66–68 However, experiments performed at higher density, as in the case of our present study, show that increased collisionality reduces the importance of this instability in the interpenetration dynamics.22,68
B. Simulation geometry
We set up the simulation to resemble the recent OMEGA laser experiment.52 This experiment was designed to attain plasma parameters (e.g., velocity, density, temperature) similar to those found in a NIF hohlraum, situated in an easily diagnosable system. The simulation setup, as shown in Fig. 1, consists of an inner rod (IR) made of high-density carbon, 3.5 g/cc, surrounded by an outer ring (OR) of graphite carbon, 2.0 g/cc. The plasmas are initialized at a 5 eV temperature. The IR has a radius of 0.6 mm and the OR has an inner radius of 1.6 mm, for a separation distance of 1.0 mm. The lasers have an angle of 30° with respect to the initial surface and have a super Gaussian spatial profile with an exponent of 4.5 and a full-width-half-maximum (FWHM) of 300 μm. The laser wavelength is 351 nm corresponding to a critical density, , of 8.9 cm−3. A total of 1.9 kJ (1.0 kJ) of laser energy is incident on the OR (IR) representing 19 (10) individual lasers of 100 J spread over the surfaces. The laser duration is a flat-top 1 ns pulse; this corresponds to a power of 1.9 (1.0) TW and an azimuthally averaged peak intensity of around () W/cm2, on the OR (IR) surface.
III. SIMULATION RESULTS
A. 2D evolution of the flows
The simulations begin to evolve when the lasers are incident on the solid targets. We show the laser intensity at 0.5 ns into the simulation in Fig. 2(a). In these images it is clear that the lasers are incident at 30° to the target surface and that there are four independent laser sources input into the simulation. The subplots Figs. 2(b)–2(f) show the ion density multiplied by the average charge state, . We note that in this simulation the carbon ions are quickly ionized, so that = 6 for almost the entire simulation. The top portions of the subplots show the density of the IR and the bottom show the density from the OR. At the earliest time of 0.5 ns, Fig. 2(b), we see the laser ablation of material expanding outward from the initial surface. This ablation is caused by the laser heating electrons via inverse bremsstrahlung. These heated electrons create a strong electron-pressure gradient, , that accelerates the ions. At later times of 0.75 and 1.0 ns, shown in Figs. 2(c) and 2(d), the plasmas expand toward each other and begin to pass through the opposing flow. At this point, the densities of the plasmas are low enough, and their velocities are fast enough, that they do not stop against each other significantly. However, as the plasmas continue to expand, each of the expansion fronts reaches deeper into the other flow and, thus, interacts with an ever-higher density. Eventually, this causes the fastest portion of the flow to slow and then stop against the opposing, higher-density flow. This effect is seen at times of 1.25 and 1.5 ns, Figs. 2(e) and 2(f), at the location of the plasma flows that are furthest from their original surface. As the opposing plasma expands, it pushes back against the incident plasma and we find that the extent of the furthest expansion is pushed back at 1.25 ns as compared to 1.5 ns. Eventually, at later times, this leads to a stagnation of the plasma flows at the midplane between the two surfaces.
While the interaction of the plasma flows occurs along the R dimension, we can see the importance of including the Z dimension from the density plots. First, the expansion in the Z direction is significant. This is due to the large ratio of the foil separation, 1 mm, to the FWHM of the laser spot, 0.3 mm. With this modest aspect ratio, the density falloff will vary between a 1D, , and a 2D, , expansion profile. The 30° angle of the laser is also important to include, given that this will influence the energy deposition and the location at which the laser is reflected. Additionally, the electron temperature will be decreased as compared to 1D, given the larger volume occupied by the plasma for a given energy. Such effects can be included to some extent using a 1.5D formulation,8 requiring some assumptions about the transverse expansion of the flow. Another 2D effect is seen as the plasma flow stagnates on and becomes entrained in the opposing flow. This is seen in the IR flow of Fig. 2(f) beginning at Z = 0 mm, R = 1.2 mm and curving out toward Z = 0.5 mm, R = 1.6 mm. Here the curvature of the IR stagnation is due to the spatial expansion of the OR flow and the spatial variation of the IR radial velocity. Even the center region of the flow at Z = 0 will be affected by the 2D nature of the interaction as at stagnation the plasma will not be confined and will be allowed to expand transversely. This will create a lower stagnation temperature in a 2D setup than would be expected from 1D due to the extra degree of freedom.
B. Evolution of velocity phase space
To study the evolution of the plasma in detail, we plot the phase space (radial velocity VR vs radius R) of the ions originating from the IR and OR as a function of time in Figs. 3(a)–3(g) and 3(h)–3(n), respectively. These plots include ions centered at 0.050 mm. Additionally, we plot a radial velocity distribution function for IR (solid) and OR (dashed) flows in Figs. 3(m)–3(s) at the diagnostic location shown in Fig. 1(b) (i.e., mm, ), which is offset by 0.1 mm toward the OR from the midplane between the two foils. The distributions are normalized to the IR peak. The OR distribution is multiplied by a factor printed on each plot.
The interaction of the two opposing flow can be understood as a transition from weakly collisional, to strongly collisional, to stagnation and entrainment of the flows.5 These stages occur roughly around , 0.9–1.0, ns, respectively.
During the weakly collisional stage, 0.8 ns, as seen in Fig. 3(a), the main feature is the increase in velocity with distance from the initial target surface. This is due to the characteristic self-similar expansion,69,70 where the plasma velocity, , is proportional to the sound speed, , and to the distance from the initial surface, x, divided by the time, t, since the laser irradiation. Here mi is the ion mass, Te is electron temperature and is the adiabatic index.
At 0.9 ns, the plasma becomes more collisional, as depicted in Fig. 3(b). The characteristic increase in velocity away from the initial surface is still evident close to the initial surface, but at distances further from the initial surface, greater than 1.3 mm for the IR flow, we see the plasma has begun to slow and heat against the opposing flow. At the furthest distance, around 1.5 mm, the outermost edge of the expanding IR plasma has slowed completely and is being pushed back with the opposing OR flow. The same features of heating and slowing are observed later in time at 1.0 ns, Fig. 3(c). At this time, at the position furthest from the initial surface, around R = 1.4 mm, the IR flow is noticeably redirected backward by the opposing OR flow and is now considerably reduced in the width of its velocity spread. This occurs as the expansion fronts reach deeper into the opposing flow and thus interact with higher density that stops and reverses the original direction of the flow, a feature we call “entrainment.” This entrained flow then equilibrates with the temperature of the opposing flow, as evidenced by its reduced velocity spread. Concurrently, ions from OR flow, shown in Fig. 3(j), undergo the same dynamics, but symmetrically flipped across the midplane at R = 1.1 mm. This symmetry is expected as OR and IR are composed of the same material and both are irradiated with similar laser intensities. The symmetry is not perfect, considering the difference in intensity and volume elements given the radial nature of the setup.
As time progresses, more lower-density plasma becomes entrained in the opposing flow. At 1.1 and 1.2 ns [Fig. 3(d) and 3(e)], this flow reversal is very clear. In Fig. 3(d) at R = 1.35 mm, for instance, there are two overlapping populations of the IR flow: one that is hot with positive velocity and one that is cold with negative velocity. Later in time, for instance at 1.4 ns, the bulk of the flow at distances far from the initial surface has turned backward, so that the flow is converging at the center between the two surfaces at R = 1.1 mm and the plasma is stagnating with the directed flow energy becoming thermalized in a stationary plasma.
At the diagnostic location of R = 1.1 mm, Figs. 3(m)–3(s), we see the OR distribution progress over time by continuously slowing down, from around 800 km/s at 0.8 ns to around 350 km/s at 1.4 ns, and by continuously becoming wider. The OR ions can be well described by a single Maxwellian–Boltzmann (MB) distribution. On the other hand, the lower-density IR ions significantly deviate from a MB distribution. In the first panel, Fig. 3(m), at 0.8 ns, the IR ions are traveling at positive velocities and are mostly MB-like. By 1.0 ns, Fig. 3(o), the IR distribution has begun to skew, developing an elongated tail toward the negative velocities of the OR flow. At later times, Figs. 3(q) and 3(r), double-humped distributions have developed that show a hotter positive velocity population, in addition to a colder, negative velocity population. Finally, Fig. 3(s), the IR plasma has been almost completely entrained in the OR flow and matches closely the velocity and temperature of that opposing flow.
C. Studying the flow evolution with fitted Maxwell–Boltzmann distributions
In our descriptions of the ion-distributions evolution, we noted that they seemed similar to Maxwell–Boltzmann (MB) distribution functions at some points, perhaps capable of being described by single or double MB distributions. To better quantify this behavior and to study the evolution of the plasma over time, we fit a summation [i.e., ] of two 1D MB distribution functions, , to the IR radial velocity flow, vR,
where mi is the ion mass, ni is the ion density, T is the ion temperature, and is the average radial velocity. We denote the two MB fits as IR+ and IR*, where we define IR+ as the population moving at a faster velocity in the positive direction and IR* as the flow moving in the more negative direction. The fitting was performed using least squares minimization to the find the parameters ni, , and T for both IR+ and IR* fits that best fit the simulation data. We show examples of these fits in Fig. 4, where the data from the simulations are shown as the solid curve, the IR+ 1D MB is shown a single-dotted dashed curve, the IR* 1D MB is shown as a double-dotted dashed curve, and the summation of the two 1D MB fits is show as the dashed curve. The times plotted illustrate deviation from a single MB distribution but are found to be well described using double MB fits. In addition to the double MB fits, we also define average values, denoted , that are defined using the moments of the 1D radial velocity distribution of the IR flow in the simulation, f(vR), where the average velocity is defined as
and the 1D radial temperature is defined as
We note that the notion of temperature may not truly be valid here as the ions deviate strongly from MB distributions. These moment averages are also used to describe the OR flow, . At this location, which is closer to the OR initial surface, the OR flow is much higher density and can, in fact, be well-described by a single MB fit.
We show the resulting parameters of these fits as they evolve over time in Fig. 5. These plots show the average values for the flow (dashed curve) and the flows (solid curve), and the double MB fits representing the IR flow, IR+ (curve with diamonds), and IR* (curve with circles). The double MB fits are only shown for times greater than 0.95 ns because at other times a single MB fit is adequate. Figure 5(a) shows ion densities. The diagnostic location chosen to examine these parameters [shown in Fig. 1(a)] is closer to the surface of the outer ring, and for this reason, the OR flow has higher density; thus, the OR density is divided by five when plotted. Both flows increase at early times, then IR flow plateaus at around cm−3, and finally, at later times the IR flow density decreases as it is pushed back by the opposing OR flow. The ratio of densities between the OR and IR flows is relatively constant for much of the evolution, with the IR flow being around 20% of the OR density. This 20% is interesting in terms of interpenetration because there are competing mechanisms pulling at the particles. On one hand, the OR flow is higher density, which causes a drag on the overall IR flow. However, as we have seen in Fig. 4, the IR plasma is able to retain two distinct MB populations that, summed together, fit quite well to the data. This is evidence that the collisions between a population are strong enough to hold the distribution together to some extent. If the species ratio was much lower, for instance, only a few percent, then the collisions within the flows would likely be completely overwhelmed by the opposing flow and it would likely not be possible to fit the data well with a double MB fit.
We look at the velocity evolution of the flows in Fig. 5(b). Here, in addition the MB fits, we plot the expected velocity from the self-similar (SS) solution, , using an average value of Te = 750 eV for Cs = 245 km/s. As shown in Fig. 6, the electron temperature is slowly varying with time, with values ranging from 650 to 900 eV in the time of interest. The self-similar solutions are plotted in Fig. 5(b) as single- and double-dotted dashed curves for the OR and IR flows, respectively. The velocity evolution shown in Fig. 5(b) shows the flows decreasing in velocity magnitude over time. At early times, less than 0.9 ns, we find that most of this decrease is due to slower-moving flows arriving later in time as expected from the SS expansion. At later times, the IR velocity decreases dramatically compared to the SS expectation. The OR flow velocity magnitude also decreases compared to the SS expectation, but this is less dramatic as it is the higher-density flow. The strong decrease in velocity of the IR flow occurs near the same time as the strong differentiation between the IR+ and IR* flows. Relatively quickly, from 1.0 to 1.2 ns, the IR* flow drops from 630 to –340 km/s, a change of 970 km/s. Over the same time, the IR+ velocity decays from 790 to 210 km/s, a change of only 580 km/s. The strong difference in the slowing can be understood by the strong dependence of slowing down on the relative velocity, , between the two flows. This strong dependence of stopping on relative velocity is one large effect that may demand the use of a kinetic-ion code as compared to fluid-ions. As compared to the average flow, the evolution of the two populations (+ and *) are dramatically different, despite both having evolved from the same single MB distribution around 1 ns. We find the IR* flow becomes relatively quickly entrained in the OR flow, while the IR+ population does not completely stagnate even by 1.5 ns.
The evolution of the temperatures is shown in Fig. 5(c). At early times, the IR* flow increases in temperature. This is due to the frictional heating of the flow, where the kinetic energy of the directed flow velocity is isotropized. This heating of the IR* continues until 1.15 ns, at which time the IR* temperature drops sharply and equilibrates with the OR temperature. The IR+ temperature follows a similar evolution, but it takes more time for the flow to heat and to subsequently equilibrate with the OR flow. These two populations deviate significantly from the average temperature at around 1.2 ns, given that both the + and * populations have a temperature of around 10 keV, while the has around 24 keV. This highlights the fact that using a single temperature fit breaks down strongly at this point; essentially the average “temperature” is measuring velocity spread between the two populations and is unrelated to a thermal MB distribution.
D. Collisional stopping and temperature equilibration rates
To get a better quantitative evaluation of the flow evolution, we examine the collisional stopping and heating rates from Jones et al.55 pertaining to MB distributions. Here the stopping rate, , for a particle α, in our case an ion, is given by
where is the relative velocity, and the summation is taken across all ion and electron species. The stopping rate, , is
where is the charge state, e is the electron charge magnitude, is the Coulomb logarithm, m is mass, is reduced mass, ϵ0 is the vacuum permittivity, n is density, and gives the dependence on temperature
Here is the thermal velocity. The value of asymptotes to unity is .
The heating equation consists of a heating term due to the conversion of flow energy into thermal energy, and a temperature equilibration term
The temperature equilibration rate, , is
with giving the dependence on the relative velocity
where, here, asymptotes to zero as .
In Fig. 7(a), we plot the ion, , (solid curve) and electron, (dashed curve), stopping rates and the ion, (single-dotted dashed curve), and electron, (double-dotted dashed curve) temperature equilibration rates for the IR* population. The same quantities for the IR+ population (e.g., ) are shown in Fig. 7(b). In both cases, it is the ion collisional terms that are strongest. As seen in Fig. 5(b), the stopping power, , of the IR* flow is much stronger than that of the IR+ flow, and thus consistent with the faster slowing observed in IR*.
Recall the abrupt peak of the temperature of the IR* flow in Fig. 5(c) that occurred at 1.15 ns. This peak can be understood as a balancing of the heating and equilibration terms in Eq. (7). At early times, less than 1.15 ns, the ion heating term dominates, causing the flow to heat. At later times, after 1.15 ns, the ion temperature equilibration term dominates, causing the IR* plasma to cool to the temperature of the OR flow. We see this dynamic in Fig. 7(a) as becomes equal to at the time of 1.15 ns. The reason for the rapid increase in the equilibration rate is its strong dependence on the relative velocity between the two flows, as seen in . This strong dependence means that only once the flows have gotten relatively close in velocity can they equilibrate; once this requirement is met, equilibration happens rapidly.
Finally, we note that this system is very dynamic and one must be aware that by observing the same location in space we are not necessarily observing the same particles. For instance, plasma moving at 500 km/s will move 0.050 mm in 0.1 ns; this is a considerable amount, being 1/20 of the total system size. Thus some of the evolution of the parameters is also due to the plasma flow.
E. Multifluid comparison
Sections III A–III D have shown different aspects of the ion distributions that have a kinetic nature. For instance, the splitting of the IR distribution into forward and backward moving distributions can be clearly seen in Figs. 3 and 4. However, previous work7 has shown that bulk hydrodynamic parameters may be similar between a multifluid and hybrid kinetic-ion approach. To understand these differences in more detail, we have used the Eulerian multifluid capacity of the code Chicago to compare with the kinetic-ion results. This method treats the particles as fluids and remaps them into a single particle per species in each cell after each time step.59
We show a pseudocolor plot of the IR (OR) ion density in Fig. 8(a) [Fig. 8(b)] at 1.2 ns. Here the multifluid and kinetic-ion results are plotted in the upper and lower portions of the figures, respectively. Noticeably, the multifluid and kinetic-ion densities are very similar throughout the majority of the simulation space. The major difference is at locations furthest from the initial target surfaces. Here the expanding plasmas do not extend as deeply into the opposing plasmas in the multifluid results. Also, there is a sharp gradient at this point in the multifluid results, whereas in the kinetic-ion case the plasma is more diffuse. Profiles of the flows at mm are shown in Fig. 9. Here the solid (dashed) curves represent kinetic-ion (multifluid) results with red and black showing the IR and OR flows, respectively. From the density profiles in Fig. 9(a), it is clear that the flows in the kinetic-ion simulations penetrate deeper than in the multifluid case. Additionally, the multifluid runs show a stagnation of material at the boundary of one flow (e.g., around 0.85 mm for the OR flow) that leads to a density increase at this point. These differences are due to the fact that the multifluid simulations must use an average stopping power for the entire ion population in a cell, whereas the kinetic-ion simulations use scattering coefficients that vary for each particle dependent on its individual velocity. As the slowing down frequency varies strongly with ion velocity, , the fastest ions in the distribution function will be slowed considerably less than the average. This feature allows the fastest ions to penetrate deeper into the opposing plasma. The same effect is responsible for the increased average velocity of the flows at deeper locations in the opposing flows as shown in Fig. 9(b), at R = 0.85 mm for the OR. In this case, the furthest reaching multifluid OR flows are shown to completely match the flow speed of the opposing IR flow, whereas the kinetic-ion flows slow down and then later increase in relatively speed. A similar aspect is seen in the temperature profiles in Fig. 9(c). The incoming OR and opposing IR flows equilibrate at R = 0.85 mm in the multifluid simulations. However, in kinetic-ion simulations, the OR temperature is reduced at R = 0.85 mm, but does not equilibrate with the IR flow. At distances closer to the initial IR surface, the temperature increases again with decreasing R. We note that here the temperature of the kinetic-ion simulations corresponds to the 3D temperature; i.e., , where is the mean velocity and is the normalized velocity distribution function.
To better understand the kinetic-ion behavior of the OR flow at this region, around R = 0.85 mm, we look at slices of the OR ion radial-velocity distribution at different radii in Fig. 10; distributions at radii of 0.90, 0.85, 0.80 and 0.75 mm are shown as cyan, red, green, and blue, respectively. None of the distributions at any radii resemble a Maxwell–Boltzmann (MB) distribution. Thus, the notion of MB temperature at these locations is not valid. Instead, the temperature is better thought of as a measurement of the variance of the velocity distribution. In this case, there are two components: a positive-velocity colder distribution, this is the colder, “entrained” flow discussed previously—and a hotter negatively flowing distribution that is able to penetrate deeper into the opposing flow. Thus the measurement of temperature is related the relative number of ions in either component. At 0.9 mm, the positive and negative flows are of similar number and thus the temperature is high. At 0.85 and 0.8 mm, the cold positive component dominates; this gives a lower temperature and a more positive velocity. Finally, at 0.75 mm the relative number of ions in the hotter negative component increases with respect to the colder positive component, and thus the temperature rises again and the velocity becomes more negative. Through this understanding of the velocity distributions, the behavior of the temperature and velocity shown in Figs. 9(b) and 9(c) can be understood.
Another method to evaluate the deviation of the velocity distributions from MB is to look at the temperature of the flow along different dimensions. Here we define the 1D temperature as used in Eq. (3). We plot the 1D temperatures in the R, θ, and Z dimensions as the magenta, green, and blue curves in Fig. 11. The standard, 3D-temperature is shown as a black curve; this value is equivalent to an average of all dimensions. The profiles show that at locations close to the initial ablation surface of a given flow, the temperatures are identical. Thus at these locations, the higher-density components (e.g., the IR flow at R < 1.0 mm), a multifluid assumption should be valid. In fact, at these locations we see the best agreement between the multifluid and kinetic-ion simulations in Fig. 9. At locations further from the initial ablation surfaces, the 1D temperatures deviate significantly from the average. Thus, at these locations, the multifluid assumption is not strictly valid. The reason is that the flow of interest is in the minority at this location, and the ion collisions are dominated by the majority flow. The multifluid assumption is based on strong collisionality within a fluid to maintain a MB distribution within that flow. This assumption breaks down when collisions are dominated by an exterior flow.
We have included this multifluid comparison in part to address the question: “are multifluid simulations good enough?” Unfortunately, the unsatisfying answer is “it depends.” First of all, the difference between single- and multifluid methods is dramatic, as shown early on.1 So in comparison with single-fluid, the multifluid results are significantly closer to the kinetic-ion simulations. Looking at the results from Fig. 9, we find, in agreement with Ref. 7, that for much of the radial extent, the multifluid and kinetic-ion simulations agree quite well. On the other hand, there are locations where the simulations disagree completely, notably where the ions undergo deeper penetration that results in a higher temperature. Such differences may be important in ICF simulations. For instance, when increased penetration depth of a high-atomic-number material leads to increased mix and potentially enhances radiative cooling. Additionally, even where the ion temperatures appear to agree (e.g., R = 0.9 mm), the kinetic-ion distribution may not be Maxwellian and thus the use notion of temperature may be misleading. As we have seen in Fig. 11, the temperature may also vary depending on the angle. Getting back to the original question, it is clear that multifluid simulations are not strictly valid in this experimental setup. However, they give results similar to the kinetic-ion results in much of the simulation. Given that multifluid simulations are much faster computationally, we believe that kinetic-ion and multifluid simulations are complementary; both can go much further in understanding interpenetration than conventional single-fluid codes.
F. Synthetic Thomson scattering spectrum
To better understand the experimental observations expected from this setup, we run a synthetic Thomson scattering calculation including the parameters from the three populations (i.e., , IR*, IR+) shown in Fig. 5. For this calculation, we assume a 4ω (263.25 nm wavelength) probe laser and a probed k-vector that is along the direction of radial flow velocity. The scattered spectrum depends on the susceptibility of each ion species and the electrons, and thus it encodes information about the density, temperature, and velocity of each species.71 Small wavelength shifts (<1%) are sensitive to the ion populations, and red (blue) shifts correspond to velocity along (counter to) the k-vector. Because the susceptibilities are complex, scattering from multiple ion populations does not in general produce a superposition of the scattered spectrum from the individual populations, so the scattered spectrum must be calculated for the entire system simultaneously.
Figure 12 shows a temporally streaked image of the synthetic Thomson diagnostic. The positive wavelength shifts, corresponding primarily to the OR flow, are dominant due to the much higher density, about 5×, of this flow. The negatively shifted component corresponds to the IR flow that merges into the OR flow over time. To see more detail, we plot 1D profiles of the signal for different times in Fig. 13. Like the 2D image, we again see the negative shifted IR component merging into the OR feature over time. In Fig. 13, we plot both a triple MB fit (i.e., , IR*, IR+) as a solid curve and a double MB fit (i.e., ) as dashed-dotted curve. We see the strongest deviation between these two at early times, with the most pronounced effect at 1.2 ns. As seen in the simulations, at this time IR* and IR+ flows have similar densities and thus give the impression of a temperature in the flow that is clearly not a MB distribution. At later times, for instance 1.4 and 1.5 ns, the differences between the two signals decrease as the IR flow itself is better described by a single MB fit. The differences between the two synthetic profiles do not seem dramatically different because the signal level of the IR flow is overshadowed by the OR signal. Nonetheless, such signals could be observed in an experiment and understanding how to fit the data is important for understanding the underlying dynamics of the system. This is specifically true when using Thomson scattering because it is generally interpreted using a forward-fitting routine that assumes a physical scenario (e.g., single flow, double flow, Te = Ti, etc.). For instance, using the fits from Fig. 5 we see that if we use a single MB fit for the IR flow, we would expect an ion temperature of 24 keV at 1.2 ns in the diagnostic location. However, in reality this temperature is better understood as being influenced by two separate flows IR+ and IR* that have ion temperatures around 12 and 6 keV, respectively. Thus, a poor choice of physical scenario can result in a gross misunderstanding of the physical parameters of the experiment. Additionally, if the two interpenetrating flow were not of the same element, it is possible that the lower-density material would scatter more strongly. This would amplify the signal from the IR flow and thus increase the difference between double and triple MB fits.
IV. SUMMARY
We have performed a large-scale, millimeter by nanoseconds, 2D RZ cylindrical simulation using a kinetic-ion, quasineutral, fluid-electron, particle-in-cell code. This setup modeled an experiment of an inner ring surrounded by an outer rod, both of which were irradiated by nanosecond-duration, 1014 W/cm2 intensity lasers. The simulation began with solid density targets and used a PIC-based raytracing method to ablate plasma to create two interpenetrating flows. The evolution of the flow parameters was studied, fit with distribution functions, and evaluated using collisional stopping rates to understand the dynamics of the interaction. In this geometry, the flows varied strongly over space and time. At early times, the fastest portions of the flows were found to strongly interpenetrate and pass through each other with minimal stopping. At later times, the fastest portions of the flows reached deeper into the higher densities of the opposing flows, causing the flow to heat strongly and then become entrained and equilibrated in temperature with the opposing flow. This originally fastest portion of the flow then travels with the opposing flow.
The strong interpenetration found in the simulations shows that the single-fluid assumption is wholly invalid for these conditions. On the other hand, the ion distribution functions can be well fit with double Maxwell–Boltzmann distributions, suggesting the potential for using a multifluid method. We have compared with a multifluid method and found that, like previous work,7 this technique captures many of the same behaviors as the kinetic-ion method. The main differences in the kinetic-ion being an increased depth of penetration and a non-MB distribution throughout the simulation. The importance of including such effects will depend on the problem at hand, and thus kinetic-ion simulations can also be used to justify the use of methods, such as multi- or single-fluid, for a given setup.
The work presented in this article is one of the largest-scale uses of the kinetic-ion, quasineutral, fluid-electron particle-in-cell method in high-energy-density science to date, providing evidence that this method can be applied to the design and interpretation of experiments at full laboratory scale. By modeling an OMEGA experiment that was specifically designed to attain similar plasma parameters to an indirectly driven NIF hohlraum, we have shown the ability of the method to capture the collisional physics required for NIF experiments. Using this method, we hope to soon model full-scale NIF experiments. The use of this method to understand the role of interpenetration and kinetic effects in NIF hohlraums, which may be responsible for the lack of capsule symmetry at low gas fills, appears promising. By understanding and modeling low gas fills, we may better optimize NIF geometries to achieve better performance.
ACKNOWLEDGMENTS
We thank S. Le Pape and L. Divol for useful discussions regarding the Thomson scattering diagnostic and A. Link for helping with numerical aspects of the simulations. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory (LLNL) under Contract No. DE-AC52-07NA27344. Support was provided by the LLNL Laboratory Directed Research and Development Program (No. 17-ERD-060). Computing support for this work came from the LLNL Institutional Computing Grand Challenge.