Rayleigh–Taylor instability (RTI)1,2 occurs when a fluid supports or accelerates another fluid of higher density. It plays a crucial role in many natural systems and engineering applications,3,4 in particular in inertial confinement fusion (ICF),5–7 which has recently reached an astonishing milestone toward the goal of ignition,8,9 and in astrophysics.10–12 RTI is an inevitable occurrence and causes severe performance degradation of ICF implosions as the laser energy irradiates the outer surface of a spherical capsule containing the fusion fuel, creating a hotter low-density plasma corona pushing onto the denser and colder shell. Modulations seeded by target roughness and laser nonuniformities can grow via RTI, thereby reducing target compression,13 mixing ablator material into the fusion fuel,14–17 and eventually degrading the neutron yield. RTI at ablative fronts (ARTI) is characterized by mass ablation driven by the laser energy deposited at the critical surface and then transported to the shell surface. It is well known that ablation stabilizes ARTI in the linear (i.e., exponential growth) phase.18–24 When the mode amplitude exceeds a threshold fraction (around 0.1 as predicted by classical RTI theory25,26) of its wavelength, the exponential growth ceases, the nonlinear phase starts, and the modulations develop into “bubbles” (lighter fluid rising through the denser fluid) and “spikes” (denser fluid falling through the lighter fluid) owing to the generation of higher harmonics. The bubble velocity transitions to a constant terminal value27,28 in the deep nonlinear phase of the classical RTI. This transition is commonly known as “nonlinear saturation.” It has been found that ablation will expand the nonlinear saturation threshold fraction from 0.1 up to 0.3 for shorter-wavelength modes according to nonlinear ARTI theory.29,30 In the highly nonlinear phase, ablation will destabilize ARTI bubble growth owing to the vortex acceleration mechanism,31,32 accelerating the bubble to a velocity well above the classical terminal value. The ARTI bubbles are accelerated by the centrifugal force provided by vortices when penetrating through the dense shell, compromising shell integrity in ICF implosions.

Accurate modeling of thermal transport in ICF plasmas is required for predicting the key processes critical to ignition target design, one of which is the growth of ARTI.33 In low-Z materials like deuterium–tritium (DT), where radiation energy is negligible, thermal transport is mainly via electrons. The electron conduction model implemented in general radiation-hydrodynamic codes is based on the classical Spitzer–Härm (SH) theory:34Q_SH = −κ_SH∇T_e, where Q_SH is the SH heat flux vector, T_e is the electron temperature, and ${\kappa }_{\text{SH}}\propto T_e^{5/2}$ is the SH thermal conductivity coefficient. The SH model starts from the following first-order expansion of the electron distribution function (EDF):$$f_e(\mathit{v})=f_0(v)+\frac{\mathit{v}}{v}\cdot {\mathit{f}}_1(v),$$(1)where f₀(v) is the isotropic part of the EDF, which is assumed to be the usual Maxwell–Boltzmann distribution, v is the velocity of electron thermal motion, and (v/v) · f₁(v) is the first-order anisotropic part of the EDF. The heat flux is determined by f₁. The SH model is local, since the heat flux is determined by the local plasma quantities. Unfortunately, this simple SH model fails to reproduce experimental data. The temperature gradients in ICF plasmas are sufficiently large that the assumption of the SH model that the mean free path of high-energy heat-carrying electrons is much smaller than the scale length of the temperature gradient $[L_T\equiv {\left|\nabla \ln (T_e)\right|}^{-1}]$ breaks down in such regions. These high-energy electrons escape steep temperature gradients, causing both classical local heat flux reduction and nonlocal preheating ahead of the temperature front. It is usually suggested that the electron nonlocal heat transport (NLHT) needs to be taken into account when the Knudsen number Kn ≡ λ_ei/L_T ≥ 0.01,35–38 where λ_ei is the electron–ion (e–i) mean free path. Because heat transport modeling affects mass ablation processes, it will also modify ARTI evolution. Discrepancies between experiments and standard SH simulations were observed in a number of pioneering works.39–43 It was shown that simulations without NLHT predicted faster RTI growths in the linear phase than were observed experimentally39–42 and that accounting for NLHT led to considerably better agreement between the simulated RTI growth rates39,40,42,43 and experiments. The NLHT linear stabilization was mostly attributed to a one-dimensional (1D) preheating effect caused by penetration of the high-energy electrons from the tail of the Maxwellian distribution into the high-density target,39,40,42–44 which causes a lower target peak density and a longer density scale length in the longitudinal hydrodynamic profiles. Although the influence of NLHT on the ARTI linear growth rates can be estimated phenomenologically by substituting the 1D hydrodynamic parameters modified by NLHT into the existing linear theories based on the SH model,18–21 its effect on ARTI evolution in the nonlinear phase is unknown, leading to risky uncertainties in ICF designs.

Modeling NLHT is challenging at ICF-relevant temporal and spatial scales. The most accurate thermal transport models solve the Vlasov–Fokker–Planck (VFP) kinetic equations,35 but are usually too computationally expensive in multiple dimensions. Instead, an ad hoc heat flux limiter45 is often employed to take account of flux inhibition effects in hydrodynamic codes to match certain aspects of experimental results. However, a flux limiter cannot capture preheating physics. More recently, the nonlocal multigroup diffusion model proposed by Schurtz, Nicolaï, and Busquet (the SNB model)46 has been considered as a method that balances physical accuracy and computational efficiency. When NLHT effects are significant, f₀ in Eq. (1) will deviate from the Maxwell–Boltzmann distribution, and f₁ will deviate from its counterpart in the SH model. Starting from the steady VFP equation with the Bhatnagar, Gross, and Krook (BGK) collisional operator,47 the SNB model gives a set of transport equations that can be solved for the deviation of f₀ from the Maxwell–Boltzmann distribution. More details of the SNB model are given in Sec. II. SNB-like algorithms have been extensively tested48–51 and used in state-of-the-art ICF codes.

In this paper, we present for the first time a numerical simulation of ARTI evolution with NLHT up to the highly nonlinear phase and uncover a new nonlinear bubble mitigation mechanism due to the reduction of vorticity generation by NLHT. NLHT suppresses vorticity generation by inhibiting the growth of spikes, which in turn suppresses the nonlinear acceleration of bubbles. This newly discovered physical mechanism allows for a better understanding of the effect of NLHT on the evolution of nonlinear ARTI in high-energy-density environments. The results described here can help to understand and evaluate the differences between experimental data and current numerical simulations and lead to better strategies for providing the necessary margin in ICF design.

Our simulations are performed using the hydrodynamic code ART,31,32 dedicated to simulating ARTI in ICF-relevant regimes.31,32,52–54ART solves the single-fluid hydrodynamic equations over a Cartesian grid, with switchable SH and SNB heat transport models being implemented. For each time step, the mass equation, the momentum equation, and the energy equation without thermal transport are first solved explicitly using a finite volume approach:$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathit{u})=0,$$(2)$$\frac{\partial (\rho \mathit{u})}{\partial t}+\nabla \cdot (\rho \mathit{u}\mathit{u})=-\nabla P+\rho \mathit{g},$$(3)$$\frac{\partial E}{\partial t}+\nabla \cdot [(P+E)\mathit{u}]=\rho \mathit{g}\cdot \mathit{u},$$(4)where ρ, u, P, and g are the density, velocity, pressure, and acceleration, respectively. E ≡ [P/(γ − 1)] + ρ|u|²/2 is the total energy, where γ = 5/3 is the ratio of specific heats. The spatial reconstruction is performed using the MUSCL (Monotone Upstream-centered Schemes for Conservation Laws)–Hancock scheme,55 and the Riemann problem at the cell interfaces is approximately solved using the Harten–Lax–van Leer–contact (HLLC)56 solver. The thermal transport part of the energy equation is then treated implicitly to avoid the strict time step requirement of explicit diffusion equation solvers:$$\rho C_v\frac{\partial T_e}{\partial t}=-\nabla \cdot \mathit{Q},$$(5)where C_v is the specific heat at constant volume and Q is the SH or SNB heat flux vector. A pressure correction due to thermal conduction is then made.

The SNB model implemented in ART divides electrons of different energies (mean free paths) into a number of groups and for each group solves a transport equation for H_g proportional to the deviation of the electron distribution from Maxwellian on the given density and temperature profiles:$$\left[\frac{1}{{\lambda }_g(\mathit{r})}-\nabla \cdot \frac{{\lambda }_g(\mathit{r})}{3}\nabla \right]H_g(\mathit{r})=-\nabla \cdot {\mathit{U}}_g(\mathit{r}).$$(6)Here, U_g is defined as$${\mathit{U}}_g\equiv \frac{{\mathit{Q}}_{\text{SH}}}{24}{\int }_{E_{g-1}/k_{\text{B}}{T}_e}^{E_g/k_{\text{B}}{T}_e}{\beta }^4{e}^{-\beta }d\beta ,$$(7)where λ_g is the mean free path, E_g is the upper energy bound of the gth group, and k_B is Boltzmann’s constant. The SNB heat flux is derived as Q_SNB = Q_SH − ∑_gλ_g∇H_g/3, where Q_SH is the SH heat flux calculated on the basis of the given density and temperature profiles. Cao et al.51 found that ∇ · Q_SNB = −∑_gH_g/λ_g, and this relation can be directly employed when solving the thermal transport equation. The SNB module is used to update the hydrodynamic quantities using a predictor/corrector implicit iterative algorithm similar to that of Cao et al.51 and has been validated with Refs. 48 and 57.

To study the growth of small-scale ARTI in the acceleration phase of ICF implosion when the aspect ratio is large, the planar-target approximation is valid and the simulation uses a rectangular computational domain in the x–z plane. For the quasi-equilibrium state in the z direction from which the ARTI perturbations grow, a profile typical of a direct-drive National Ignition Facility (NIF) target is used, as illustrated in Fig. 1(a). The initial ablation front (the interface between the dense and the ablated plasma) is located at z_a and the peak density is ρ_a = 4 g/cm³. The quasi-equilibrium state with the SH model is obtained by integrating the 1D hydrodynamic equilibrium equations in the frame of reference of the shell from the ablation front toward both sides. Our simulations do not include the low-density region where the lasers are interacting with the plasma, and therefore we do not directly handle laser absorption. Instead, the laser energy transported toward the ablation front is simulated by a constant bottom-boundary heat flux calculated self-consistently with the SH model on the basis of the 1D hydrodynamic profiles to ablate the target at an ablation velocity V_a = 3.5 µm/ns (i.e., the penetration velocity of the ablation front into the heavy shell material), and the corresponding ablation pressure (the pressure at the ablation surface) is P_a = 130 Mbar. The quasi-equilibrium state with the SNB model is initialized using the SH hydrodynamic profiles and evolves self-consistently with the SNB model. The SNB and SH heat fluxes are kept the same at the bottom boundary by setting the reflective boundary condition ∂_zH_g = 0. To adjust the strength with which NLHT modifies the 1D hydrodynamic profiles, our simulations are designed with different values of Kn. The characteristic Kn in each simulation is chosen as the maximum local Knudsen number in the system occurring at the bottom boundary, and Kn can be increased by extending the bottom boundary of the simulation box toward the critical surface. The distance L from the bottom boundary (z = 0) to the ablation surface (z = z_a) is set to {38, 60} µm, corresponding to Kn = {0.010, 0.011}. Expanding L from 38 µm to 60 µm not only increases Kn by 10%, but also significantly expands the spatial region in which the nonlocality condition is satisfied (nonlocal effects are usually considered important35–38 when Kn ≳ 10⁻²), since Kn increases slowly with L [see the local Knudsen number profile in Fig. 1(a)]. The larger-Kn case is expected to exhibit more significant NLHT modifications of the 1D hydrodynamic profiles near the ablation front, as shown in Fig. 1(a). A simulation is also performed for another case with shorter L (18 µm, yielding Kn = 0.009), and NLHT is expected to be even weaker in this case.

The gravity is initialized as g₀ = 110 µm/ns². Since the shell mass decreases owing to ablation, the effective acceleration g(t) is slowly and automatically adjusted in time during the simulation to keep the ablation front approximately fixed in space, i.e., $g(t)=[{(P+\rho u^2)}_{\text{lower}}-{(P+\rho u^2)}_{\text{upper}}]/M_s$, where the subscripts “lower” and “upper” indicate the integral values at the bottom and top boundaries, respectively, and M_s is the mass of the remaining plasma in the computational domain. This is equivalent to studying the ARTI growth in the frame of the shell. The RTI is seeded via velocity perturbations at the ablation front in the form V_pz = V_p0 cos(kx)exp(−k|z|), V_p0/V_a = 0.286. The simulations are performed for a series of wavelengths λ ≡ 2π/k. The grid size is Δx = Δz = 0.1 µm. Periodic boundary conditions are applied along the x direction, and inflow and outflow boundary conditions are applied on the upper and lower boundaries, respectively, in the z direction.

The linear growth rates are found to be reduced by NLHT through its enhancement of the ablation velocity V_a due to stronger preheating. The ARTI amplitudes η in the simulations are measured as half of the peak-to-valley heights. The linear growth rate is fitted with data in the exponential-growth stage when η ≤ 0.1λ (below which ARTI is usually considered in the linear phase) for a series of wavelengths [see Fig. 1(b)]. The theoretical growth rate curves for comparison are calculated using $\gamma ={(A_{T}kg-A_T^2{k}^2{V}_a^2/r_d)}^{1/2}-(1+A_T)k{V}_a$,20 where r_d = ρ_l/ρ_h is the density ratio between the effective light (ρ_l) and heavy (ρ_h) fluids, and A_T ≡ (ρ_h − ρ_l)/(ρ_h + ρ_l) is the effective Atwood number. Since the density varies continuously, it is necessary to define appropriate densities for the light and heavy fluids, according to Eqs. (5) and (6) in Ref. 20: ρ_h ≈ ρ_a and ${\rho }_l\approx {\rho }_a{\mu }_0{(k{L}_0)}^{1/\nu }$, where L₀ = L_mv/(ν + 1)⁺¹ is the characteristic thickness of the ablation front, L_m = min[|ρ/(∂ρ/∂x)|] is the minimum density gradient scale length, μ₀ = (2/ν)^1//Γ(1 + 1/ν) + 0.12/ν², ν is the thermal conduction power index, and Γ(x) is the gamma function. The theoretical curves for the SNB cases are calculated using the same Atwood number A_T and density scale length L₀ as in the SH cases, together with the measured V_a in each individual SNB simulation. The fairly good agreement between theory and simulations that can be achieved by adjusting V_a alone is an indication that enhanced V_a is the key factor lowering the linear growth rates. The enhancement of V_a is attributed to the enhanced heat flux and reduced peak density at the ablation front caused by preheating. The growth rate differences between SNB and SH cases are more significant for larger Kn. The linear growth rate in the Kn = 0.009 case is found to be closer to that in the SH case.

In the Kn = 0.01 case, although the linear growth rates are reduced by NLHT, the slowdowns of bubbles and spikes are different. As shown in Fig. 2(a), the SNB spike trajectory deviates from the SH case earlier than the bubble trajectories. By t = 2 ns [i.e., Z_t = 253 µm, where $Z_t\equiv {\int }_0^t{\int }_0^{t^{\prime \prime }}g(t^{\prime })d{t}^{\prime }d{t}^{\prime \prime }$ is the displacement function commonly used instead of time to compensate for slight variations of gravity], when the ARTI amplitude is large enough to be considered to be in the nonlinear stage, the bubble trajectories from the two models are still close, indicating nonuniform ablation of bubbles and spikes by the nonlocal heat flux. In the Kn = 0.011 case, which is not shown in Fig. 2(a), it is found that stronger preheating reduces bubble growth in the linear stage. Therefore, the Kn = 0.01 cases allow us to focus on the NLHT effects on nonlinear bubble evolution while minimizing the impact on bubble linear growth caused by 1D preheating effects on the hydrodynamic profiles.

The ARTI bubble velocity U_b, defined as the speed with which the bubble vertex penetrates through the slab of dense fluid, is plotted in Fig. 2(b) for two different unstable perturbation wavelengths close to the linear cutoff (λ > λ_cutoff ≈ 6 µm) and two different thermal transport models. The bubble velocities start to grow from the values of V_a in each individual case. As the bubble amplitude becomes larger, the low-density plasma filling the bubble is cooled down by the cold walls containing the bubble, and the mass ablation from the bubble walls becomes negligible. At this stage, the bubble behaves as in the classical RTI case without ablation, and the bubble velocity saturates at the two-dimensional (2D) classical value $U_b^{\text{cl}2\text{D}}={[g(1-r_d)/(3k)]}^{1/2}$ given by potential-flow theory,27,28 as shown in Fig. 2(b). The ratio $U_b/U_b^{\text{cl}2\text{D}}$ approaches 1 as ARTI saturates. There is little difference between the values of U_b from the SH and SNB models up to this stage, which is consistent with the fact that the bubble trajectories in the linear phase have not been significantly affected by switching from SH to SNB in the Kn = 0.01 case.

Switching heat transport models alters the bubble reacceleration processes after the first saturation of U_b to $U_b^{\text{cl}2\text{D}}$, as plotted in Fig. 2(b). With both thermal transport models, the simulations show that the vortex carried by the material ablated from the spike is transported into the bubble and accumulates near the bubble vertex as the ablated plasma fills the bubble volume. The vortex then provides a centrifugal force that raises the bubble vertex and accelerates U_b beyond $U_b^{\text{cl}2\text{D}}$, i.e., reacceleration occurs. NLHT causes a reduction in bubble reacceleration. It is found that smaller-wavelength modes are more susceptible to NLHT, because the reacceleration stage is more strongly affected by NLHT in the 7 µm case [Fig. 2(b)]. Simulations have also been performed for larger perturbation wavelengths, and it has been found that the difference between the bubble velocities from the SNB and SH models is much smaller for λ greater than 20 µm.

Reduction of bubble reacceleration with the SNB model is consistent with the weaker vorticity accumulation inside the bubble plotted in Fig. 2(c). The vorticity ω₀ in Fig. 2(c) is the volume average of ∇ × u inside the bubble between 1/k below the bubble vertex and the vertex, as shown in Fig. 3(a). With both thermal transport models, the nonlinear evolution of the bubble velocities can be approximated well by the asymptotic bubble velocity formula given in Ref. 31:$$U_b^{\text{rot}2\text{D}}={\left[\frac{g(1-r_d)}{3k}+r_d\frac{{\omega }_0^2}{4k^2}\right]}^{1/2}.$$(8)A comparison of the simulations with the result from Eq. (8) is shown in Fig. 2(b). This indicates that the vorticity is also the key reason for the reacceleration in the SNB cases. The vorticity growth with the SNB model is slower than that with the SH model, and this effect is more significant for smaller-wavelength cases.

The vorticity is generated in the region where ablation is significant. As the bubble and spike grow, the hot plasma inside the bubble is cooled down by the shell plasma. As a consequence, the majority of the heat flux can hardly enter the inside of the bubble, and mass ablation is predominant at the spike tip. Figure 3(a) shows the long bubble-spike structure. The mushroom-like structures usually seen in classical RTI near the spike tip are suppressed in ARTI owing to intense ablation. Figure 3(b) shows the temperature contours, and it can be seen that the heat flux vectors illustrated by arrows are mainly around and below the spike tip. The evolution of vorticity obeys Dω/Dt ≡ ∂ω/∂t + (u · ∇)ω = −∇P × ∇ρ/ρ², with compressibility and viscosity effects neglected. The right-hand side, known as the baroclinic term, is recognized as the major source generating vorticity in ARTI. The vorticity filling the bubble [Fig. 3(c)] is originally driven by the baroclinic source [Fig. 3(d)] concentrated near the spike tip and is affected by the ablation process there.

The ablation process near the spike is affected by heat transport modeling and influences the falling speed of the spike. The SNB model produces larger heat fluxes in the ablation region than the SH model. In Fig. 4(a), the SNB correction to the heat flux along the axis of the spike in the λ = 7 µm SNB case is plotted. With the SNB model, there is a heat flux increase as a preheating effect in the coronal plasma. The SNB heat flux tunnels further through the high-density spike and causes a 35-fold increase compared with the SH heat flux near the sharp-density-gradient spike interface. Larger heat fluxes near the interface at the spike tip cause greater ablation of the spike and slow down the spike speed. This spike slowdown, which is visible from the linear stage to the deeply nonlinear stage [Fig. 2(a)], is found to be more significant for λ = 7 µm than for 10 µm. This trend is consistent with stronger ablation stabilization of smaller-wavelength perturbations in the linear phase of ARTI.

The baroclinic source is found to be largely correlated with the dimensionless quantity kL_spike in our simulations, where L_spike is the spike height defined as the vertical distance between the bubble and spike interfaces. The volume-averaged baroclinic term ${\dot{\omega }}_{\text{baro}}$ is found to have a linear dependence on kL_spike until the late stage of the simulations, as shown in Fig. 4(b). It can also be seen from Fig. 4(b) that this dependence is insensitive to the heat transport model used. Therefore, a lower spike speed leads to weaker vortex generation, which further reduces the bubble velocity via the vortex acceleration mechanism.

In our 2D simulations with the SH model, both the bubble velocity and the vorticity inside the bubble saturate in the highly nonlinear phase,31,32 as also shown in Figs. 2(b) and 2(c). It can be seen that the NLHT effect in the SNB cases leads to a slower vorticity accumulation inside the bubble, while the bubble velocity may still approach a vortex-accelerated terminal velocity similar to that in the SH case. Given adequate growth time in the SNB case, the slowed-down vorticity generation ${\dot{\omega }}_{\text{baro}}$ grows larger as L_spike grows [see Fig. 4(b)], and thus the vorticity accumulated inside the bubble is likely to increase to a similar saturation value as in the SH case [see Fig. 2(c)].

In summary, it has been shown that NLHT can mitigate ARTI growth in both its linear and nonlinear stages. The strength of the effect of NLHT on the ablation front is varied by adjusting the characteristic Kn. For Kn ∼ 10⁻², the key factor reducing the linear growth rate is the enhanced V_a. The reduction in linear growth rate is found to be due more to slowdown of spikes than to that of bubbles, since NLHT causes greater ablation of spikes than of bubbles. For Kn = 0.01, the linear phase of spike growth is inhibited by NLHT, but the linear phase of bubble growth is minimally affected, and thus study of this regime can help isolate the NLHT effects on bubble nonlinear growth, which is a 2D effect that cannot be predicted by a 1D hydrodynamic simulation. It has been found that NLHT enhances the ablation near the side of the spike tip and slows down the spike. Since the spike height is correlated with the vorticity generation rate, the nonlinear bubble reacceleration due to the vortex generated near the spike tip and transported into the bubble is reduced by NLHT. This newly discovered mechanism is more significant for smaller-wavelength modes and can be utilized to suppress the nonlinear bubble growth of short-wavelength ARTI seeded by the laser imprinting phase in direct-drive laser fusion. The controlled use of NLHT has already attracted some research interest. Otani et al.58 utilized “cocktail color irradiation” to mitigate RTI by moderately enhancing NLHT. Our findings may help optimize those strategies utilizing NLHT to mitigate RTI in ICF-relevant experiments.

In the presence of intense vorticity generation accompanied with ARTI, the self-generated magnetic field is also expected to be intense,59,60 because the self-generated magnetic field is closely related to the vorticity (and the baroclinic term). The combined influence of self-generated magnetic fields and electron NLHT61 on ARTI evolution is an interesting topic. Important related topics, such as bubble merger and competition, convergent geometry, three-dimensional effects, and coupling with radiation transport should be studied in future work.

Acknowledgment. This research was supported by Science Challenge Project No. TZ2016001, by the Strategic Priority Research Program of Chinese Academy of Sciences under Grant Nos. XDA25050400, XDA25010200, and XDB16000000, by the National Natural Science Foundation of China (NSFC) under Grant Nos. U1530261, 12175229, and 11621202, and by the Fundamental Research Funds for the Central Universities. Some of the numerical calculations in this paper were performed on the supercomputing system at the Supercomputing Center of the University of Science and Technology of China.