Ion kinetic effects on the evolution of Richtmyer–Meshkov instability and interfacial mix

Plasma effects, such as the multi-component kinetic diffusion and self-generated electromagnetic fields, are recognized as a pivotal key to understanding the physics of interface evolution in inertial confinement fusion and supernova remnants. In this work, a two-dimensional hybrid fluid-PIC code is used to investigate the ion kinetic effects of the single-mode Richtmyer–Meshkov instability (RMI) at the interface between hydrogen plasma and carbon plasma. After an electrostatic shockwave passing through the perturbed interface, the RMI, which reshapes the interface, grows via the vorticity depositing as well as the self-generated magnetic field. After scaling the growth of the interfacial mix region with time, the density transition layer has been found to exceed the disturbance wavelength and lead to a suppression of the instability evolution finally.


Introduction
The Richtmyer-Meshkov instability (RMI), which was primitively investigated by Richtmyer [1] and Meshkov [2] in the 1960s, arises when a powerful shock passes through an interface between different materials and resulting in the initial interfacial perturbations growing continuously. Over the past decades, the RMI has been proven to play significant roles in astrophysics and inertial confinement fusion (ICF) [3,4]. The generation and evolution of this interfacial instability, which refers to baroclinic torque and deposited vorticity, is usually reducing an interfacial mix and resulting in a transition to turbulence at the late time. In the previous studies on the interfacial mix induced by the RMI, attention is mainly paid to hydrodynamic effects, though some important questions remain unsolved. In supernova remnant problems [5], a strong shock generates due to the abruptly halted core collapse rebounding and travels outward through different shells, after the powerful and luminous stellar explosion. In ICF cases, after the soft x-ray interacts with the surface of the fuel pellet, the shock wave generates and transmits the interface ablator and the code fuel [6,7]. In both scenarios, the evolution of the RMI and the interfacial mix under extreme conditions is understood incompletely by now.
A lot of numerical and analytical works have been carried out to examine this issue. However, in single-fluid simulations, the mixing area is simply assumed to be proportional to the amplitude of the turbulent interface because the mass transfer across the interface of different materials is forbidden [8]. Then, the multi-fluids code, which allows for the interfacial transfer of mass, momentum, and energy, is used to analyze the interface statistics in the RMI-related problems [9,10]. However, the hydrodynamic simulations cannot include the ion kinetic effects, which have been found to play important roles in ICF studies [11]. For instance, when a strong shock passes through the interface between the ablator and the fuel, it leads to a far higher level of ion mix via kinetic processes, in comparison with that predicted by the hydrodynamic perturbative transport model [12]. As a consequence of ion kinetic effects, the density transition layer at the plasma interface, which is larger than that given by hydrodynamic simulations, might bring some important influence on the evolution of hydrodynamic instabilities.
In this article, using the hybrid fluid-PIC code, the influence of the ion kinetic effects on the properties of the strong shocks and the evolution of the RMI are studied. In this model, a strong electrostatic shock wave is generated via a fast plasma flow colliding with the hydrogen plasma. When this shock passes through the interface between the light plasma and the heavy plasma, it leads to a high level of mix as well as the RMI growth. It is also found that the growth of the mixing region is scaled by different functions at two phases of the evolution of the RMI. Our hybrid simulations indicate that the interfacial mix, which is known as the density transition layer, could suppress the RMI at a late time.

Simulation setup and basic parameters
One of the main difficulties in RMI studies with ion kinetic effects is that the complex physical processes are hard to be modeled. The Knudsen number, expressed as K n = max(λ i /L i ) with the plasma gradient L i = (d log n i /dx) −1 and the mean free path λ i = v ti /v i for the ion species, is a key parameter to evaluate how important the associated kinetic effects are. More recently, quite a lot of studies show the evidence of the large Knudsen number observed in the sharp interface front where the fast kinetic diffusion of multi-component plasma dominates before a transition of traditional hydrodynamic diffusive process [13,14]. To investigate the ion kinetic effects and the evolution of these sharp interfaces, a series of simulations were performed with a hybrid fluid-PIC code called ASCENT-H [15]. In the ASCENT-H code, the ions with large mass are described as the particles while the electrons with tiny mass are considered to be the fluid. Such a model combined the PIC codes and MHD codes is appropriate for the subsequent assumptions [16]: (a) the dimensions L of the simulation box are much larger than the Debye length (L λ De , λ De = 0 T e /n e e 2 ), which is regarded as a quasi-neutrality system; (b) the time of a simulation is much longer than the relaxation time of electron-electron collisions (T τ ee , τ ee = (4π 0 ) 2 2n e e 4 ln Λ ), which is recognized as a quasi-electron-thermal-equilibrium system, where T e is the electron temperature, n e is the electron number density, m e is the mass of electron and ln Λ is the corresponding Coulomb logarithm. In the ASCENT-H code, the ion-ion collisions are handled by a binary collision model with the Monte-Carlo method [17,18], and the energy exchange between ions and electrons is calculated from the Jones' model [19]. The framework of this hybrid model has been shown in reference [15] in detail.
The hybrid simulation parameters have been shown in figure 1(a). These simulations are performed on a box with a length of x L = 50 μm and a width of y L = 8 μm. At the same time, each micron is divided into 50 meshes while 49 computational particles per simulation cell are used to represent ions. Both the field and particle boundary conditions are periodic, either in the x or the y direction. The temporal evolution of energies of electron fluid and ion particles has been shown in figure 1(b) clearly. Generally, the simulation box is regarded as an energy-conserving system, because the loss of total energy is negligible comparing to the energy exchange between ions and electrons. The left interval x < I 1 contains the piston light plasma (H + ) with an initial electron density n e = 4.9 × 10 23 cm −3 and a drift velocity u 0 = 0.01c. The light plasma (hydrogen H + with n e = 4.9 × 10 23 cm −3 ) is located at the interval 0 < x < I 2 and the rest interval I 2 < x < x L contains the heavy plasma (carbon C + with n e = 4.9 × 10 23 cm −3 ). In this case, the Atwood number is about 0.85. The initial planar interface of the shock is located at x I 1 = 20 μm while the second cosinusoidal, single-mode perturbation on the interface is located at x I 2 (y) = x L /2 + a cos(2πy/y L ) with the amplitude a = 2 μm. The pressure for both sides of the disturbance interface I 2 is at equilibrium with the same initial temperature T = 0.1 keV of both the electrons and ions. In such an initial condition, the Debye length of electrons is λ De = 1.1 × 10 −4 μm which is much smaller than the simulation box, and the relaxation time of electron-electron collisions is τ ee = 7.0 × 10 −2 fs which is much shorter than the simulation time. Thus, electrons are reasonable to be assumed as a fluid. In the present simulation, the run time is about 992τ 0 , where τ 0 = 2π/ω pi,H = 6.75 fs is the plasma frequency of hydrogen ions at the initial condition.

Evolution of RMI
At the early stage, the shock is generated due to the counter-streaming plasma flows. As described by the theoretical work of the shock-generation problem, the shock is affected by the properties of the plasma flow injection [20,21]. With the drift velocity u d = 0.01c in our simulation, the Mach number is about M = u d /c s 21.7, where the ion acoustic speed is about c s,H + = 4.6 × 10 −3 c for hydrogen plasma. In this process, the mean-free-path of counter-streaming plasmas [22] is cm = 82.7 μm which is much larger than the thickness of the shockwave. In other words, the ion-ion collisions between upstream and downstream plasmas are negligible while the electromagnetic fields interacting with ions can be the dominant mechanism of dissipation for super-sonic flow [23]. As a result, this shock can be classified as an electrostatic shockwave which allows for the possibility of particle trapping and reflection effects [20]. Besides, the ratio of the temperature gradient at the interaction region and the mean-free-path is not large enough, leading to that the classical Spitzer-Harm heat flow mode is no longer valid. It should be treated as the free streaming limit heat flow for electrons q FS e = −f q f [24], where q f = n e T

3/2 e /
√ m e and f = 0.09 in our simulation according to some related hydrodynamics simulations and experiments [25].
The main physical quantities of this strong shock have been shown in figure 2 after being averaged in y-direction at t = 158τ 0 . It is noticed that the electron thermal flux always propagates faster than the shockwave, leading to two significant gradients of electrons pressure (P e = k b n e T e ) and the strong electric fields as shown in figures 2(c) and (d) (the black solid line). As it is known, in plasma, thermal conductivity is determined by the electrons, whereas viscosity is due to the ions [26]. Namely, the spatial density distribution of shock is mainly dependent on ion properties. However, the thermal energy of electrons transports faster than the shockwave due to the high thermal conduction of electrons, then the ions in front of the shockwave are preheated by the electron fluid via electron-ion resistivity and friction force. Except that, the super-sonic up-stream ions can transform to the sub-sonic down-stream ions according to the electric potential, and the kinetic energy of up-stream flows would transform to the thermal energy of down-stream flows at the same time. The distribution of total energy of electrons and ions has been shown in figure 2(b), where the energy of electron fluid is calculated by e = ρ e u 2 e /2 + (γ − 1)P e and the energy of ion particles is calculated by As shown in this sub-graph, some super-thermal ions are passing through the shock, which are named as 'precursor particles'. According to figures 2(c) and (d), it is easy to find that the 'precursor particles' are mainly made of drifting hydrogen plasma. Furthermore, in our simulation, the magnetic field observed is trivial in this stage.
At the moment t = 175τ 0 , the shock wave moved from interface I 1 and arrived at the initial perturbation I 2 while the RMI started to grow. The exposition for the generation of the RMI can be explained as follows. At the initial instant, there is a planar shock wave that travels from the light plasma to the heavy plasma. The transmitted shock wave reshapes as the turbulence after passing through interface I 2 .
Due to the reshaping, the pressure of the shock front is not at equilibrium no longer. The transmitted wave, with the same direction as the inject shock wave, leads to a converging pressure at the peak of turbulence (y ∼ 4 μm) and a diverging pressure at the valley of turbulence (y ∼ 2 μm) in the simulation. On the contrary, the reflected wave, with the opposite direction of the inject shock front, leads to lower pressure at the interfacial crest and higher pressure at the interfacial trough. As a result, there is a trend that the heavy ions are more likely to diffuse into light plasma at the trough while the light ions are more susceptible to the mixing in heavy plasma at the crest.
The ions trajectories have been exhibited in figure 3 with the 2D contoured ion vortex plot. The upper subfigure shows hydrogen ions' orbits, and the lower subfigure records the carbon ions' trajectories while the transparency of these colored lines presents the changing of time from 0τ 0 to 992τ 0 . The auxiliary dashed line marked with I 2L and I 2H represent the interfaces of light plasma and heavy plasma, respectively. In comparison with MHD simulations, there is no clear interface in our hybrid fluid-PIC simulations because of the self-consistently ion mixing at the interfaces. Thus, an artificial interface is defined as a ten-percent cut-off of initial ion number density. At the same time, the area between I 2L and I 2H is defined as the mixing region. According to the tracks of particles, the hydrogen ions travel from the crest to the trough and pass the interface, then they are deflected by the vortex. In comparison with light hydrogen particles, the carbon ions are much heavier and hard to be moved or deflected. After traveling a short trip toward the crest, they are swept by the interface I 2L due to the smaller velocities. In order to evaluate the sensitivity of RMI evolution on the Atwood number, we performed a simulation with an initial smaller Atwood number (A = 0.5), in which the carbon plasma is replaced by the tritium plasma. It is found that the RMI grows slower and the mixing area becomes larger (not shown in the figure). The physical reason is simple: the linear growth rate of RMI is proportional to the Atwood number and the diffusion coefficient is proportional to 1/ρ 0 . On the contrary, RMI grows faster with some macro-structures in the large Atwood number case shown in [15].
Another advantage of the hybrid fluid-PIC code is to aid in understanding the resulting actions of self-generated electromagnetic fields on plasmas with the self-consistent calculations. Particularly, the magnetic field, which is generated by the Biermann battery effect, has been found in ICF hohlraum and attracts many interests recently [27][28][29][30]. In figure 4, we plot the physical quantities of RMI at different moments to indicate the behavior of electromagnetic fields and their effects on the evolution of RMI.
At the initial instant, the magnetic fields generated by the source term ∂B/∂t = − e c ∇n e ×∇T e n e which is called Biermann battery effect. Compared to two-fluid MHD simulations, the gradient of electron pressure is in a dominant position of generalized Ohms's Law which is expressed in reference [15] as where η is the resistivity, R is the momentum transfer term and ∇q ν is the von-Neumann viscosity term. In our simulations, the electric field intensity can reach a peak of 4 × 10 10 V m −1 in the early time. The curl of electromagnetic fields, especially contributed by the electric field, is the source of the vortex for ions. However, the last order of the vortex function of electrons has an adverse direction compared to the second order, which means the fields would suppress the RMI by the effects on electron fluid. In fact, as Bond et al pointed out in their two-fluid MHD work [9], the self-generated magnetic field with the peak intensity at about 70 ∼ 100 T is considered to be dynamically ignorable for unmagnetized 2D RMI problem in the high electron pressure cases where the plasma beta is β = 4πn e T e /B 2 max 2 × 10 3 in our studies.

Statistics of interfacial mixture
Since the vast physical phenomena are determined by the spatial distribution of density and the component of different species, primary attention should be paid to ions interface statistics. For example, in ICF hohlraum, the mixing between the shell and fuel is supposed to be one of the influential causes of degradation and ignition failure [31]. As is discussed in the previous section, there is no traditional interface between the heavy plasma and the light plasma, thus the mix region is recognized as the sum of simulation cells that are filled by both the light and heavy ions. As shown in figure 5, the growth rate of the mixing region is statistically analyzed and can be divided into two different phases. The first phase is the shock disturbing phase, which started at 175τ 0 and ended at 475τ 0 in our simulation. It certainly follows the rule that the mixing area is affected by the length of the interface swept by the shock front. As Keenan et.al proposed in reference [32], the mix layer width is assumed as , where u m ∝ v s is the mixing rate in the x-direction, d 0 is the width and t 0 is the time when the shock front arrives at the disturbed interface. Then, in the two-dimensional case, the growth rate of mix area S can be expressed as ∂S/∂t = u m L I (t). The length of disturbance with the cosine shape can be integrated though it could not be expressed by the elementary function. The numerical result of the mixing growth has been shown in figure 6(a). Thus, a triangle-shaped interface was set up as a contrast simulation which leads to a linear growth at the shock disturbing phase (phase 1). As figure 6(b) shows, at the early time of phase 1, it fits the numerical data pretty well which is proportional to the interfacial length. After integrating u m L I (t) with time, the area of the mixing region can be fitted as S ∝ t 2 at phase 1 for a triangle-shaped interface. However, the growth rate would slow down after the shock passed over the interface at the end of phase 1.
The second phase is the instability growing phase (phase 2), which started at the moment t 475τ 0 when the shock front has passed the whole interface in the present simulation. As it is known, the RMI, which is always being considered as a specific case of Rayleigh-Taylor instability, is driven by shockwaves in an infinitesimally short time on the interface and grow up even in the absence of shock [26]. So, the mixing process is assumed to be due to the particles' diffusion across the interface with the effects of shock front disappearing. In these views, the mixing area S is given by where t 1 is the end time of phase 1, S 1 is the mix area at t = t 1 and ζ s is a fitted parameter to describe the mixing rate. According to the detailed 1D gas-metal diffusion research in reference [33], the coefficient D However, for the complex plasma conditions during the instability growing, it is hard to evaluate the exact diffusion coefficient. So we simplified it as and plot the fitted result of mixing grow rate as the red dash line in figures 6(a) and (b). As we mentioned earlier, the ion kinetic effects can be evaluated by the Knudsen number, which is related to the mean free path as shown [34]  where the coulomb logarithm ln Λ = 23 − ln{Z 3 In other words, the ion kinetic effects could be unimportant in high-density cases. In our hybrid code, the ions are described by the PIC method, and the kinetic effects are self-consistently considered. The higher plasma density could reduce the kinetic effects of the ions and verify the classical RMI. These simulations with different densities verify the classical RMI and show the effects of Knudsen number on the characteristic of RMI [35]. The higher (n e = 4.9 × 10 24 cm −3 ) and lower (n e = 4.9 × 10 22 cm −3 ) density cases with the same Atwood number and piston bulk velocity have been performed to support these views. For the high-density case, some micro-structure grows at both sides at the interface, which is not found in large Knudsen number cases. For the low-density case, the distorting of the interface is not obvious and the diffusion across the interface is dominant. Besides, as shown in figure 7, at the simulation time t = 992τ 0 , the length of the widest mix layer has been increasing to d crest = 3.8 μm. It is much larger than the wavelength of the disturbance interface λ cos = 4 μm in our simulation. At the same time, the mass density transition layers are shown as black solid lines which vary with the locations of the interface in sub-figures of figure 7. In reference [36], the continuous density profile has been applied at the initial interface to lessen the density gradient and the baroclinic generation of vorticity. Here, the plasma interface can be regarded as continuous interface because of the self-consistently included ion kinetic effects as shown in figures 7(a) and (b). As considered in reference [37], the growth of interfacial perturbations can be expressed as where δu is the shock-induced velocity of the interface, ψ is called as the growth reduction factor determined by the density transition layer width and Atwood number. And ψ is calculated from the eigenvalue function [38] If the density distribution with transition layer could be expressed as the growth reduction factor is given by which is larger than the unit for the continuous interface. Specifically, the continuous interface would lead to a suppression of RMI. As a result of the larger width of shock-induced kinetical mix (d L ∝ t), the growth reduction factor is also larger than that in the hydrodynamic scheme (d L ∝ t 0.5 ). In order to show the suppression of the RMI growth by the ion kinetic effect more clearly, we carried out the hydrodynamic simulations with the FLASH code (requested from flash.uchicago.edu), and find that the gradient of mass density becomes larger and many micro-structures appears at the interface in the hydrodynamic simulations. In comparison, the micro-structures of RMI are suppressed in the hybrid-PIC results with the self-generated transition layer.

Conclusion
In summary, 2D3V hybrid fluid-PIC simulations have been carried out to investigate the evolution of RMI and the interfacial mix with ion kinetic effects. It has been found that the deflection of shock front after passing an interfacial disturbance causes the vortex disposing at the interface and that the Biermann battery source term triggers the self-generated magnetic field. With the particle tracking method, the light fluid is dynamically proved to pass through the interface of heavy plasma at the crest location. Furthermore, the self-generated magnetic field only has negligible impact on the motion of ions in our simulations. In order to evaluate the interfacial mixing, the mixing process is divided into two phases, which consists of the shock disturbing phase and the instability growing phase. In the first phase, the mixing rate is proportional to the length of the interface swept by the shock front. In the second phase, the mixing rate is fitted by the square root of time. In our views, the interfacial mix induced by strong electrostatic shocks at the early stage might lead to the suppression of RMI growth and should be considered in ICF-related studies. Furthermore, the influence of initial conditions on the ion kinetic effects of RMI is also an interest topic that we shall be investigating in the future.