Simulation and assessment of material mixing in an indirect-drive implosion with a hybrid fluid-PIC code

Hybrid fluid-PIC simulations aimed at a better understanding of the implosion physics and the material mixing into the hot spot are described. The application of a hybrid fluid-PIC code is motivated by the difficulty of modeling the material mixing by the commonly used radiation hydrodynamic simulations. Hybrid fluid-PIC techniques, which treat the ions with the traditional particle-in-cell method, and electrons with a massless fluid, are more adaptable to handle the heating of DT fuel through PdV work and the material mixing near the DT ice-gas interface and ablator-fuel interface of a compressed capsule. During implosion shock convergence, significant reactant temperature separation and a noticeable amount of material mixing are observed, both of which have important consequences for estimating neutron yield and the understanding of implosions. Physical explanations for these phenomena are discussed, with the non-equilibrium effect in the hotspot and hydrodynamic instabilities at the interface as the likely explanation, respectively. The hybrid fluid-PIC method would be helpful to test the phenomenological fluid model describing the material mixing in ICF implosion.


Introduction
Recently, inertial confinement fusion (ICF) [1,2] experiments have reached a record fusion yield of nearly 3.15 MJ [3]. For the first time, the laser fusion experiments have passed the unity fusion conditions that the energy generated by fusion reactions is larger than that was put into the system [4]. However, experiments showed that the implosion performance of capsules is limited by the material mixing inside the hotspot, which is usually seeded by capsule imperfections amplified by hydrodynamic instabilities during acceleration and convergence. The hotspot mix degrades yield by increasing radiative losses and conductive, as well as displacing fuel [5]. Furthermore, the determination of the transport coefficients, e.g., viscosity and diffusivity, is also sensitive to the distribution of ion mix [6]. Therefore, it is crucially important to understand the amount of material mixing in the hotspot region, which is challenging to model and has not been demonstrated yet.
ICF experiments usually generate a broad range of plasma conditions, ranging in temperature of 0.01-10 keV, ranging in density of 10 19 -10 26 cm −3 . Consequentially, ICF experiments often contain a wide range of physical phenomena [7]: non-degenerate to degenerate, kinetic to fluid. The Knudsen number K n , which is a metric to determine whether the kinetic physics plays a role, is usually defined as the ratio of ion mean-free path (λ i ) to a typical plasma scale length (L), K n ≡ λ i /L. In low-temperature and high-density plasmas, K n ≪ 0.01, collisions dominate the physics in the system and it is suitable to describe the plasmas by the hydrodynamic model. The modeling and design of inertial fusion targets mainly rely on radiation hydrodynamic (RH) codes nowadays. However, the conditions are not always satisfied, ion kinetic effects and material mixing outside the scope of RH simulations have been considered as the possible reason to explain the anomalies in fusion yield in recent ICF experiments [8]. While kinetic treatment (e.g., Vlasov-Fokker-Planck methods [9,10] and Monte-Carlo kinetic particle code [11]) may improve our understanding of material mixing and ion kinetic effects in ICF experiments, it falls short of offering shell dynamics, leaving uncertainties compared to recent ICF experimental observations.
Here, we use hybrid fluid-PIC techniques to model the material mixing and ion kinetic effects in an indirect-drive implosion. Its development is motivated by the kinetic plasma effects that are found to play important roles in several regions inside ICF hohlraums. For example, for shock wave propagation across the inner surface of the implosion capsule and in the rebound of this shock wave at the capsule center [12,13], etc. In our hybrid code [14], ions are treated using the traditional particle-in-cell method, while electrons are carefully modeled as a massless fluid. The kinetic treatment of multi-species ions allows an assessment of the material mixing and ion kinetic effects in an indirect-drive implosion. And, it is more reliable and adaptable to handle the heating of deuteriumtritium (DT) fuel through PdV work. However, our code does not have enough ICF capabilities since some required physics input is yet to be implemented at present, e.g., radiation ablation-driven shock wave. To solve this problem and save computational resources, the implosion physics in an early stage of compression is still simulated by the radiation hydrodynamic code RDMG [15,16], which models well the whole ablation process. The remaining deceleration and stagnation phases of the implosion are then simulated with our hybrid fluid-PIC code Ascent-H. Here, the conditions from the RDMG simulation at the deceleration time are linked to the hybrid code Ascent-H.

Simulation setup
This work addresses ion kinetic effects, and hydrodynamic instabilities in the hot spot of an indirect-drive implosion target, using the state-of-art hybrid fluid-PIC code Ascent-H [14]. In distinction to most simulation tools that have been used to simulate the ICF implosion processes, Ascent-H treats the ions with large numbers of computation particles, and the collisions between different particles are modeled with the binary Monte Carlo method [17]. For this sake, the ion viscosity and typical ionic kinetic effects are self-consistently included. Because we mainly focus on the time scales of ions, it is reasonable to treat the electrons as a massless fluid. In Ascent-H, the ion motion equation, electron pressure equation, and electromagnetic fields can be described by: F ie (x i ) in Eq. 1 is the collision force [18], which solves the ionelectron momentum exchange. p e p e (n e , T e ) is defined as the electron pressure term, and n e , T e are the electron number density and electron temperature, respectively. The right-hand side of Eq. 2, P e ∇ · V e denotes the "pdV" term, ∇ · q e for the electron heat flux, and Q e the electron heating. Eq. 3 is the so-called generalized Ohm's law, which provides an equation for the low-frequency electric field in the plasma. Here, η (η ⊥ ) is the resistivity parallel (transverse) to the magnetic field, and R T is the thermal force due to the electron temperature gradient. Note that the self-generated electromagnetic fields in plasma can also be self-consistently considered in our hybrid simulations. To consider the charge separation due to the electric field, the electron density n e can be derived from the Poisson equation: ∇ · E −4πe(n e − Zn i ). The detailed model, denotations, and formulations can be found in Ref. [9].
Here, we discuss 2D indirect-drive implosion simulations following the typical parameters of ICF experiments on the Shenguang laser facility. For the sake of simplicity, it is a cylindrical symmetric CH shell (2D disk) with a radius of 435 μm, filled with deuterium (D) and tritium (T) mixture with an initial ratio of 50:50. The initial thicknesses of the ablator and DT ice are 60 μm, 35 μm, respectively. The initial densities of the CH ablator, DT ice, and DT gas are 1.07, 0.255, 0.0003 g/cm 3 ,

FIGURE 1
Evolution of the hydrodynamic fronts represented by the fluid velocity gradient (|∇u|) for an indirect-drive implosion simulated by the RH code RDMG. The black dash-dotted line represents the neutron yield at the Bangtime (arb. unit). The red solid and dashed lines represent the ablator surface and ice-gas interface, respectively. t 1 , t 2 , t 3 , t 4 represent the beginning of deceleration stage (t = 2.22 ns), shock convergence time (second), Bangtime (t = 2.57 ns), and stagnation time (t = 2.63 ns), respectively.

Frontiers in Physics
frontiersin.org respectively. The peak radiation temperature is about 245 eV. The first stage of the target implosion is simulated by the 1D multi-group radiation transfer hydrodynamic (RH) code RDMG. The shell is imploded by ablation pressure driven by the soft x-rays generated by the interaction of laser beams with the inside wall of Au hohlraum. X-rays heat the shell exterior and vaporize it, creating expanding ablative blowoff outside and an ablator pressure peak that pushes the shell inward. This process can be described by a so-called rocket model [19]. Figure 1 shows the trajectories of the hydrodynamic fronts in the time-space plane, one can see a clear view of the shock dynamics from the shock convergence time to bangtime. The merged strong shock, which propagates through the gas, travels to the center and rebounds at t = 2 ns. The reflected shock then bounces back and forth inside the gas between the target center and the inner surface of the shell. At t = 2.22 ns, the reflected shock wave collides with the imploding shell. This is the beginning of the shell deceleration phase. During this stage, the target converges significantly, and the shell acts like a piston which compresses the central gas and creates a hot spot. Meanwhile, the kinetic energy of the shell is then gradually converted to the internal energy of the hot spot. During the deceleration stage, the ice-gas interface is Rayleigh-Taylor (RT) unstable [20]. RT instabilities might result in an undesirable ion mix at the ice-gas interface which increases energy loss and tends to create a warm-spot instead of a designed hot-spot. Thus, the combination of hydrodynamic instabilities and ion mixing would lead to complex shock behavior whose impact is inadequately captured by the usual RH simulations. Here, we use the hybridfluid PIC code to simulate the implosion physics starting from the deceleration stage.
The initial conditions (density profiles, temperature, and velocity) are extracted from the RDMG simulation and remapped onto the Ascent-H grid at the beginning of the deceleration stage (t = 2.22 ns), as shown in Figure 2. At this time, the DT ice has already been compressed to a high-density spike of 6 × 10 23 cm −3 , corresponding to an increase in density of over a factor 40 near the ablator-fuel interface. The maximum pressure is about 400 Mbar near the DT ice-gas interface. The Ascent-H initial conditions matched the RDMG output with high fidelity. The size of the 2D simulation box is L x 600 μm, L y 600 μm, with 900 × 900, or 1200 × 1200 meshes in the x and y directions in different simulation runs. Note that the width of meshes can greatly exceed the limits of electron Debye length because the electrons in the hybrid scheme are treated as a massless fluid. The simulations were run from 2.22 to 2.8 ns with the time step being about dt = 1.7-2.3 fs. All ions are fully ionized, with the C 6+ ion mass m C 12m p , D ion mass m D 2m p , and m T 3m p , where m p is proton mass. The periodic boundaries are used in the simulation, and we checked that the boundary conditions do not affect the implosion physics in the present cases.

Hybrid fluid-PIC simulation results
An overview of the implosion dynamics from Ascent-H kinetic simulations is shown in Figure 3, which shows the density ( Figures 3A-C), ion temperature ( Figures 3D-F), and pressure (Figures 3G-I) at three typical times. One can see that the outer edge of the dense carbon shell has abundant filament-like structures ( Figure 3A) that are not shown in usual RH simulations. These filament-like structures result from RTIs that are associated with a rarefaction wave. In Figure 3, the first row shows the density, temperature, and pressure at shock collision time t = 2.41 ns. We observe that a peak ion temperature of 3~6 keV is formed at the center, but the density and pressure are still low. This ion temperature approximately matches the measured integrated implosion observables (hot-spot temperature, neutron yield) in the Shenguang laser facility [21,22]. Usually, the ion temperature in the hot spot would be overestimated by the RH simulations due to the mishandling of conductivity, ion kinetic effect, self-generated electromagnetic fields, and so on. Here, in our hybrid simulations, the ion kinetic effect and electromagnetic fields are self-consistently included, which reduced the discrepancies between the simulated and measured ion temperatures.
The second row shows the density, temperature, and pressure at bang time t = 2.57 ns. One can see that a hot spot is created in the center, which is surrounded by the cold DT fuel ( Figure 3E). And, the pressure reaches its peak value although the ion temperature drops down ( Figure 3H). At this moment, the implosion reaches a peak neutron production rate. Although the ablator is affected by the RTI and ion mixing, the temperature and isobaric pressure equilibrium conditions are maintained well in the hot spot region. More importantly, the ice-gas interface is also Rayleigh- Taylor (RT) unstable [23,26] during the deceleration stage of the implosion, which is shown in Figure 3B. The seeds for this RT instability growth include the engineering features, such as a fill tube, shell roughness, and bulk inhomogeneities [24]. RT instabilities result in an undesirable mixing of the cold and hot fuel at the ice-gas interface, which increases energy loss and tends to produce a smaller hot-spot volume. In our hybrid simulation, the radius of the hotspot is only 50 μm, which is much smaller than the RH simulations with less RTI. It is known that the shell then continues to stagnate, and finally disassemble. At the disassembly phase of the implosions, one can see that the target has mode 4 asymmetries at t = 2.69 ns (as shown in the third row in Figure 3), which is affected by the boundary condition. Fortunately, the boundary conditions do not affect the formation of the hot-spot at the bangtime. Figure 4 shows the evolution of the density and ion temperature at the target center of the hotspot. Figure 4A shows that the density of DT plasma increases with time after the shock wave convergence. It increases quite abruptly after around the Bangtime, and eventually decreases again in the disassembly phase. Interestingly, ion temperature reaches a peak value of T i,D~5 keV when the shock wave converges at the center at time t = 2.41 ns. After the rebounding of the shock at the center, the ion temperature drops quickly by approximately 50%. It is known that the reflected shock bounces back and forth inside the gas between the target center and the inner surface of the shell, decreasing in strength after every reflection. At t~2.57 ns, the ion temperature increases again around the Bangtime. By comparing the ion temperature of D and T in Figure 4B, we find a non-equilibrium effect in the hotspot. At the center, the temperature of T is higher than that of D. The physics of this non-equilibrium effect can be explained by the formula below. In the laboratory frame, one can estimate the shock speed and Mach number for different ion species by [25], where γ is the specific heat ratio, α denotes different ion species, and v p is the piston velocity. v 0,α (Z + 1)γk B T/m α denotes the ion acoustic velocity. Thus, one can estimate that v 0,α is higher for the lighter D component. From Eq. 5, we get to know that the shock is Frontiers in Physics frontiersin.org also faster for the lighter D component. From Eq. 6, we get to know that the shock Mach number M α is higher for the heavy T component. Thus, the temperature of T ions behind the shock will also be higher from the thermal equilibrium Rankine-Hugoniot relations [26] here, T 1 and T 2 denote the ion temperature in the upstream and downstream regions. During the shock convergence phase, the temperature separation of D and T in the hot spot volume is most significant, which is consistent with Le's simulation [27]. Assuming Maxwellian reactivity-weighted averages, one can find T T / T D 1.2 (for 50/50 D/T) at the shock convergence time, and T T / T D 1.09 at the Bangtime. In the range of temperatures 0.3-100 keV, the reactivity of DT reactions can be simply estimated by [28], 〈σ]〉 DT 9.1 × 10 −16 exp −0.57 ln T 2.13 64.2 cm 3 s (8) Hence, the total evaluated fusion yield would be overestimated by~15% with 〈σ]〉 DT based on T T compared to T D at the Bangtime.
In order to show the development of the hydrodynamic instabilities during the implosion, we plot the snapshot of the carbon and tritium number densities at the shock collision time (t = 2.41 ns) and the Bangtime (t = 2.57 ns) in Figure 5. The hydrodynamic instabilities at the ablator-DT interface are seen. For the hybrid fluid-PIC studies, the instability can be linked to the Richtmyer-Meshkov instability (RMI) that is produced during shock breakout from the ablator into DT layer. The RMI is barely visible at t = 2.41 ns but, by t = 2.57 ns, has been amplified via interaction with the reflected shock. This instability is also shown in the snapshot of the density of T (see Figures 5B, D). As a result, during the deceleration of the DT fuel layer, ablator material jets can Frontiers in Physics frontiersin.org penetrate the hot spot, which might cool the hot spot and degrade the neutron yield. Note that the hybrid fluid-PIC simulations are prone to creating hydrodynamic instabilities seeds due to ion thermal fluctuation. This fluctuation could be technically reduced by setting more meshes and particles but is limited by computational capability. Since the PIC approach tracks all ion motions, all non-equilibrium physical processes such as diffusion, and dissipation, which are difficult to be accurately captured by the traditional hydrodynamic simulations, are involved naturally in the evolution of RMI. At the late time, the short-wavelength hydrodynamic instabilities interact with each other. While in the hybrid fluid-PIC simulations, the bubble/spike amplitude could be lower due to the ion mixing. The physics of the suppression of RMI could be explained by the increase of viscosity due to the ion mixing. One can estimate the kinematic viscosity by the formula of Braginskii [29], If the DT ions mix into the carbon plasma, the viscosity of carbon will be significantly increased due to the decreasing of number-density-weighted average charge [30] Here, x C and x DT denotes the fraction of the number densities of carbon and DT ions. The effects of viscosity and ion mixing on the growth of RMI have been studied by several researchers. Mikaelian [31], Brouillette, and Sturtevant [32] have shown the influence of a finite interfacial density gradient on the growth of an RM unstable perturbation. The combined effects of mass diffusion and plasma viscosity have been given by Robey et al. [33], Here, _ η is the growth rate, A is the Atwood number, u c is the initial velocity at the interface, k is the perturbation wave number. Note that ψ(k, D, Δt) denotes the RMI growth rate reduction factor, which is due to a temporally increasing diffusional layer thickness of the interface. ψ(k, D, Δt) is obtained by deriving the eigenvalue of the velocity perturbation equation at the interface [33]. According to theoretical estimates, the ion mixing can cause an increase in the viscosity and growth rate reduction factor, both of which can reduce the growth rate of the RMI. For a discontinuous interface, ψ 1. In the present case, our calculations suggest that the growth rate reduction factor can be as large as ψ~5 for k~10 rad/μm. In the ablator-fuel interface, simulations showed an ion temperature of 100eV, and a density of 1~2 g/cm 3 . Given Δt~0.5 ns and x C 0.1~0.9, the kinematic viscosity can be estimated as v s 0.001~0.03 cm 2 /s, resulting in a growth reduction factor by R 1 ≡ exp[−2k 2 v s Δt] 0.996~0.873 for the perturbation wavelength of λ s 1 μm. However, for longer perturbation wavelengths λ s > 5 μm, this growth reduction factor can be negligible R 1~0 .99. Therefore, we get to know that ion mixing tends to reduce the short wavelength modes. In ICF implosion, the mixing of the ablator material mass into the fuel is also an important factor in reducing the neutron yield. Ion mixing between ablator material and fuel is usually not uniform and incomplete. How much mix occurs and how this mix is distributed are crucial to evaluate what the effects are on the neutron yield. To achieve a better quantitative evaluation of the mixing, we draw the product of the T + and C 6+ number densities, n T n C , at the shock collision time and the Bangtime in Figure 6. This product shows the atomically mixed region at the ablator-DT interface which will degrade the performance of capsule implosions. One can see in Figure 6A, the width of the mixing region is very narrow (several microns) since RMI is barely visible at the shock collision time. However, at the Bangtime, as the RMI develops, non-linearity becomes significant and asymmetric spike and bubble structures are generated gradually, which indicates the RMI dominates the ion mixing at the ablator-fuel interface.
The slice of the density profiles also indicates that the width of the mixing region increases over time. The full width of the mixed region, which is defined by 10% and 90% Carbon contours in Figure 7, depends on the density profile. A crude estimate of the width of the mixed ablator-fuel region can be obtained by finding the typical perturbation amplitude, which grows linearly as time goes on: R min ≃ 55.34(t − 2.26)μm before the stagnation time. It indicates that the ion mixing is mainly due to the RMI phase of the perturbation evolution. Furthermore, we checked the roles of electromagnetic fields on the implosion, and found that the impact of these electromagnetic fields is not obvious due to the strong collision processes in the implosion target.
It should be pointed out that the hybrid fluid-PIC simulations are prone to creating hydrodynamic instability seeds due to the statistical noise (finite particle number and grid effect). The mixing ratio of the materials can be related to the statistical noise since the hydrodynamic instabilities are dependent on the instability seeds. Thus, the hydrodynamic instability seeds can be reduced by using a larger number of particles and higher resolution in hybrid simulations, however, most likely they cannot be completely eliminated. This is the characteristic of all particle-based methods, in which the statistical noise is present in the thermodynamic properties of the kinetic simulations. In our hybrid simulations, to reduce the statistical noise, we adopted a fourth-order spline interpolation for the currents and spatial smoothing for the physical quantities.

Summary
In summary, we modeled an indirect-drive cylindrical implosion using a hybrid fluid-PIC model with a kinetic treatment of both the plastic shell and DT fuel ions. The primary advancement of our work is that the injection and mixing between the ablator and fuel are self-consistently included. It is found that RTI at the ice-gas interface leads to an undesirable mixing of the cold and hot fuel, which tends to produce a smaller hot-spot volume and reduce the neutron yield. Furthermore, it is found that the RMI can be developed at the ablator-fuel interface. Although the growth rate reduction is self-consistently included in the simulation, the instability still brings significant ion mixing at the ablator-fuel interface. It is also demonstrated that reactant temperature separation ( T T / T D ≃ 1.09~1.2) in the hotspot region, has important consequences for estimating neutron yield and our understanding of implosions. The hybrid simulation results suggest that there is a margin to enhance the neutron yield in ICF implosion by improving the fabricated target quality to reduce the magnitude of defects that seed the hydrodynamic instabilities. As physics models of hybrid fluid-PIC code become more developed, it could be helpful to further explore plasma effects for ICF implosion.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions
HC, SYZ, and SPZ conceived of the presented paper. HC constructed the hybrid-PIC simulations and wrote this article, FG and ZD constructed hydrodynamic simulations. WZ and BD assisted with the data analysis. All other authors contributed to the discussion of the results.

Funding
This work was supported by the National Natural Science Foundation of China (Grants No. 11975055), and the National Key Programme for S&T Research and Development (Grant No. 2022YFA1603200).