The Richtmyer–Meshkov instability of thermal, isotope and species interfaces in a five-moment multi-fluid plasma

Abstract

The Richtmyer–Meshkov instability (RMI) results from the impulsive acceleration of a density interface where either it or the acceleration is perturbed. Density interfaces may arise due to a change in gas species, isotope, temperature or a combination of these. We computationally investigate the effect of interface type on the plasma RMI, which is relevant for a range of applications, including inertial confinement fusion. We simulate the evolution of single-mode perturbed thermal, species and isotope interfaces in an ideal ion–electron plasma using the multi-fluid plasma (MFP) model. We find that, in the MFP model, the evolution of different types of interface differs significantly, in contrast to single-fluid models where they behave similarly if the Atwood number is matched. The thermal and species interfaces produce the most severe response to shock acceleration, experiencing the secondary instabilities and enhanced primary mode growth. The isotope interface evolution is restrained in comparison with the former cases, resembling the response predicted by single-fluid models. The determining factor in the severity of the MFP RMI is the density ratio across the initial interface in the electron fluid, which is unity for an isotope interface. We observe that, as the density ratio across the electron interface decreases, so do the magnitudes of the self-generated fields and consequently the severity of the growth amplification. Generally, the evolution of the RMI with different types of interface becomes more similar as the level of coupling between the ion and electron fluids is increased, characterised by reducing the plasma non-dimensional skin depth.

1. Introduction

Research in thermonuclear fusion has been pursued for over a generation in the hopes of developing a carbon-neutral and plentiful power source that would solve several existential challenges facing humanity. Significant progress towards power generation has been made with many milestones passed and others on the near horizon. Inertial confinement fusion (ICF) is one fusion technology which is predicted to be the first to achieve the next milestone, ‘ignition’. During ignition, self-sustaining fusion reactions consume all inertially confined fuel in the target, this has never been achieved. The ICF experiments have made significant progress towards this goal of ignition, notably achieving a ‘burning’ plasma (Zylstra et al. 2022), where more heat is released from fusion reactions than is input to the fuel target. There are many obstacles to achieving the ignition condition but the presence of hydrodynamic instabilities within the fuel target is a leading problem that must be overcome (Lindl et al. 2014; Nagel et al. 2017; Remington et al. 2019).

The most advanced fusion experiments, in the context of progress towards ignition, are conducted by the National Ignition Facility (NIF). Experiments at the NIF operate in the region between a burning plasma, where fusion has begun but is not self-sustaining, and the ignition point. The ablative shock in these experiments propagates with implosion velocities of 300–400 km s$^{-1}$ (Betti et al. 2015) and produces ablation front pressures in excess of 1 TPa. Peak pressures experienced in the hot-spot due to heating and spherical convergence can reach $1\unicode{x2013}10 \times 10^{15}$ Pa (Remington et al. 2019). The implosion results in a final fuel assembly comprised of a hot-spot density and temperature of 30–100 g cc$^{-1}$ and 3–7 keV (34.8–81.2 MK), respectively, and a compressed fuel-shell density and temperature of 300–1000 g cc$^{-1}$ and 200–400 eV (2.32–4.64 MK) (Betti et al. 2015). Peak self-generated magnetic fields are of the order of $10^{3}$ T (Srinivasan & Tang 2012), which is strong enough to affect the electron heat conduction. Any physical imperfections within the fuel assembly can act as an initiation point for the Richtmyer–Meshkov instability (RMI), the instability of an impulsively accelerated density interface (DI) where the DI and/or flow fields are perturbed (Richtmyer 1960; Meshkov 1969). In ICF the characteristic length scales of hydrodynamic instabilities are of the order of 1–1000 $\mathrm {\mu }$m (Remington et al. 2019).

The highly energetic and violent processes created during the ICF implosion make it difficult to receive telemetry on many of the parameters of interest within the fuel assembly. The problem is further complicated by the requirement for non-invasive diagnostic techniques, so as not to disrupt the implosion symmetry and fusion process. The burn-averaged total areal density, neutron-averaged hot-spot ion temperature and neutron yield are typical measurable quantities (Zhou & Betti 2008). The areal density is a measure of the amount of compression achieved and is a fundamental ICF parameter given from the Lawson criteria. All other parameters are related to the measurable parameters, assuming that the physics of energy confinement and losses is well understood and correctly described by existing theories and simulation tools. The significant limitations on experimental diagnostics mean that researchers have turned to numerical simulations to gain insights otherwise out of reach.

The hydrodynamically equivalent physics program (Lindl et al. 2004) studied the hydrodynamic instabilities in preparation for the first NIF ICF ignition campaign. The study concluded from early experiments on planar, cylindrical and hemispherical geometries that ‘ICF ignition capsules remain in the linear or weakly nonlinear [growth rate] regime’ (Lindl et al. 2004). Lindl et al. (2004) concluded that the ablation front of the target fuel pellet is Rayleigh–Taylor unstable and interface perturbations can seed the RMI within the capsule interior during acceleration and then deceleration of non-ablated material. These instabilities, if severe enough, ‘can cause ablator material to mix into the core and radiatively cool the hot spot, decreasing the hot spot temperature and nuclear yield’ Park et al. (2014).

Despite the awareness of hydrodynamic instabilities as an obstacle to fusion experiments, their significance was underestimated in the lead up to the national ignition campaign of 2009–2012. The instabilities were modelled using simplified analytical theory combined with empirical relations based on experimental data from previous investigations. A review of the first three years of ignition experiments concluded that ‘Current evidence points to low-mode asymmetry and hydrodynamic instability as key areas of research to improve the performance of ignition experiments on the NIF and are a central focus of the Ignition Program going forward’ (Lindl et al. 2014). Since this time, considerable effort has been made to suppress the Rayleigh–Taylor instability (RTI) instability and to a lesser extent the RMI, with a focus on the laser pulse characteristics.

The hydrodynamic instabilities, which are well understood for conventional fluids, are complicated by electromagnetic (EM) effects in the plasma which constitutes the ICF fuel target after the laser drive commences. The RMI is one of the primary instabilities that plagues ICF efforts. The literature on the plasma RMI is limited and lacks a consensus on how to appropriately model the ICF RMI. The common models for the plasma RMI are the ideal magnetohydrodynamic (MHD) model (Wheatley, Pullin & Samtaney 2005; Wheatley, Samtaney & Pullin 2012; Wheatley et al. 2013; Mostert et al. 2015, 2017), the Hall magnetohydrodynamic (HMHD) model (Srinivasan & Tang 2012; Shen et al. 2019) and the multi-fluid plasma (MFP) model (Bond et al. 2017b). The motivation for using a MFP model is its superior grasp of the physics involved, compared with the ideal and HMHD models which are simplified in their treatment of the physics. In essence, the MFP theory captures more accurately the effects of charge separation, Lorentz forces, self-generated EM fields, fluid interactions (electron fluid exciting ion fluid) and high frequency phenomena. Furthermore, it allows for investigation of the fundamental types of interfaces present in a plasma RMI which cannot be distinguished by MHD and HMHD due to their simplifications.

The first theoretical and experimental characterisations of the RMI was by Richtmyer (1960) and Meshkov (1969), respectively, giving the instability its namesake. At the time of writing, the applications driving new research on the RMI are the mitigation of hydrodynamic instabilities in ICF (Hohenberger et al. 2012; Lindl et al. 2014; Bond et al. 2017a,b; Nagel et al. 2017; Remington et al. 2019; Bender et al. 2021), and enhancing mixing in supersonic combustion (Yang, Kubota & Zukoski 1993; Yang, Chang & Bao 2014). Furthermore, as indicated by the reviews of RMI (and RTI) by Zhou (2017) and Brouillette (2002), the RMI is important in many other natural and engineered formats. These formats include: many astrophysical phenomena (Arnett et al. 1989; Arnett 2000), atmospheric sonic boom propagation (Davy & Blackstock 1971), driver gas contamination in reflected shock tunnels (Stalker & Crane 1978; Brouillette & Bonazza 1999), combustion wave deflagration-to-detonation transition (Khokhlov et al. 1999a; Khokhlov, Oran & Thomas 1999b), laser–material interactions including but not limited to micro-fluid dynamics (Lugomer 2007) and micron-scale fragment ejection (Buttler et al. 2012), high energy density turbulent mixing (Bender et al. 2021) and many more fundamental studies investigating solid–liquid and solid–solid medium interactions with lasers and fluid flows.

Bond et al. (2017b) investigated the MFP RMI of a thermal interface for a single-mode perturbation. Simulation results showed the MFP thermal RMI (TRMI) differed significantly from the single-fluid hydrodynamic case, which in the absence of an applied magnetic field represents the same solution of an ideal MHD simulation. This investigation is the precursor to the work presented in the current paper with the same methods implemented by Bond et al. (2017b) reproduced here but for differing applications.

The key phenomena observed in the investigation were

• • formation of a precursor electron shock that accelerates the electron interface prior to arrival of the ion shock;

• • Lorentz force bulk oscillations;

• • electromagnetically driven RTI (ERTI);

• • electron fluid exciting the ion-fluid interface;

• • and decreasing severity of MFP effects for reduced skin depth.

The oscillation of the electron fluid about the ion interface is caused by the electron precursor shock. An EM force attempts to neutralise the plasma charge distribution, however, the electron fluid continually overshoots the distribution that would neutralise the plasma because it lacks physical dissipation. The overshoot generates EM fields, dominated by the electric field, that drive a reversal in motion and further drive the oscillations. At early times the oscillation is dominant and drives a variable acceleration RTI that alternates between stable and unstable configurations.

The electron fluid, owing to its lesser particle mass, experiences orders of magnitude greater acceleration that produce high-wavenumber instabilities much sooner than the fluid comprised of the more massive ion particles. The passage of the ion shock at later times produces a region of positive charge between the shock and the interface that results in EM forces driving the electromagnetically drive RTI in the unstable direction until simulation end. The RTI in the ion fluid results in growth of the low-mode ion interface perturbation to substantially enhance the interface growth. The distortion of the ion shock as it traverses the interface in conjunction with rapidly propagating electron vortices results in triple points.

The works by Bond et al. (2017b) varied the plasma Debye length ($\lambda _D$), a measure of the coupling between the ion and electron fluid. Upon refining the Debye length, many of the MFP effects are reduced. However, a result of the tighter coupling is an increase in high-wavenumber features on the ion interface. All cases of Debye length exhibit self-generated EM fields, the source of the electromagnetically driven RTI and Lorentz force bulk fluid oscillations. These EM fields are more intense in the smaller Debye length case (tighter coupling). Bond et al. (2017b) concludes that the MFP effects have the potential to produce significant deviations from the behaviour of the plasma RMI reported by MHD models.

The methodology employed in the current study is driven by an effort to simplify the ICF RMI to a point where its effects can be understood without the superposition of other physical phenomena confusing the investigation. The ‘other’ physical phenomena we refer to are the radiation transport, laser interactions, multi-phase modelling to include the dynamics of the shell, nuclear reactions, converging geometry, ablation and the process of ionisation. The isolation of the RMI dramatically reduces the complexity of the problem, investigative tools and analyses required for the works proposed. The removal of other physics allows the computational resources on hand to be directed to the physics of the plasma RMI. Therefore, the far more physically accurate, but also more computationally expensive, MFP model is used instead of the ideal and HMHD models to investigate the plasma RMI.

This work studies the three fundamental types of DIs, within an ideal MFP, that can host the RMI and compares them with each other. Section 2.1 will discuss the relevance of the three interface types to the ICF problem. Section 2.2 introduces the ideal MFP model and § 2.3 outlines the numerical tool and non-dimensionalisation used for implementation of the model. Section 2.4 specifies the exact scenarios simulated in this study. The code is validated for the hydrodynamic RMI in § 2.3.2 Some key phenomena which occur in the two-dimensional (2-D) MFP RMI base flows will be isolated in the discussion of the analogous one-dimensional (1-D) base flows explored in § 3.1. A comparison of the 2-D planar RMI interface types will be made in § 3.2 and with a discussion of the key parameters influencing the RMI evolution in § 4 and some comparisons with the hydrodynamic RMI in § 3.2.4. A conclusion of the works completed is presented in § 5.

2. Problem definition and methodology

2.1. Fundamental interfaces

The ICF fuel capsule is a complicated engineered assembly of ablator, dopant and fuel layers that provide many seeds for hydrodynamic instabilities and general asymmetries. Figure 2 shows the key types of perturbations that may seed the RMI and RTI within the ICF fuel capsule. Native interface perturbations and bulk defects in the capsule material are consequences of the processing techniques used for fabrication e.g. the surface roughness due to inhomogeneous solidification of the deuterium–tritium ice layer. Energy coupling non-uniformities and asymmetries are caused by unequal irradiation (indirect or direct drive) of the fuel target. Intentional engineering features, for example the fill tube used in constructing the fuel target that remains embedded, are at significant risk of instability growth initiation. The varying sources of seeds can each excite instabilities at different locations within the fuel capsule. The situation is further complicated by the ionisation of the fuel capsule constituent materials during implosion as this process augments the perturbations.

The fundamental types of interface that have unique responses to the RMI are the: (i) thermal, (ii) species and (iii) isotope interfaces. Figure 1 shows the non-dimensional initial conditions for all interface types studied. To elucidate the effect of each interface type we enforce mechanical equilibrium (pressure through $P=nk_{B}T$ (where P is pressure, n, number density, k _B, Boltzmann constant, T, temperature) and charge neutrality ($n_eq_e + n_iq_i$) across the interface. The (i) thermal interface possesses a discontinuity in thermal energy within both the ion and electron fluids, producing a DI in both fluids. The TRMI may occur at the solid–gas fuel interface, shown in figure 2, seeded by native surface roughness on the solid fuel. The (ii) species interface possesses a is discontinuity in ion particle mass and ion charge, also producing a DI in both fluids. The species RMI (SRMI) exists at the interfaces formed by layers of ablator–dopant and dopant–solid fuel and can be especially detrimental to ignition by introduction of non-fusible material into the hot-spot. Finally, the (iii) isotope interface possesses a discontinuity in the ion particle mass only, i.e. different isotopes of an element, that results in a DI in only the ion fluid since the isotopes on each side of the interface have the same number of electrons in their shells. An interface discontinuity in charge density and thermal energy, and whether these discontinuities exist in both or one of the ion and electron fluids, will produce unique responses to the RMI, despite having the same mass-density ratio. The isotope RMI (IRMI) does not have an immediately obvious place in an ideal ICF fuel target, however, it provides insight into targets with inhomogeneous fuel layers, intended or defective, as well as providing completeness to elucidate key influences in (i) and (ii).

Figure 1. Illustration of the initial conditions for the four interface conditions investigated. Columns form left to right show the electron mass density, ion mass density, electron temperature and ion temperature.

Figure 2. Sources of hydrodynamic instabilities for cryogenic ICF targets of deuterium–tritium (DT) fuel, reproduced with permission from Smalyuk et al. (2017).

2.2. The MFP model

This study is conducted using the ideal MFP model, modelling the plasma as comprised of only the ion and electron fluids, neglecting any neutral species. It is appropriate to use a fluid modelling approach for the ICF plasma because of the high density and temperatures at which the plasma exists (Bellan 2008; Chen 2016). It is appropriate to apply the fluid description of plasma when the hydrodynamical time scale ($\tau _H$) is much slower than the thermal relaxation time scale ($\tau _e$, $\tau _i$ and $\tau _{eq}$) of all charged particles comprising the plasma (Goedbloed, Goedbloed & Poedts 2004). Furthermore, for ideal MFP to be valid we require the dissipative diffusion time scale ($\tau _D$) to be greater than the hydrodynamical time scale. For a typical ICF plasma, the characteristic time scales are $\tau _e \approx 10^{-17} \ll \tau _i \approx 10^{-15} \ll \tau _{eq} \approx 10^{-14} \ll \tau _H \approx 10^{-11} \ll \tau _D \approx 10^{-7}$, and therefore the use of an ideal MFP model is satisfied. We note that, while not the case in this work, the presence of an exceptionally strong magnetic field produces a perpendicular diffusion time scale that is of the same magnitude as the hydrodynamical time scale increasing the error in fluid description.

In the ideal MFP model, the fluid equations are derived by taking moments of the Vlasov equation, where the result ends as the familiar Euler equations with EM forcing terms. A set of fluid equations, conservation of mass; momentum; and energy, exist for each fluid modelled. In our case we have one ion fluid and one electron fluid, yielding two sets of fluid equations. The following assumptions for the plasmas investigated are used unless otherwise noted. The ions and electrons are each in thermal equilibrium with themselves. The velocity distribution for each fluid is Maxwellian and this distribution gives an accurate representation of currents within each species. The plasmas, as per ideal MFP theory, neglect viscous and resistive terms in the momentum equation and conductive heat transfer in the energy equation. This reduction can be appropriate when the time scale of diffusive terms is much larger than that of the advective terms (Goedbloed et al. 2004), which is generally the case in our problem of interest. Furthermore, the inclusion of these terms (that are often significant in effect) may obscure the ideal MFP phenomena that are the focus of the study. The thermal relaxation (resulting from collisions between species) is neglected for this reason, despite the expectation that they will produce a significant heat exchange between species.

We follow the non-dimensionalised ideal MFP model formulation of Bond et al. (2017b) based on the work by Loverich (2003). The non-dimensionalisation reduces the variance of the magnitudes of variables which must be resolved. This procedure is important when the disparate variables must be solved together in a single system. The chosen non-dimensionalisation results in a less stiff numerical system and conveniently allows the plasma regime to be set by only two parameters, the skin depth, $d_S$, and the plasma ratio of thermodynamic and magnetic pressure, $\beta$. The non-dimensionalisation is shown below, note the ‘hat’ symbol and subscript zero indicate a non-dimensional and reference parameter, respectively,

(2.1)\begin{equation} \left.\begin{array}{l@{}} \displaystyle \hat{n} =\dfrac{n}{n_0} \quad \hat{m} =\dfrac{m}{m_0} \quad \hat{q} = \dfrac{q}{q_0} \quad \hat{\rho} = \dfrac{\rho}{\rho_0} \\ \displaystyle \hat{\boldsymbol{u}} = \dfrac{\boldsymbol{u}}{u_0} \quad \hat{p} = \dfrac{p}{n_0m_0u_0^{2}} \quad \hat{\varepsilon} =\dfrac{\varepsilon}{n_0m_0u_0^{2}} \quad \hat{x} = \dfrac{x}{x_0}\\ \displaystyle \hat{c} = \dfrac{c}{u_0} \quad \hat{t} = \dfrac{t}{t_0} \quad \hat{\boldsymbol{B}} =\dfrac{\boldsymbol{B}}{B_0} \quad \hat{\boldsymbol{E}} = \dfrac{\boldsymbol{E}}{cB_0} \\ \displaystyle \widehat{\psi_E} = \dfrac{\psi_E}{B_0} \quad \widehat{\psi_B} =\dfrac{\psi_B}{cB_0} \quad \widehat{d_S} = \dfrac{d_S}{x_0} \end{array}\right\}, \end{equation}

where $n$ is the number density, $m$ is the particle mass, $q$ is the charge, $\rho$ is mass density, $\boldsymbol {u}$ is the velocity vector, $p$ is thermodynamic pressure, $\varepsilon$ is the thermal and kinetic specific energy, $x$ is a length, $c$ is the speed of light, $t$ is time, $\boldsymbol {B}$ is the magnetic field vector, $\boldsymbol {E}$ is the electric field vector and $\psi _E$ and $\psi _B$ represent the wave speed of the divergence cleaning. Further variables to be used temperature $T$, Boltzmann's constant $k_B$, atomic number of a species $Z$, ratio of specific heats $\gamma$, vacuum permittivity $\epsilon _0$, the permeability of free space $\mu _0$ and the hydrodynamic Mach number of the propagating shock $M$. Note that the subscript $\alpha$ indicates some unique fluid i.e. ions or electrons.

The choice of reference values for each property is problem specific. The default choices used in the works completed and proposed are as given next unless otherwise stated. The reference length $x_0$ is set to the typical wavelength of the RMI perturbation of about 1 $\mathrm {\mu }$m, reference charge is set to the elementary charge, reference velocity is set to the electron thermal velocity and the reference number density and reference mass are set to achieve the required $d_S$ and $\beta$. The reference magnetic field strength is set according to $\beta$ (2.14a,b), reference density is set to the product of reference number density and mass, reference time is set according to the ratio of reference length and reference velocity and reference mass is set to a hundred times an electron mass – an artificially large electron mass relaxes the problem stiffness allowing for more tractable solutions. A comparison of typical ICF conditions with the current simulation properties is given in table 1. The current conditions do not exactly mirror those of ICF but are close enough to give insight into the physical phenomena occurring. It is vital to understand that the purpose of these works is to provide detailed qualitative insight on the influence of the RMI and secondary instabilities which affect ICF efforts, not the exact quantitative extent of which these phenomena occur.

Table 1. Comparison of typical ICF parameter orders of magnitude (Betti et al. 2015) and those simulated in the present work.

The non-dimensional forms of the fluid equations, with the hat symbol dropped for brevity, are given in (2.2)–(2.4). The fluid equations are recognisable as the inviscid Euler equations with EM source terms included (gravity is negligible) and exist for each fluid

(2.2)$$\begin{gather} \frac{\partial \rho_\alpha }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot}{\rho_{\alpha}\boldsymbol{u}_\alpha} =0, \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial \rho_\alpha \boldsymbol{u}_\alpha}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot}{(\rho_{\alpha}\boldsymbol{u}_\alpha\boldsymbol{u}_\alpha + p_\alpha \boldsymbol{I})} =\sqrt{\frac{2}{\beta_0}}\frac{n_\alpha q_\alpha}{d_S}(c\boldsymbol{E} + \boldsymbol{u}_\alpha\times \boldsymbol{B}), \end{gather}$$

(2.4)$$\begin{gather}\frac{\partial \varepsilon_\alpha }{\partial t} +\boldsymbol{\nabla} \boldsymbol{\cdot}{(( \varepsilon_\alpha +p_\alpha )\boldsymbol{u}_\alpha)} = \sqrt{\frac{2}{\beta_0}}\frac{n_\alpha q_\alpha c}{d_S}\boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{u}_\alpha. \end{gather}$$

Here, the subscript $\alpha \in (i,e)$ represents the species modelled. We define mass density, pressure and energy density by

(2.5a–c)\begin{equation} \rho_{\alpha} = n_\alpha m_\alpha, \quad p_\alpha = n_\alpha k_B T_\alpha, \quad \textrm{and} \quad \varepsilon_\alpha = \frac{p_\alpha}{\gamma - 1} + \frac{\rho_\alpha|\boldsymbol{u}_\alpha|^{2}}{2}. \end{equation}

Maxwell's equations govern the evolution of the EM fields and are given in non-dimensional form

(2.6)$$\begin{gather} \frac{\partial \boldsymbol{B}}{\partial t} + c\boldsymbol{\nabla} \times{\boldsymbol{E}} = 0, \end{gather}$$

(2.7)$$\begin{gather}\frac{\partial \boldsymbol{E}}{\partial t} - c\boldsymbol{\nabla} \times{\boldsymbol{B}} ={-}\frac{c}{d_S}\sqrt[]{\frac{\beta_0}{2}}\sum_\alpha n_\alpha q_\alpha \boldsymbol{u}_\alpha, \end{gather}$$

(2.8)$$\begin{gather}c\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{E}} = \frac{c^{2}}{d_S}\sqrt{\frac{\beta_0}{2}}\sum_\alpha n_\alpha q_\alpha, \end{gather}$$

(2.9)$$\begin{gather}\text{and} \quad \boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{B}} = 0. \end{gather}$$

Discretisation of the EM divergence constraints, (2.8) and (2.9), introduces error, that is mitigated by implementation of the perfectly hyperbolic form of Maxwell's equations, as per Munz, Ommes & Schneider (2000), resulting in the new Maxwell's equations

(2.10)$$\begin{gather} \frac{\partial \boldsymbol{B}}{\partial t} + c\boldsymbol{\nabla} \times{\boldsymbol{E}} + c\varGamma_B\boldsymbol{\nabla}{\psi_B}= 0, \end{gather}$$

(2.11)$$\begin{gather}\frac{\partial \boldsymbol{E}}{\partial t} - c\boldsymbol{\nabla} \times{\boldsymbol{B}} +c\varGamma_E\boldsymbol{\nabla}{\psi_E} ={-}\frac{c}{d_S}\sqrt[]{\frac{\beta_0}{2}}\sum_\alpha n_\alpha q_\alpha \boldsymbol{u}_\alpha, \end{gather}$$

(2.12)$$\begin{gather}\frac{\partial\psi_E}{\partial t} + c\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{E}} = \frac{c^{2}\varGamma_E}{d_S}\sqrt{\frac{\beta_0}{2}}\sum_\alpha n_\alpha q_\alpha, \end{gather}$$

(2.13)$$\begin{gather}\text{and} \quad \frac{\partial \psi_B}{\partial t} + c\varGamma_B\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{B}} = 0. \end{gather}$$

We define the reference values of skin depth and $\beta$ ratio as

(2.14a,b)\begin{equation} d_{S,0} = \frac{1}{q_0}\sqrt{\frac{m_0}{\mu_0 n_0}} \quad \textrm{and} \quad \beta_0 = \frac{2\mu_0 n_0 m_0 u_0^{2}}{B_0 ^{2}}.\end{equation}

The Lorentz force acting on a single particle is given by $\mathcal {L}=q_\alpha (\boldsymbol {E} + \boldsymbol {u}\times \boldsymbol {B})$, and the $X$ and $Y$ components (important for the 2-D simulation results) are given by

(2.15a,b)\begin{equation} \mathcal{L}_x=q_\alpha (E_x + u_yB_z - u_zB_y) \quad \textrm{and} \quad \mathcal{L}_y=q_\alpha(E_y + u_xB_z - u_zB_x). \end{equation}

2.3. Numerical solver

2.3.1. Description

The numerical implementation used for this work was developed by Bond et al. (2017b). A finite volume method is used to allow for a development of a 3-D modelling capability and to leverage existing knowledge in computational fluid dynamics of EM systems, such as Loverich (2003) and Abgrall & Kumar (2014). The implementation uses a 2-D three-vector architecture, therefore no x- and y-magnetic field components are generated by the fluid motions in simulations, however, there is generation of the z-component of magnetic field. The solution requires a very high degree of spatial and temporal refinement due to the wide range of length scales and advective speeds associated with the plasma regimes modelled. The spatial refinement is satisfied by implementing the system of equations in the adaptive mesh adaptive mesh refinement (AMR) framework, AMReX (Zhang et al. 2019). The trigger for the cell refinement is the relative ion and electron density gradient

(2.16)\begin{equation} \frac{|\boldsymbol{\nabla}{\rho_\alpha}|{\rm \Delta} x}{\rho_{\alpha}} > \frac{5}{100}. \end{equation}

We implement a general consideration for effective spatial resolution such that the smallest length scale is refined by at least two cells, $\delta x \leq {\lambda _{min}}/{2}$, although this is only enforced when the normalised charge density exceeded the small threshold of 5/1000, $n_\alpha q_\alpha > \frac {5}{1000}$ (an emphasis is placed on high resolution in areas of EM features). The temporal refinement is satisfied in some part by manipulating the resolved value of speed of light in the non-dimensionalised parameter space. The manipulation is in choosing a large value for the reference velocity so the non-dimensional speed of light, $\hat {c} = c/u_0$, is smaller, thereby, lessening the temporal refinement required by the Courant–Friedrichs–Lewy (CFL) condition than what is typical in most plasmas. Care must be taken to ensure that the non-dimensional speed of light is still the greatest characteristic speed in the system, otherwise non-physical behaviour due to interaction of fluid and EM waves may occur.

As discussed by Loverich (2003), the homogeneous part of each of the three distinct systems of equations (ion fluid, electron fluid and electromagnetic fields) can be solved independently because the systems are only coupled through the source terms. This allows the option, which is exercised by Bond et al. (2017b), to solve each system using its own Riemann solver. The source terms are then solved using an implicit solver. The time integration is done by a two-stage second-order accurate Runge–Kutta scheme. Piecewise linear spatial reconstruction is used for both cell centred and face centred variables with limiting by a min–mod limiter in the primitive variables. As per the discussion above, the homogeneous solution of the fluid and EM system of equations is solved using Harten–Lax–van Leer–contact (HLLC) solver and Harten–Lax–van Leer–Einfeld (HLLE) solver (Einfeldt et al. 1991; Toro, Spruce & Speares 1994) and HLLE (Harten, Lax & van Leer 1983; Einfeldt 1988) approximate Riemann solvers. The source terms are solved locally using an implicit method as per Abgrall & Kumar (2014).

There exists no analytical test case and solution for the two-fluid model. There are, however, numerous plasma solvers that have been produced and documented in the literature, allowing for comparison of numerical experiments. The standard test for MFP solvers is the two-fluid electromagnetic plasma shock which is an extension of the Brio and Wu shock tube (Brio & Wu 1988). Our numerical solver has been verified by comparing the solution of the EM plasma shock with Loverich (2003).

Verification of sufficient solution convergence is completed by comparing flow statistics of varying effective resolutions and performing a Richardson extrapolation to estimate convergence error. A 2-D refinement study was conducted in the preceding work by Bond et al. (2017b). Flow statistics are used instead of pointwise convergence because the latter is not possible without physical dissipation being modelled to set a minimum physical length scale. The chosen flow statistic is the ion interface width. The average L2 (the square root of the sum of the squared errors) norm calculated for the relative error of the fine grid solution to the Richardson extrapolated solution showed that the calculated error for relevant length scales, $d_D=0.1$ and $d_D=0.01$) was significantly below 1 % when using effective resolutions of 2048 and 4096, respectively. It was observed that these resolutions correspond to those required to resolve the high frequency plasma wave packets in the 1-D base flow. This is sufficient to capture the bulk flow characteristics of the interface evolution. One of the characteristic Debye lengths in the present study, $d_D=0.002$, is not covered by these results, but the 1-D convergence study in § 2.3.3 indicates that the high frequency waves generated are well resolved. This meets the criteria for resolving the interface statistics from (Bond et al. 2019) and gives the authors confidence that bulk flow characteristics are captured in two dimensions at these resolutions. Effective resolutions of 2048, 4096 and 4096 are used for the Debye lengths of 0.2, 0.02 and 0.002.

2.3.2. Hydrodynamic RMI validation

The experiments of Motl et al. (2009) were simulated with Cerberus, run hydrodynamically (gas is uncharged), to verify that the combination of numerical methods implemented is able to adequately capture the dynamics of the RMI. Good qualitative and quantitative agreement was found. The experimental conditions of the scenario (case 8 in the work by Motl et al. 2009) simulated were an interface of helium and sulphur hexa-fluoride gasses at temperature $295$ K and pressure $98$ kPa, initial perturbation amplitude $1.36$ cm, perturbation wavelength $16.74$ cm and Atwood number $\mathcal {A}=0.95$. Figures 3 and 4 show our simulation and the experimental result by Motl et al. (2009), respectively. The simulation reproduces all the qualitative features from the experiment, where both figures show a mushroom has developed at the head of the spike and behind the head, KHI rollers develop behind the interface, and a reverse jet is formed at the base of the spike. Figure 5 shows the comparison of the interface width, $\eta$, for the experiment and simulation, which show good agreement. There are two differences to note between simulation and experimental conditions (i) the vertical shock tube arrangement employed by Motl et al. (2009) and (ii) and the lack of viscosity in the ideal hydrodynamic simulation. The ideal simulations did not model the gravitational force or viscosity which both retard the growth of the instability. These two differences account for the faster growth rate in the simulations. The comparison with the experiment shows Cerberus is able to accurately model the hydrodynamic RMI and secondary instabilities in a hydrodynamic fluid.

Figure 3. Cerberus simulation RMI interface mass-density contours at dimensional times 0.67 ms and 1.77 ms for experimental conditions Motl et al. (2009) for the case 8 of a $M=1.95$ shock.

Figure 4. Time sequence of experimental images for He/$\textrm {SF}_6$, $M=1.95$ case (a) initial condition $t\approx 0.00$,(b) $t\approx 0.67$ ms and (c) $t\approx 1.77$ ms. Reproduced from Motl, B., Oakley, J., Ranjan, D., Weber, C., Anderson, M. and Bonazza, R., 2009. Experimental validation of a Richtmyer–Meshkov scaling law over large density ratio and shock strength ranges. Phys. Fluids, vol. 21(12), p. 126102, with the permission of AIP Publishing.

Figure 5. Comparison of the interface width from Motl et al. (2009) case 8 and Cerberus simulation results for a $M=1.95$ shock in a He/$SF_6$ interface.

The scenarios from Dell et al. (2017) were simulated with Cerberus and compared with the simulation results and empirical relation by Dell et al. (2017). Dell et al. (2017) use a smoothed particle hydrodynamic numerical scheme (Stellingwerf 1991) and derives an empirical relation for the nonlinear and linear initial growth rates. Our simulations reproduced the conditions from Dell et al. (2017) with shock $M=3$ and Atwood number $\mathcal {A}=0.8$. The initial perturbation wavelength was $\lambda =3.3\times 10^{-3}$ m, light and heavy densities were $\rho _l = 1\times 10^{-3}$ kg m$^{-3}$ and $\rho _h = 9\times 10^{-3}$ kg m$^{-3}$, respectively. The simulations were run with reference parameters of $x_0=3.3\times 10^{-3}$ m, $P_0=2.494\times 10^{6}$ Pa, $T_0 = 300$ K, and reference mass equal to mass of the light fluid particle. The particles in the light and heavy regions were assigned non-dimensional masses of $1$ and $9$, respectively, and both were ideal monoatomic gases with $\gamma = 5/3$. The domain length was $40$ in $x$ and $2$ in $y$ with zero-gradient and periodic boundary conditions in $x$ and $y$, respectively.

The ratios of amplitude to wavelength, $a_0/\lambda$, simulated were $\{0.06$, $0.25$, $0.3$, $0.4$, $0.5$, $1.0\}$ all with the conditions specified above. The maximum grid refinement in each case (for the adaptive mesh refinement used by Cerberus) was set to ensure the interface amplitude was well resolved (${\rm \Delta} x \approx a_0/100$) , likewise, the time-step constraint was enforced (more restrictive than the CFL condition) such that the interaction of the initiating shock wave and interface was accurately resolved ($\textrm {d} t\approx a_0/10u_{shock}$). The initial growth rate was calculated as prescribed by Dell et al. (2017), as the slope of the linear regression of interface amplitude data from the time the reflected and transmitted shocks are emitted from the bubble, until the time $0.3\tau$, where $\tau =\lambda /v_\infty$.

After shock traversal of the interface, the interface itself and surrounding fluid have a velocity that is characteristic of the problem. This velocity for a planar interface is $v_\infty$ and is used to define the background motion. It can be derived from linear theory (Richtmyer 1960; Wouchuk 2001) or simulations. Values for $v_\infty /c_l$, where $c_l$ is the sound speed in the light fluid, at the conditions described above from our simulation and that of Dell et al. (2017) are $1.0391$ and $1.0328$, respectively. The value of $v_\infty /c_l$ from linear theory (Wouchuk 2001) is $1.0391$, showing excellent agreement with our simulation and good agreement with that of Dell et al. (2017). Figures 6–8 show the mass density and x-velocity (offset by the background motion $v_\infty$) of the fluid at late time for $a_0/ \lambda =0.06$, $0.25$ and $0.4$. In the frame of reference moving with speed equal to $v_\infty$, theory (Abarzhi, Nishihara & Glimm 2003; Abarzhi 2008) predicts that the fluid motion away from the interface is small, there is significant fluid motion close to the interface, and the nonlinear RMI bubbles flatten at late times. We observe that for low values of $a_0/ \lambda$ the fluid motion away from the interface is very small but as $a_0/ \lambda$ increases, the increased presence of transverse waves disrupts the velocity field away from the interface more significantly, however, compared with the velocity field about the interface, it is still small. Our simulation results reproduce the qualitative prediction from the nonlinear theory.

Figure 6. The mass density (a) and $x$-velocity (b), offset by the background motion $v_\infty$, for the $M=3$, $\mathcal {A}=0.8$ and $a_0/\lambda =0.06$ case. A contour line of the DI is shown in purple.

Our results indicate that as the ratio of initial amplitude to wavelength increases, the growth rate does not increase monotonically. We find that the growth rates increase, reaching a maximum between $a_0/\lambda =0.25\textrm { and } 0.3$, then decreasing such that at $a_0/\lambda =1$ the growth rate is comparable to the $a_0/\lambda <0.1$, in qualitative agreement with the results from Dell et al. (2017). The simulations that enter the nonlinear regime exhibit the characteristic flattening of the bubble. This indicates that the pressure fluctuations around the interface are weak at the time of flattening (Pandian, Stellingwerf & Abarzhi 2017). Results also capture the small-scale dynamics of the problem, such as reverse jets which occur when two flows collide at small angle of attack (Pandian et al. 2017) and are fully immersed in the heavy fluid. Figures 9 and 10 show the bubble speed (in the frame of reference of $v_\infty$ and normalised by $v_\infty$) for $a_0/\lambda =0.06$ and $0.25$. The figures show the bubble speeds and growth after the reflected and transmitted shocks pass the bubble. The bubble speed decreases but is periodically influenced by the passage of transverse reflected shocks. The influence of these transverse shocks is greater for increased $a_0/\lambda$ as inclination of the interface increases, aligning the reflected shocks more with the y-direction, therefore, interacting with the interface for longer time due to the periodic boundary conditions in y. The $a_0/\lambda =0.06$ case does not flatten during simulation run, but we see from figures 7 and 10 that the $a_0/\lambda =0.25$ case does, with a feature developing on the bubble. This feature forms in a similar way to the reverse jets at the base of the spike, except here the curvature of the fluid density profile (in the bulk of the fluid) is opposite, leading to a jet in the negative x-direction. The growth of this jet at the bubble is what produces the increase in bubble speed at late time $a_0/\lambda =0.06$ case, it is has already developed by the times shown in the $a_0/\lambda =0.25$ and $0.4$ cases. We see from figures 6–8 that the x-velocity across the interface when the bubble has flattened is continuous but still discontinuous in the $a_0/\lambda =0.06$ where the bubble still has finite curvature, in agreement with the nonlinear theory for interface with no mass flow (Abarzhi et al. 2003). Overall, our simulation results show good qualitative agreement with the theory, reproducing all the expected features.

Figure 7. The mass density (a) and x-velocity (b), offset by the background motion $v_\infty$, for the $M=3$, $\mathcal {A}=0.8$ and $a_0/\lambda =0.25$ case. A contour line of the DI is shown in purple.

Figure 8. The mass density (a) and x-velocity (b), offset by the background motion $v_\infty$, for the $M=3$, $\mathcal {A}=0.8$ and $a_0/\lambda =0.4$ case. A contour line of the DI is shown in purple.

Figure 9. The growth rate and bubble speed for the $M=3$, $\mathcal {A}=0.8$ and $a_0/\lambda =0.06$ case. The velocity is in the initiating shock direction, the reference frame of the background motion, and normalised by the background motion ($(v - v_\infty )/v_\infty$.

Figure 10. The growth rate and bubble speed for the $M=3$, $\mathcal {A}=0.8$ and $a_0/\lambda =0.25$ case. The velocity is in the initiating shock direction, the reference frame of the background motion, and normalised by the background motion ($(v - v_\infty )/v_\infty$.

The results and comparison with Dell et al. (2017) are given in table 2. Some significant deviation from the empirical relation is found in the small $a_0/\lambda$ range and for the $a_0/\lambda =1.00$ case but good agreement is found for the remainder. The comparison with the simulation results is mostly in good agreement except for the $a_0/\lambda =0.06$ and $1.00$ cases. The simulations results by Dell et al. (2017), which form a subset of the group of simulation results that the empirical model was based on, exhibited deviations from the empirical model that are similar to our simulations. The empirical relations for the linear and nonlinear regimes (Dell et al. 2017) were found to be a useful predictor of behaviour, despite some deviation, and is attractive due to its simplicity and universality, being dependent on $\mathcal {A}$ and $a_0/\lambda$, and being calibrated from simulation data.

Table 2. Comparison of Cerberus results with the empirical relation and simulation results of Dell et al. (2017). *In the case of $a_0/\lambda =0.06$, the empirical relation by Dell et al. (2017) for the initial growth rate $v_0$ was used, whereas the empirical relation for the nonlinear $v_0$ was used for all others.

Good qualitative and quantitative agreement was found when comparing with experiments (Motl et al. 2009) and reasonable agreement was found when comparing with the simulation results and empirical relation derived from a different numerical scheme (Dell et al. 2017), giving the authors confidence that the Cerberus code accurately captures the RMI and secondary instabilities.

2.3.3. The 1-D convergence study

The results of a 1-D test case have been used to determine an appropriate resolution in the streamwise direction. The ideal MFP model, lacking physical dissipation, does not converge to a single solution. As the grid is refined, progressively smaller flow features will emerge meaning no minimum length scale limit exists in this regard. Therefore, the criterion for appropriate resolution was taken to be the resolution of all plasma waves and to pragmatically minimise the relative change between simulation results. The results of successive refinements of resolution is shown in figure 11 for a skin depth of 0.1 (Debye length of 0.002) unit lengths. We see significant improvement between resolution of 1024 to 2048 cells per unit length and continued but diminished improvement between 2048 and 16384 cells per unit length. The expanded view especially shows that a reasonable resolution that captures the flow features well is the 4096 resolution.

2.4. Problem description

The generic problem follows that of Bond et al. (2017b) which was designed to be comparable to previous RMI studies (Samtaney 2003; Wheatley et al. 2005, 2013). Figure 12 shows the initial configuration that forms a three zone Riemann problem consisting of zones zero ($S_0$), one ($S_1$) and two ($S_2$) from the left boundary. The interface between $S_0$ and $S_1$ generates a shock that travels to the right where it encounters the interface formed between fluids in region $S_1$ and $S_2$. The perturbation of the interface between $S_1$ and $S_2$, is chosen to consist of a single period of a sinusoid with an amplitude of 0.1 non-dimensional length ($\frac {1}{10}$ domain width). The perturbation mean location in the x-dimension is 0.2 non-dimensional lengths from the initial shock. The DI studied in this work is a material interface. The initial interface is a stationary contact discontinuity without any mass or heat flux (we do not consider, for example, the kind that would be involved in the RMI of a travelling flame front Ilyin & Abarzhi 2022). Furthermore, there are no phase changes associated with the interface in this work.

As an initial condition, charge neutrality is enforced everywhere and mechanical equilibrium is enforced between $S_1$ and $S_2$ as minimum requirements for initial stability and clarity of results. These two conditions ensure there is no movement of the interface $S_1$ to $S_2$ prior to the shock interaction due to EM or hydrodynamic effects, excluding the shocks involved in the simulation. The interface is not likely to be stable in actual ICF experiments due to an X-ray preheat preceding the shock, shown to exist in shock tube experiments (Keiter et al. 2002; Yamada, Kajino & Ohtani 2019), and drive asymmetries. However, it is useful to model the interface as stable so we may elucidate the fundamental physical phenomena.

Figure 11. (a) Ion and electron number density from a 1-D simulation of a discontinuous interface with incident shock at non-dimensional time $t=0.3$, illustrating grid convergence. (b) Zoomed view of results showing only the electron number density.

Figure 12. An example of initial the initial conditions and developed evolution of the RMI.

The boundary conditions are periodic in the $y$-dimension and zero gradient in the $x$-dimension. The domain is taken to be one reference length in the $y$-dimension, and in the $x$-dimension it is two times the distance travelled at the non-dimensional speed of light during simulation time, to avoid any possible numerical artefacts reflecting off the boundaries and interacting with the observed solution space. The expansive $x$-dimension of the simulation space is made practical by the adaptive mesh refinement used in the numerical solver, where a very coarse base grid of only eight cells across the domain allows a grid to be inexpensively extended away from the region of interest. The boundary between $S_1$ and $S_2$ in figure 12 represents the RMI DI which is chosen to be in the light-to-heavy configuration to avoid complications from a phase inversion (heavy-to-light configuration). A hyperbolic tangent transition function is used between $S_1$ and $S_2$ to impose a smooth transition that makes the solution less susceptible to numerical artefacts and to ensure a consistent interface thickness at different resolutions. This function has the form

(2.17)\begin{equation} f(x) = f_R + \frac{f_L-f_R}{2} \left(1+\tanh\left(\frac{2x}{\eta}\text{arctanh} \left[\frac{9f_R-10f_L}{10(f_L-f_R)}\right] \right)\right), \end{equation}

where $f_L$ and $f_R$ are the variables of interest on the left and right of the interface, and $\eta$ is the width containing 90 % of the transition, chosen as 0.01 non-dimensional lengths.

2.5. The MFP interface modelling

A significant benefit of the MFP model is the ability to investigate fine details of the fluids that comprise the RMI interface types. Therefore, the subtle differences in the thermal, isotope and species interface types can be distinctly simulated. If a single-fluid model, such as ideal MHD or HMHD were used, these details would not be captured and all three interfaces would behave similarly. Consider figure 1, the unique differences across the cases presented would not be captured, in MHD only one homogeneous case of a DI could be used to investigate all cases. The thermal case models a discontinuity in the temperature of elemental hydrogen, ${{}^{1}_{1}\textrm {H}}$. The isotope case models a discontinuity in particle mass in the ion fluid that is produced by the interface of hydrogen, ${{}^{1}_{1}\textrm {H}}$, and tritium, ${{}^{3}_{1}\textrm {H}}$. The first species case is modelled with an interface of hydrogen, ${{}^{1}_{1}\textrm {H}}$, with helium (isotope), ${{}^{3}_{2}\textrm {He}}$, and the second species case is modelled with an interface of hydrogen with a fictitious isotope of lithium, ${{}^{3}_{3}\textrm {Li}}$. We reiterate once more that the scenarios simulated are not accurate quantitative representations of the phenomena occurring in ICF, however, the physical phenomena which are elucidated through these simulations are representative of those which would be exhibited in an ICF experiment.

2.5.1. Parameters

The parameter space in the simulations is large, therefore a subset must be chosen to make the investigation tractable. The parameters used to define the simulated regime were taken from previous investigations (Samtaney 2003; Wheatley et al. 2005, 2013) and are the ion-fluid species mass densities either side of the interface, the ion partial pressure for $S_1$ and $S_2$, species hydrodynamic and EM properties e.g. ratio of specific heats and particle charge, and the shock Mach number. Electron-fluid parameters such as fluid mass densities, pressures, number densities and temperatures (kT); are set according to the ideal gas equation of state, normal shock relations and physical properties of the species involved. The following parameters were set:

(2.18)\begin{equation} \left. \begin{array}{@{}c@{}} m_e = 0.01, \quad m_{i0} =m_{i1} =1.0, \quad \rho_{i1} = 1, \\ q_e ={-}1.0, \quad q_{i0} =q_{i1} = 1.0, \quad \rho_{i2}=3 ,\\ \gamma_e = 5/3, \quad \gamma_{i0}=\gamma_{i1}=\gamma_{i2}=5/3, \quad p_{i1} = 0.5 \\ \textrm{and} \quad M_0 = 2.0. \end{array}\right\} \end{equation}

The non-trivial relations are the normal shock relations and scalar pressure for a gas obeying a Maxwellian (ideal gas law)

(2.19)\begin{equation} \left. \begin{array}{@{}c@{}} \rho_{i0} = \dfrac{\rho_{i1}}{1 - (2/(\gamma+1))(1 - 1/M_0^{2})},\\ p_{i0} = p_{i1}\left(1 + \dfrac{2\gamma}{\gamma+1}(M^{2} - 1)\right),\\ p_\alpha = n_\alpha kT_\alpha , \end{array}\right\} \end{equation}

and the requirements of charge neutrality resulting in the relations

(2.20)\begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle n_{i0} = \dfrac{\rho_{i0}}{m_{{{}^{1}_{1}{\rm H}}}}, \quad n_{i1} = \dfrac{\rho_{i1}}{m_{{{}^{1}_{1}{\rm H}}}}, \quad n_{i2} = \dfrac{\rho_{i2}}{m_{i2}},\\ \displaystyle n_{e0} = n_{i0}Z_{{{}^{1}_{1}{\rm H}}}, \quad n_{e1} = n_{i1}Z_{{{}^{1}_{1}{\rm H}}}, \quad n_{e2} = n_{i2}Z_{{i2}} \end{array}\right\}. \end{equation}

All other parameters are set as a result of those above and the requirements of charge neutrality and mechanical equilibrium. The remaining relations are

(2.21)\begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle P_{i2} = P_{i1}, \\ \displaystyle P_{e0} = P_{i0}, \quad P_{e1} = P_{i1}, \quad P_{e2} = P_{i1},\\ \displaystyle kT_{i0} = \dfrac{P_{i0}}{n_{i0}}, \quad kT_{i1} = \dfrac{P_{i1}}{n_{i1}}, \quad kT_{i2} = \dfrac{P_{i2}}{n_{i2}}, \\ \displaystyle kT_{e0} = \dfrac{P_{e0}}{n_{i0}Z_{i1}}, \quad kT_{e1} = \dfrac{P_{e1}}{n_{i1}Z_{i1}}, \quad \textrm{and} \quad kT_{e2} =\dfrac{P_{e2}}{n_{i2}Z_{i2}}. \end{array}\right\} \end{equation}

The unique parameters for the three interface types are given in table 3.

Table 3. Distinguishing parameters for the thermal, isotope and species interface cases.

Ideally, a single shock front would be maintained across both fluids until interaction with the DI. However, the greater sound speed and shock breakdown in the electron fluid mean an electron shock wave co-incident to the ion shock cannot be produced and maintained after initialisation. The resulting general Riemann problem in the electron fluid propagates a single shock and multiple waves that traverse the interface prior to the arrival of the ion shock.

The volume of fluid tracer quantity used to track the interface is also used to calculate the mixture properties of the ion fluid. The same hyperbolic tangent transition function is used for the tracer as for the conserved properties on the interface, and the tracer has a value varying $\varrho \in [0, 1]$ from left to right states. The single ion fluid modelled in this work is technically a mixture, where the tracer property is equivalent to a mixture fraction. The tracer value is conserved and is convected as a passive scalar. A mixture equation, shown below, is used to find the particle properties when the tracer value is between zero and one. In the equation below, $\varrho$ is the tracer value, $\phi$ is some property which varies as a linear mixture of the two values

(2.22)\begin{equation} \phi_{effective} = (1-\varrho)\phi_0 + \varrho\phi_1. \end{equation}

3. Results

3.1. A 1-D exemplar

The 1-D base flow simulation results for the case of an unperturbed thermal, species and isotope interface, as per table 3, are useful for appreciating the simpler MFP mechanisms that exist in the 2-D flows. Figure 13, charge density contour plots for displacement–time axes, shows the key phenomena observable from the 1-D simulations which are the self-generation of EM field, the Lorentz force influence on the fluids, electron shock degradation and oscillatory electron motion. Immediately obvious from figure 13 is the breakdown of the general Riemann problem in the electron fluid, producing a backwards travelling rarefaction and forward travelling contact discontinuity, shock and multiple waves that are not visible in the figure. The shock and multiple waves traverse the electron-fluid density interface (EFDI), except for in the Isotope case which has no EFDI, prior to the arrival of the ion shock. The differential speed of the ion and electron shocks (and wave packets) produce a non-zero charge density that produces the first self-generated EM fields. As the electron shock processes the EFDI, the charge separation is increased as the denser side is compressed.

Figure 13. Charge-density $x{-}t$ diagrams for $d_S = 10$ scenarios. (a) Thermal discontinuity, (b) species Li-3 discontinuity, (c) isotope H-3 discontinuity and (d) species He-3 discontinuity.

Following the EFDI compression, an oscillatory behaviour emerges in the electron fluid (viewed through the charge density). This oscillation is the result of the EFDI, compressed after the shock traversal, moving through the homogeneous and mostly undisturbed ion-fluid background. As the EFDI moves to the right, it depletes the negative charge within its former location producing a net positive charge there. This deviation from charge neutrality generates an electric field that attracts the EFDI, opposing its shock imparted motion. The reciprocal of this process happens after the EFDI is accelerated to the left, causing the oscillation to the right. This oscillation only abates once the ion shock has processed the ion-fluid density interface (IFDI). The DI experiences a positive x-electric field for a long duration which is produced from a line of non-zero charge density that results from the interaction of EFDI and IFDI. This x-electric field endures until simulation end. In the two-dimensional case, this x-electric field contributes to a Lorentz force that drives the electromagnetically driven Rayleigh–Taylor instability (ERTI) in an unstable configuration, contributing to growth.

The effects discussed above occur for all interface types with differences in the extent. The thermal interface, figure 13(a), and Li-3 (fictitious element for the purpose of comparison), figure 13(b), show a very similar result with a minor difference in the interaction of the electron- and ion-fluid DI (EFDI and IFDI). After the initial shock traversal, the EFDI moves to the right with greater speed in the Li-3 species case than in the thermal case. However, the IFDI has the same speed in both, meaning that the thermal case experiences an interaction between the two fluid interfaces sooner. For the He-3 species case, figure 13(d), the severity of charge separation achieved is lesser than the thermal case and Li-3 species case scenarios because of the lower EFDI density ratio that the initiating shocks encounter. This reduction in charge density ultimately affects the magnitude of the self-generated fields in the plasma and driving of the secondary ERTI in the 2-D cases. The isotope case, figure 13(c), results highlight the significance of the initial EFDI. The cases with a lesser EFDI ratio (or none at all), figures 13(c) and 13(d), produce weaker charge separation and consequentially less energetic responses in the fluid motions. This result is especially clear in the isotope case results, where the peak fluid velocities and electromagnetic field generation is almost an order of magnitude lower than the thermal and Li-3 species cases.

Lower magnitude EM fields are generated in the isotope case, and any cases with lower EFDI ratio such as the He-3 species case, because of reduced current densities. The relative motion between the ion and electron fluids combined with different charge density across the interface, within each fluid, will drive a varying current density depending on the type of interface. The differing wave speeds in the fluids causes a significant relative motion at early times, see the small magnitude negative charge density at early time. Secondly, the oscillatory behaviour of the EFDI, being pulled back by the region of positive charge density in its wake after processing the EFDI at early times adds some complexity to the charge density. Thirdly and most significantly, different compression and acceleration of the DI in each fluid after processing by their respective shocks cause the greatest magnitude charge density and current densities. This is shown in figures 13(a), 13(b), and 13(d) by the bright red (high positive charge density) between the IFDI and EFDI after shock processing. Crucially, this region of high positive charge density is absent from the isotope case.

Examining the results for skin depths of 10, 1 and 0.1 (figure 14) demonstrates the increased coupling and subsequent reduction in multi-fluid effects for the 1-D simulations. Figure 14(a) corresponds to the charge density plot for the thermal case shown in figure 13(a). Figure 14(b) shows the reduction in skin depth has caused a reduction in the oscillation amplitudes as the electron fluid is no longer able to achieve as large separation and consequent counteracting Lorentz force that would otherwise drive the greater oscillation amplitudes. The oscillations, although reduced in amplitude, increase in frequency, figure 14(b), as the electron fluid is pulled back sooner than it would be in the loosely coupled case. Within the $d_S = 0.1$ case, figure 14(c), no oscillation is visible. The reduction in skin depth allows the electron fluid to transfer a greater portion of its accumulated energy to the ions but the electrons accumulate less energy due to the reduced differential motion from the increased coupling. In two dimensions, however, the reduced skin depth cases produce more significant MFP effects than in one dimension because (i) seeding of interface perturbations, by the electron fluid, triggers the KHI and (ii) shock refraction at oblique angles produces destabilising features.

Figure 14. Electron density distribution $x$–$t$ diagram contour plot for Thermal for varying skin depths.(a) $d_S=10$, (b) $d_S=1$ and (c) $d_S=0.1$.

3.2. Comparison of 2-D fundamental interface types

This section compares the RMI simulation results for each of the fundamental interface types under varying values of non-dimensional skin depth. The key phenomena influencing the RMI are identified for each interface type. The MFP simulations are compared with the analogous single-fluid hydrodynamic simulations. In the absence of an applied magnetic field, as we simulate, single-fluid hydrodynamic and ideal MHD simulations are identical when the same density ratio is used. We provide two reference cases for the single-fluid hydrodynamic scenarios, pressures of $p=1$ and $p=0.5$. In the case of very loose coupling, the ion and electron fluids are permitted to experience large relative motions and act as if almost independent from each other, thus the partial pressure of 0.5 in the single-fluid case is an appropriate comparison with the ion fluid. In the case of very tight coupling, the converse is true and the two fluids behave as one, thus a partial pressure of one in the hydrodynamic case is the most appropriate comparison. In this manner, the $p=1$ and $p=0.5$ hydrodynamic results create an envelope of theoretically expected behaviour for the RMI in the ion fluid. This section closes with a discussion of the relative importance of the key mechanisms.

The simulation data are presented throughout this section with the format of figures 16 and 17. The former shows the fluid mass-density distributions at simulation end for each scenario and the latter shows time series of interface statistics. The interface statistics used are the circulation ($\varGamma _z$, on a half period of the interface), generation of circulation ($\dot {\varGamma _{Z, total}}$, total) and the ion-fluid DI amplitude width ($\eta$) and growth rate ($\dot {\eta }$). The vorticity equation

(3.1)\begin{align} \frac{\partial \boldsymbol{\omega}}{\partial t} ={-}\boldsymbol{\omega}(\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{u}}) -(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}{}) \boldsymbol{\omega} + \frac{1}{\rho_\alpha^{2}}(\boldsymbol{\nabla}{\rho_\alpha}\times \boldsymbol{\nabla}{p_\alpha}) + \boldsymbol{\nabla} \times{\left(\sqrt{\frac{2}{\beta}}\frac{q_\alpha}{d_S m_\alpha} (c\boldsymbol{E}+\boldsymbol{u}_\alpha \times \boldsymbol{B})\right)},\end{align}

shows the terms which combine to give the total time rate of change of vorticity, related to circulation by

(3.2)\begin{equation} \frac{\partial \varGamma } {\partial t} = \frac{\partial } {\partial t} \iint_S \boldsymbol{\omega}\boldsymbol{\cdot} \, {{\rm d}} \boldsymbol{S} = \iint_S \frac{\partial \boldsymbol{\omega}}{\partial t} \boldsymbol{\cdot} \, {{\rm d}} \boldsymbol{S}. \end{equation}

For a small finite area in the $x$–$y$ plane we have ${\partial \varGamma _{Z}}{\partial t} = A({\partial \boldsymbol {\omega _z}}/{\partial t})$ where $A$ is the area of fluid enclosed by a loop with vorticity. If we take a cell to be the area enclosed by the enclosing loop, we can estimate the instantaneous circulation generated by multiplying terms from (3.1) by the area of a cell. Summing the contributions from the cells that comprise the interface is how the circulation generation terms are calculated. Likewise, the circulation is found from the sum of vorticity multiplied by cell area for all cells comprising the interface

(3.3)$$\begin{gather} \frac{\partial \varGamma_{Z}}{\partial t} = \sum_i A_i\frac{\partial \omega_Z}{\partial t}, \quad \text{for } i \in \text{ Interface cells}, \end{gather}$$

(3.4)$$\begin{gather}\varGamma_{Z} = \sum_i A_i \omega_{Z,i}, \quad \text{for } i \in \text{ Interface cells}. \end{gather}$$

The interface cells are determined by their volume of fluid tracer value, $\varrho$, that is convected with the flow. The interface width is calculated from the tracer average across the y-dimension where the X coordinates of $\varrho _y \in [0.05, 0.95]$ define the interface width

(3.5)\begin{equation} \varrho_y = \frac{\displaystyle \int \varrho {{\rm d} y}}{\displaystyle \int {{\rm d} y}}. \end{equation}

The maximum and minimum transition points are then used to define the interface width. Numerical derivatives are then used to find the temporal growth rate of the interface from the interface width values.

The Canny shock detector proposed by Fujimoto, Kawasaki & Kitamura (2019) was used to identify shocks travelling across the interface. The method is based off the Canny edge detection algorithm for image analysis. This detector was used rather than conventional methods because it relied on a simple principle that would be applicable despite strong EM influences. The algorithm steps through (i) computing pressure gradient magnitudes and direction, (ii) collection of candidates – local maximum search and filter, (iii) and assessment of candidate points (comparing with shock jump conditions with an error tolerance). The Canny shock detector produced reasonable results that are used in the analysis for tracking the traversal of transverse shocks across the interface.

The interface types represented on the time series plots are denoted as follows: H-p = 1 – single-fluid hydrodynamic reference solution with pressure of 1, H-p = 0.5 – single-fluid hydrodynamic reference solution with pressure of 0.5, I-$d_S$-10 – IRMI with skin depth of 10, S-$d_S$-10 – hydrogen-helium SRMI with skin depth of 10, L-$d_S$-10 – hydrogen-lithium SRMI with skin depth of 10 and T-$d_S$-10 – TRMI with skin depth of 10. Figure 18 shows the density contour plot for the reference single-fluid solutions as discussed in the preceding paragraph.

When interpreting the results of the 2-D cases terminology introduced in figure 12 is used: (i) whenever the y-dimension is discussed in context of the DI, the lower x-symmetry plane is considered, (ii) the spike refers to a region of high density fluid penetrating lower density, (iii) the opposite case to point 2 is referred to as a bubble, (iv) the region transitioning between the spike and bubble is referred to as the neck and (v) the bottom left numbers within brackets on contour plots are the minimum and maximum values for that figure.

It is useful to consider which regime the plasma RMI is in through the analogous hydrodynamic problem. The $p=1$ and $p=0.5$ hydrodynamic results for width and growth rate are compared with Richtmyer's impulsive model (Richtmyer 1960) in figure 15. Considering also figure 18, we expect the hydrodynamic cases and likely the plasma RMI will be transitioning between the linear and nonlinear growth phases.

Figure 15. Hydrodynamic simulation results for (a) interface width $\eta$ and (b) growth rate $\dot {\eta }$, for the hydrodynamic cases. Richtmyer's impulsive model (Richtmyer 1960) is provided for comparison.

3.2.1. Loosely coupled fluids $d_S=10$

A large non-dimensional skin depth indicates a loose coupling between the electron and ion fluids that ultimately enhances secondary instabilities. The EM forces couple the two fluids together but at this skin depth significant relative motion between them is permitted. This relative motion is critical in the development of secondary instabilities because it allows the evolution of expansive EM fields and allows the electron to perturb the ion interface. Figure 16 shows the ion and electron fluids at simulation end along with the Lorentz force in the x-direction. Charged particles within EM fields experience a Lorentz forces that provides a continuous forcing to drive the RTI. Characteristic RTI features (mushrooms) can be seen on the ion-fluid bubbles of the SRMI-Li case in figure 16. EM fields can similarly generate torques due to the gradient of Lorentz forces which acts as a source of vorticity that can amplify or diminish the KHI depending on the sign or the vorticity. The KHI rollers are seen along the neck of the SRMI-He, SRMI-Li and TRMI cases, where it is more severe in the latter two due to the greater EM forces and consequent torques. The relative motion also produces interactions between the fluids that introduce additional interface perturbations, exacerbating the primary RMI and secondary instabilities. The secondary instabilities augment the primary interface perturbation growth rates, increasing them, as well as contributing to the emergence of high-wavenumber secondary interface perturbations.

Figure 16. Electron- and ion-fluid mass density and x-Lorentz force contours at simulation end time with Debye length ($\lambda _D = 2\times 10^{-1}$). Bracketed tuples indicate the minimum and maximum values on each contour.

Problem symmetry is obviously broken in the electron fluid and some ion-fluid features also indicate this. The root cause is the interaction of complex systems of many vortices and inhomogeneous machine precision error. The small numerical deviations across the domain due to rounding error are not homogeneous, these differences are then amplified by the complex vortex field.

The reader will note that the isotope interface is distinctly restrained despite the permitted relative motion between the fluids. This relative motion is important for generation of EM fields that drive the secondary ERTI and can augment the evolution of the secondary KHI instability. The Secondary instabilities occurring in two dimensions introduce translational motion in directions other than the shock propagation direction as well as rotational motion. These effects on relative motion are in addition to the influences discussed for one dimension i.e. oscillatory motion, differing shock speeds and fluid charge-density discontinuities across the interface, see § 3.1. Furthermore, the development time scales of the secondary instabilities in electron and ion fluids are more disparate than the shock time scales, therefore relative motion with the instabilities is greater than without. Since the electron time scales for instabilities are much smaller than ions, the presence of a significant (or indeed non-zero) initial electron-fluid density ratio produces much greater relative motion than when only the ion fluid possesses a density ratio.

Figure 16 shows a dramatic difference in both evolution and severity of the IRMI compared with the other interface types. The primary factor in this difference is the absence of an initial electron-fluid interface in the IRMI. This absence results in greatly reduced EM field generation and effective removal of ion interface excitation by the electron fluid prior to the ion shock arrival. The effect of the ERTI is subsequently reduced due to its driving Lorentz force being proportional to the EM fields. The KHI is also reduced due to decreased shear and fewer seed points along the interface which would have been excited by energetic electron-fluid perturbations. All these changes produce a slower amplitude growth rate that aligns closely with the single-fluid reference case (green dashed line labelled H-p = 1), as seen in figure 17.

Figure 17. The evolution of the IFDI (a) cumulative circulation $\varGamma _Z$, (b) total instantaneous generation of circulation $\dot {\varGamma _{Z, total}}$, (c) interface width $\eta$ and (d) growth rate $\dot {\eta }$, for all scenarios with $d_S = 10$ ($\lambda _D = 2\times 10^{-1}$).

Figure 18. Reference single-fluid hydrodynamic solution for the RMI showing the simulation end time ion density distribution. Bracketed tuples indicate the minimum and maximum each contour.

An interesting aspect of the isotope case is the bubble that forms on the x-symmetry line in the electron-fluid DI, where all other interfaces produce a spike. The bubble is more clearly seen in figure 19 in the $\rho _e$ contour plot at $t=1$. Initially, in the electron fluid there exists only a shock and a homogeneous and stationary density field. As the simulation begins, the breakdown of the electron shock and the differing propagation speed of the ion shock and the remaining electron waves produces a narrow-planar band of higher density in the electrons, figure 19 up to $t=0.005$. The acceleration of this narrow band of electrons at early times is due to the x-electric field. This acceleration is aligned with the band's density gradient and therefore does not produce any baroclinic torque. However, after the ion shock processes the perturbed ion-fluid DI, figure 19 $t=0.165$, appreciable y-electric fields are produced, consequently misaligning the directions of the pressure and density gradients, generating baroclinic torque. The electric field, which dominates the Lorentz force in these simulations, is opposite signed either side of the x-symmetry line, producing opposite signed baroclinic torque either side of the x-symmetry line. The result is positive and negative baroclinic torque above and below the x-symmetry line, acting to contort the narrow band of electrons into a bubble. By $t=0.4$, after the initial contortion of the narrow band by the y-electric field, the x-electric field is now misaligned with the density gradient along the interface. The long lived x-electric field continues to drive continual baroclinic torque generation that reinforces the bubble. The bubble continues to grow until simulation end since there is no significant counteracting torque.

Figure 19. Evolution of the IRMI (a) electron density contour, (b) x- and (c) y- electric fields, left to right, for $\lambda _D = 2\times 10^{-1}$. Bracketed tuples indicate the minimum and maximum each contour.

The electron-fluid flow field that evolves in the IRMI is not produced in the other interface types and may provide a unique seed for instabilities at later times. The formation of a bubble or spike on the x-symmetry line depends on the relative magnitude and orientation of the circulation produced by the electron shock passage and the ERTI (note the baroclinic torque is one source of circulation). In the TRMI and SRMI, the circulation deposited by the ERTI (at early simulation time) is lesser in magnitude and opposes the electron shock deposited circulation, and so it has no effect on the direction of the spike. However, in the IRMI the ERTI circulation does not encounter an appreciable opposing circulation in the electron fluid and so it generates a bubble where all other interface types form a spike. The IRMI bubble then grows and develops features along its extent which will evolve the opposite manner to the analogous cases in the TRMI and SRMI for the same continuous or impulsive forcing due to the direction of the density gradient. This distinctly opposite evolution may lead to unique and fundamentally different behaviours at late time.

Figure 16 shows the Li-3 SRMI final state is closer to the TRMI results of Bond et al. (2017b) than the He-3 SRMI case. The Li-3 species case exhibits three distinct qualities in the ion fluid when compared with the He-3 species case; a more severe RMI in terms of amplitude growth rate and final amplitude; a much more perturbed interface exhibiting more developed features (KHI rollers and ERTI spikes); and a narrower spike and neck region. To reiterate, the electron-fluid density ratio greatly affects the overall plasma RMI since it governs the charge separation induced by the relative motion of the interfaces. The increase in electron density ratio in the Li-3 case to the levels of the thermal interface results of Bond et al. (2017b) and subsequent aligning of RMI amplitude and growth rate suggests that the electron density ratio may play a dominant role in the plasma RMI development at large Debye lengths. This is a physically plausible mechanic when understood from the perspective of the secondary ERTI. The ERTI is driven by EM fields to which the electron-fluid dynamics contributes with electron motion and separation from the ion fluid. Both of these factors are amplified in the case of a greater electron-fluid density ratio. Likewise, the increased interface disturbance is due to the more energetic electron fluid exciting the ion fluid. However, the third distinct feature – the narrowing of the spike in comparison with the thermal case – must be a consequence of varying ion properties of temperature, number density and charge between the initial parameters of the TRMI and Li-3 SRMI.

The thermal and lithium species cases are very similar in their evolution owing to their identical electron-fluid and ion-fluid density ratios. The TRMI does, however, produce a greater amplitude growth during simulation time, figure 17 $\eta$ time series. Figure 16 shows the SRMI-Li has a narrower neck than the TRMI at final time, the neck width of the TRMI (0.52) is 30 % greater than the SRMI-Li (0.4) at their narrowest points. Figures 20 and 21 show the symmetry of the SRMI-Li case is coherent whereas the TRMI case has broken symmetry at large scale. The narrowing effect appears to be influenced by Lorentz forces and transverse shocks. The broken symmetry in the TRMI affects shock refraction and EM field evolution, thereby influencing the shock strength and direction, and Lorentz forces driving the fluid flow. In regard to a difference in shock refraction, the SRMI-Li (figure 21), shows two coherent shocks having just traversed the interface (from both y-boundaries) and moving towards the line of symmetry, however, the TRMI (figure 20) shows only one shock. Additionally, the y-Lorentz force acting along the spike is positive (negative) when above (below) the axis of symmetry, leading up to and after the narrowing of the spike working to widen the spike. Beyond the Lorentz force, the only linear drivers of motion within the fluid are the existing pressure gradients. Inspection of the magnitudes of Lorentz forces show that the Li-3 SRMI experiences weaker y-Lorentz forces, which act to widen the spike, than the TRMI. Furthermore, there is a lower pressure behind the ion-fluid DI in the Li-3 SRMI and greater pressure gradients, which mostly act to narrow the spike. These factors result in the much more pronounced narrowing in the Li-3 case, from the spike to the neck, as opposed to only the spike in the TRMI.

Figure 20. The TRMI (a) Canny shock detector results, (b) low filtered pressure gradient magnitude,(c) pressure, (d) x-Lorentz force, (e) y-Lorentz force for $\lambda _D = 2\times 10^{-1}$ at $t-0.9$.

Figure 21. The SRMI-Li (a) Canny shock detector results, (b) low filtered pressure gradient magnitude,(c) pressure, (d) x-Lorentz force, (e) y-Lorentz force for $\lambda _D = 2\times 10^{-1}$ at $t-0.9$.

The instability evolution will change significantly when moving from the present 2-D 3-vector (2D3V), i.e. two spatial dimensions but all three vector components, work to a full 3D3V simulation. The coherence observed in the 2-D case e.g. KHI rollers will be subject to 3-D instabilities and lose their 2-D coherent structures. This type of interaction is likely for the primary instability spikes and bubbles also. In addition, turbulent mixing zones will likely evolve in the longer time scales. These effects will likely culminate in reduced instability growth.

It would be practical, using the same effective resolution, to produce a 3-D simulation of the ideal MFP RMI with Cerberus. It is estimated that the 3-D simulation would require 33 M cells, increasing from approximately 2 M cells for the current work. This is only possible due to the AMR framework provided by AMReX. However, simulations with transport processes (collisions) seem to be out of reach for 3-D simulations due to the severe time step constraint imposed by diffusive and conductive time scales, and the current configuration of the code (explicit and only second order accurate in time). Investigation of multi-mode perturbations are possible and could be easily implemented through designing the initial density interface geometry accordingly.

3.2.2. Moderately coupled fluids $d_S=1$

The results of the reduction in skin depth, figures 22 and 23, show the SRMI-He was the only scenario to produce an increase in interface amplitude when compared with the loosely coupled fluid case, figures 16 and 17. The loosely and moderately coupled cases develop a similar interface circulation deposition profile until $t\approx 0.7$ where the moderately coupled case ($\lambda _D=2\times 10^{-2}$) sees an increased deposition of circulation on the interface, where the loosely coupled case drops, compare figures 17 and 23. The late time rally in circulation produces a final interface amplitude width that is greater than the preceding Debye length case, 0.412 compared with 0.402 non-dimensional lengths. The combination of increased amplitude width and development of high-wavenumber interface perturbations make the $\lambda _D=2\times 10^{-2}$ SRMI-He case more severe than the loosely coupled case in terms of the detrimental effects for ICF experiments.

Figure 22. Electron- and ion-fluid mass-density contours at simulation end time $d_S = 1$ ($\lambda _D = 2\times 10^{-2}$). Bracketed tuples indicate the minimum and maximum each contour.

Figure 23. The evolution of the IFDI (a) cumulative circulation $\varGamma _Z$, (b) total instantaneous generation of circulation $\dot {\varGamma _{Z, total}}$, (c) interface width $\eta$ and (d) growth rate $\dot {\eta }$, for all scenarios with $d_S = 1$ ($\lambda _D = 2\times 10^{-2}$).

The IRMI in the ion fluid produces a reduced interface amplitude, when compared with the loosely coupled case, but shares the same fluid and EM features. The electron fluid, however, is different, producing a greater number of ERTI features on the electron-fluid DI. Figure 24 shows the inverted primary electron-fluid bubble centred on the x-symmetry line is still present (observed in the larger Debye length case for the IRMI) although diminished in size. Additionally, there is a secondary spike and bubble adjacent to the primary bubble on either side of the x-axis of symmetry. The secondary spike and bubble are caused by the same mechanism as discussed in § 3.2.1 but with the smaller and alternating EM field regions seeding the alternating pattern of bubble and spike that is amplified by the imprint of the ion fluid transverse-reflected shocks (TRS) on the electron fluid. These ion-fluid TRS imprints are not significant in the $\lambda _D=2\times 10^{-1}$ IRMI case because of the weaker ion–electron coupling and the weak strength of the TRS. While the $\lambda _D=2\times 10^{-2}$ case has some distinctions from the $\lambda _D=2\times 10^{-1}$ case, these do not appear to be caused by a unique mechanism. In fact, had there been a greater density gradient in the electron fluid it is likely the inversions along the imprinted DI would not occur because its source, the TRS, does not seem strong enough in the cases studied thus far to produce enough circulation to overpower the initial deposition of circulation by incident ion and electron shocks.

Figure 24. The IRMI $\lambda _D=2\times 10^{-2}$ ion and electron density contours at simulation end showing the ERTI and transverse-reflected shock generated perturbations. Bracketed tuples indicate the minimum and maximum each contour.

A reduction in skin depth shrinks the spatial scale of EM features and increases their intensity, overall reducing the influence on EM driven secondary instabilities. The ion x-Lorentz forces (directly proportional to the EM fields) shown in figure 20 and figure 25 are indicative of the EM fields for the $\lambda _D=2\times 10^{-1}$ and $\lambda _D=2\times 10^{-2}$ cases, respectively. In the $\lambda _D=2\times 10^{-1}$ case (§ 3.2.1), the ion and electron fluids behave as if they were almost independent of each other, with their large separations giving rise to coherent and expansive regions of EM fields and forces. These EM forces compound the effects of the RMI to produce greater growth rates and amplitudes. In the $\lambda _D=2\times 10^{-2}$ case, the self-generated EM fields exist in much smaller regions of alternating and intense strength. The net effect on interface growth with fine-grained regions of alternating EM fields and consequent Lorentz forces is lesser than the case of larger alternating regions of EM fields, even though the magnitudes of the Lorentz forces are greater. The ERTI for the moderately coupled fluids is thus reduced when compared with the tightly coupled case but is still a significant contributor to the interface growth rate.

Figure 25. An example of the intense and alternating regions of Lorentz force causing the erratic interface statistics for a smoothly evolving interface area; (a) $\lambda _D=2\times 10^{-2}$ He-3 SRMI x-Lorentz force at $t=1$ with interface contour whited out and (b) interface area for all interface and skin depth scenarios.

Increased coupling allows the electron fluid to excite the ion fluid interface with high-wavenumber disturbances that are then amplified by the shock refraction and secondary instabilities. The electron fluid transfers (relative to the loosely couple case) more of its energy to the ion fluid, imparting high-wavenumber content along the interface that is then amplified through deposition of circulation by the secondary ERTI and KHI. The high-wavenumber content is further amplified by shocks created from a refraction process not observed in the loosely coupled case. The numerous refracted shocks are observable in figure 22 behind the transmitted shock. These shocks result from the refraction of the incident ion-shock with the electron-excited ion-fluid DI and the formation of triple points along the transmitted ion shocks. The interaction of the transmitted shocks are clearly observable in the cellular patterns downstream of the ion-fluid DI. The interface is multiply shocked, initially from the incident shock and then the shocks emanating from the triple point sweeping across the ion-fluid DI.

The source of interface circulation deposition and growth rate appear to be highly transient phenomena. This is suggested by the time rate of change of circulation and interface amplitude in figure 23. The 2-D circulation contributions are given by the vorticity multiplied by relevant fluid area and the vorticity itself is governed by the equation

(3.6)\begin{align} \frac{\partial \boldsymbol{\omega}}{\partial t} ={-}\boldsymbol{\omega}(\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{u}}) -(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}{}) \boldsymbol{\omega} + \frac{1}{\rho_\alpha^{2}}(\boldsymbol{\nabla}{\rho_\alpha}\times \boldsymbol{\nabla}{p_\alpha}) + \boldsymbol{\nabla} \times{\left(\sqrt{\frac{2}{\beta}}\frac{q_\alpha}{d_S m_\alpha} \left(c\boldsymbol{E}+\boldsymbol{u}_\alpha \times \boldsymbol{B}\right)\right)}.\end{align}

The sources of vorticity in the ideal case, ignoring stretching of vorticity due to velocity gradients and flow compressibility, are the EM fields and fluid pressure gradients (imposed by the initial shock and subsequent refracted shocks). Therefore the sources of vorticity and/or the interface area are highly transient. Figure 25(b) shows the time series of the interface area – the interface area is defined as the sum of each cell's area where the cell contains between 0.05 and 0.95 volume of fluid tracer and is used to calculate the interface statistics. The simulation produces an interface area that varies smoothly with time, therefore, an incorrect identification of the interface is not the cause of the observed volatility in circulation. Rather, the cause is the highly energetic and variable nature of the $\lambda _D=2\times 10^{-2}$ scenarios with regard to the shock refraction and EM field development. Figure 25(a) shows the highly variable x-Lorentz force acting on the ion fluid for the SRMI-He scenario, contributing to the variability in circulation generation on the interface. Figure 22 shows the numerous shocks generated from triple points, after refraction through the interface, which then interact with the interface once more. These numerous and incoherently propagating shocks provide a variable pressure and density along the interface that contributes to the highly variable generation rates of circulation in figure 23.

3.2.3. Tightly coupled fluids $d_S=0.1$

The smallest Debye length cases investigated appear to be approaching the single-fluid limit. Figure 26 shows a high degree of resemblance between the ion and electron density contours in the TRMI and both SRMI cases. The IRMI also shows an imprint of the IFDI on the electron fluid. Indeed, the two fluids appear to be acting almost as one in all cases. The $\lambda _D=2\times 10^{-3}$ produces TRMI and SRMI growth rate and interface amplitudes significantly below the reference single-fluid case, a significant difference from preceding $\lambda _D$ cases. The EM effects also appear to be reducing in magnitude across the interface.

Figure 26. Electron- and ion-fluid mass-density contours at simulation end time $d_S = 0.1$ ($\lambda _D = 2\times 10^{-3}$). Bracketed tuples indicate the minimum and maximum each contour. Note the contour plot scale is saturated due to the strong Lorentz forces surrounding the original transmitted shock front.

No new phenomena occur in the IRMI evolution in the case of tightly coupled fluids. The only modification, which is a common element of all interface cases, is what seems to be an increase in the refraction frequency of the TRS. As mentioned above, the shock refraction process with a small Debye length produces many reflected waves at high frequency, in terms of emission from the interface rather than wave frequency, and in close spacing to each other. The combination of the ERTI and the imprint of the many closely spaced ion TRS on the electron fluid seems to produce many more perturbations that grow along the ion-fluid DI spike and neck. Interestingly, the relative growth of high-wavenumber perturbations decreases from large to small Debye lengths, suggesting a decrease in the effectiveness of the ERTI and KHI. The decrease in secondary instabilities supports the idea that reducing differential movement between the fluids reduces the expanse of the EM fields, exposing the interface to many alternating regions of influence (consequent vorticity either reinforcing and opposing instabilities) instead of fewer regions of a definitive action. Overall the IRMI does not produce a large difference in ion fluid DI amplitude or growth rate across the Debye length trials, as seen from a comparison of figures 17, 23 and 27.

Figure 27. The evolution of the IFDI (a) cumulative circulation $\varGamma _Z$, (b) total instantaneous generation of circulation $\dot {\varGamma _{Z, total}}$, (c) interface width $\eta$ and (d) growth rate $\dot {\eta }$, for all scenarios with $d_S = 0.1$ ($\lambda _D =2\times 10^{-3}$).

The $\lambda _D=2\times 10^{-3}$ TRMI and both SRMI cases produce a significant MFP shock refraction effect despite approaching the single-fluid limit (p = 1 reference case in figure 18). The shock refraction of the incident shock on the DI produces a reflected shock, transmitted shock, and high frequency waves. The high frequency waves, observed in figure 28, are generated on the inclined portions of the interface, above and below the x-symmetry line, and propagate along the interface. The waves show strong structural resemblance in the pressure and density gradient contour plots. As the high frequency waves pass along the interface they seed high-wavenumber instabilities. These waves then pass the interface and interact with the contact discontinuity – generated from the spontaneous breakdown of the initial electron shock of which an imprint is made on the ion fluid due to high coupling of the fluids – and excite the contact discontinuity. The contact discontinuity later impacts the interface, conforming around it but not penetrating during simulation time. The new interface perturbation development mechanism appears much less severe in size and number of features than the mechanism observed in the larger Debye length cases. However, the perturbations left behind are of a higher frequency. Note this interface perturbation development mechanism is different to that described for the moderately coupled case ($\lambda _D=2\times 10^{-2}$) as that was driven by shocks formed from triple points after the initial refraction. Here in the $\lambda _D=2\times 10^{-3}$ case, waves are generated directly from the refraction of the incident ion shock. The absence of this mechanism in the IRMI, see figures 28 and 29, and comparison with the single-fluid results suggests this is a phenomena unique to the two fluid model.

Figure 28. The SRMI-He (a) shock detector, (b) low filtered pressure gradients and (c) x- and (d) y-Lorentz force components for $\lambda _D = 2\times 10^{-3}$.

Figure 29. The IRMI (a) shock detector, (b) low filtered pressure gradients and (c) x- and (d) y-Lorentz force components for $\lambda _D = 2\times 10^{-3}$.

3.2.4. Qualitative hydrodynamic RMI features in the plasma RMI

The plasma RMI is subject to significant electromagnetic effects not present in the hydrodynamic case. Despite the additional Lorentz forcing terms, that provide additional sources of bulk fluid motion and direct circulation contributions, the typical qualitative features remain. We have observed already the plasma RMI instability evolves similarly to the hydrodynamic case, exhibiting the secondary KHI, but with the addition of Lorentz forcing and circulation generation, and the secondary ERTI. In the smallest skin depth case, an interesting parallel between the plasma and hydrodynamic RMI is that between the electromagnetic effects and the background motion. In the hydrodynamic RMI the fluid motion is most intense near to the interface and small away from it. A similar observation is made for the electromagnetic activity which is most intense near to the interface and negligible away from it (between the reflected and transmitted shocks), see figure 26 x-Lorentz force. The behaviour of the background motion, in the small skin depth plasma RMI, follows that predicted by the hydrodynamic theory, figure 30.

Figure 30. The ion-fluid x-velocity field in the stationary reference frame (a) and the background motion reference frame (b) for the TRMI small skin depth, $d_S = 0.1$.

The large skin depth behaviour of the plasma RMI diverges more from the hydrodynamic case due to the interaction of the electron and ion fluids which have loose electromagnetic coupling. The effect is akin to two miscible and separate RMI evolutions superposed with each other. The loose electromagnetic coupling allows expansive regions of electromagnetic fields to develop that affect both fluids. There is significant electromagnetic activity away from the interface, see the x-Lorentz force in figure 16, and there is significant background motion away from the interface, figure 31. In the SRMI-Li large skin depth case, shown in figure 16, we see a fluid evolution that resembles bubble flattening that is much more pronounced than for the other interface types, raising doubts that the mechanism for this is the same in the hydrodynamic case (Wouchuk 2001; Abarzhi et al. 2003; Pandian et al. 2017).

Figure 31. The ion-fluid x-velocity field in the stationary reference frame (a) and the background motion reference frame (b) for the TRMI large skin depth, $d_S = 10$, case.

The apparent bubble flattening in the SRMI-Li case is a result of the electromagnetic effects rather than the conventional RMI nonlinear effect. We see from the hydrodynamic reference solutions (the $p=0.5$ being relevant due to the loose coupling), figure 18, and the isotope case, figure 16, that the RMI is exhibiting early signs of bubble flattening. One might expect that the additional circulation generation from electromagnetic terms, in the plasma RMI with appreciable electron fluid-density interface, would drive the RMI to more quickly develop nonlinear effects such as bubble flattening. However, the Lorentz forces within the plasma resist changes to the charge distribution (such as that from a changing ion DI) thereby slowing the RMI evolution as seen in the TRMI and SRMI-He cases. The SRMI-Li case behaves differently, appearing to proceed with an evolution that resembles bubble flattening. The difference (consider a comparison of the TRMI and SRMI-Li cases) is the location of strongest vortex cores in the electron fluid, indicated by the low density spots in the mass-density contours of figures 32 (TRMI) and 33 (SRMI-Li). The strongest vortex core in the SRMI-Li and TRMI cases are closely behind the spike and far behind the spike, respectively. The location of these strongest vortices produces a significantly different electron-fluid evolution and consequently a different electromagnetic field evolution and subsequent resisting Lorentz forces. Comparison of figures 32 and 33 (which have the same contour levels) at times before the bubble flattening illustrates these different evolutions. At $t=0.68$ the SRMI-Li case has large regions of ion y-Lorenz forcing toward the x-symmetry line (reinforcing flattening) in the bulk fluid adjacent to the bubble, whereas the TRMI has fewer of these regions and more regions of ion y-Lorenz forcing away from the x-symmetry line. At $t=0.8$ as the bubble flattening is starting to become significant in the SRMI-Li case, the TRMI case has comparatively little change. At this point the SRMI-Li case shows weaker resisting y-Lorentz forces along the DI in comparison with the TRMI case. These Lorentz forces in the SRMI-Li facilitate the flattening and diverge from the TRMI and SRMI-He cases. At later times, in response to the narrowing in the SRMI-Li case the Lorentz forces intensify to resist further narrowing of the region between the spike and bubble.

Figure 32. The electron mass density, charge density and ion y-Lorentz force for the TRMI large skin depth, $d_S = 10$, case. The ion interface contour is shown in purple.

Figure 33. The electron mass density, charge density and ion y-Lorentz force for the SRMI-Li large skin depth, $d_S = 10$, case. The ion interface contour is shown in purple.

There is potential for future work to investigate more closely the relation between the plasma RMI and the HRMI theory, particularly around the non-monotonic relationship between growth rate and initial amplitude. Direct quantitative comparison between simulation results of the plasma RMI to hydrodynamic theory would not be meaningful as the influence of electromagnetic effects has been shown to be substantial. A useful but challenging theoretical analysis would incorporate equations that can account for the coupling effect between the electron and ion fluids, for finite skin depth, that lead to variable generation and evolution of electromagnetic fields that drive the secondary ERTI. A possible approach is to model simplified (incompressible) sets of equation for each species, extract charge-density and current-density equations, and relate these to the electromagnetic fields, but this is beyond the scope of the present work.

4. Discussion of dominant parameters

Each interface type produces a unique evolution of the RMI that manifests as different interface amplitudes, growth rates and imposed interface perturbations. These differences raise concern for the simplifications made in single-fluid models for specific plasma regimes. The key parameters influencing the evolution of the RMI under different interface types are those which significantly affect the severity of the RMI and secondary instabilities in regard to primary perturbation amplitude and growth. These key parameters are (i) the density ratio across the initial electron-fluid DI and (ii) the plasma Debye length. The ion charge may have a unique effect on the evolution of the RMI but this is not confirmed. The density ratio across the electron-fluid DI and Debye length greatly influence the manifestation of the ERTI, TRS, and shock refraction process. All of these phenomena are uniquely MFP effects and are absent in the single-fluid case.

The influence of the EFDI density ratio is particularly observed through comparison of the IRMI and other scenarios. The IRMI simulations demonstrate a negligible interaction of the initial shock and the EFDI. Consequently the well-behaved electron-fluid motion does not generate significant charge separation and EM fields, hence it does not convert significant amounts of energy to the IFDI via the EM fields. The electron fluid, in the case of low or zero EFDI density ratio, is unable to influence the evolution of the RMI and secondary instabilities. The electron fluid is unable to (i) perturb the IFDI prior to ion shock arrival and (ii) generate significant EM fields that would drive Lorentz force oscillations and the secondary ERTI in the ion fluid. Without an initial EFDI, the single-fluid and IRMI cases, a distinctly single-fluid RMI evolution is observed. Conversely, the cases with significant EFDI density ratio, TRMI and SRMI cases, produce an energetic and erratic electron fluid, greater magnitude EM fields (EM energy density as well), and MFP effects that augment the ICF RMI.

The expectation of approaching the magnetohydrodynamic limit with vanishingly small Debye length has been shown comprehensively in the literature (Bellan 2008; Chen 2016; Shen et al. 2018) and is supported by the simulation results. Generally, decreasing only the Debye length via skin depth results in significantly reduced amplitude and growth rates. The one exception is the SRMI-He $d_S=1$ where more expansive EM fields, relative to the $d_S=1$ SRMI-LI and TRMI, lead to a greater reinforcing of the RMI by the ERTI. However, reducing the Debye length further aligns the SRMI-He case with the single-fluid amplitude. There are many plasma regimes, both within ICF and other applications, that do contain a finite Debye length which requires the use of MFP modelling to capture relevant phenomena.

Interestingly, all scenarios with an initial EFDI experience a maximum in interface perturbations at $d_S=1$. This suggests that excitation of the IFDI interface by the electron fluid relies on competing mechanisms. It seems that the scenario of $d_S=1$ may be close to a local maximum of multi-fluid effects, where the coupling of the ion and electron fluids is tight enough to allow significant excitation but the coupling is also flexible enough to allow for the necessary charge separation to produce strong enough fields to pump the electrons with energy to excite the ions. This effect is, however, diminished in the case of the $d_S=0.1$ scenarios where the excitation of the ion fluid by the electron fluid is visible in the high-wavenumber perturbations of the interface, however, the ERTI is necessarily diminished by the reduced strength and expanse of EM fields.

A comparison of the SRMI-Li and TRMI cases suggests the ion charge has a non-trivial effect on the interface evolution. Both cases have the same EFDI density ratio and despite having differing ion properties for S2 most fluid properties are identical. The only differing bulk property is the ion temperature – 0.5 for SRMI-Li and 0.167 for the TRMI – but shock propagation speeds are identical in the fluids due to the ratio of molar masses. The differing ion charge effect manifests as increased narrowing of the neck region and increased number of secondary perturbations along the bubble region of the ion interface in the loosely coupled case ($\lambda _D=2\times 10^{-1}$). Some differences also manifest in the smaller Debye length cases but to a lesser extent. The SRMI-Li and TRMI show different shock refraction processes, which would explain the different effects seen as each can be attributed to the TRS interacting with the IFDI.

5. Conclusions

The MFP model produces a different RMI evolution for the isotope, species, and thermal DI cases that are indistinguishable in the ideal MHD model and other single-fluid models. MFP models are able to distinguish between ions, isotopes, and electrons allowing the modelling of charge separation and physical evolution of EM fields that is not possible with single-fluid models. The uniquely MFP effects that result affect the evolution and severity of the RMI in the interface cases studied which would otherwise be missed in the application of single-fluid models such as the ideal MHD model. The key phenomena influencing the interface scenarios are the

1. (i) primary RMI;

2. (ii) electromagnetically driven RTI;

3. (iii) Kelvin–Helmholtz instability;

4. (iv) electron-fluid excitation of the ion-fluid interface;

5. (v) Lorentz force bulk fluid accelerations and vorticity deposition;

6. (vi) transverse-reflected shock-wave interaction with the ion-fluid interface;

7. (vii) a MFP shock refraction process.

Points (i), (iii) and (vi) are observable in the hydrodynamic models of the RMI and MHD reductions of the plasma RMI, however, the remaining phenomena are beyond their grasp because of the required capacity to model self-generated EM fields and the evolution of electron fluid. All three interface types studied (isotope, species and thermal RMI) experience the phenomena above to varying degrees. The isotope case, lacking an initial DI in the electrons is the closest to the single-fluid limit. The TRMI and SRMI scenarios with their initial electron-fluid mass-density ratios experience significant multi-fluid effects and amplification of the RMI in the larger Debye length cases. As the Debye length is reduced (increasing coupling between the ion and electron fluids ) all three novel interface types tested experience a significantly reduced RMI growth rate and amplitude, approaching the single-fluid limit but retain multi-fluid phenomena. The results of the study show that MFP effects should be considered important when (i) Debye lengths are large, and (ii) distinct DIs are formed from isotope, species and thermal discontinuities.

In relation to ICF experiments, the phenomena discussed here present several points that pose a risk of degrading capsule implosion. As seen in the results, the species and thermal interfaces experience the worst instability growth (in comparison with hydrogen isotopes). A high atomic number species, such as carbon in the ICF fuel casing material, that is ionised during the implosion will produce a large density ratio in the electrons (when forming an interface with the ionised fuel). This could produce all the detrimental effects of a highly chaotic and energetic electron-fluid RMI discussed in this paper. The shock refraction process at small Debye lengths results in significant interface perturbation. This increased surface perturbation is a concern for maintaining the carefully designed implosion geometry and coherence during implosions. The results represent only the early time evolution of the instability and under ideal conditions, introduction of transport processes is a vital next step in assessing the effect on ICF experiments.