Alfvén number for the Richtmyer-Meshkov instability in magnetized plasmas

Figure 1(a) illustrates the initial configuration for the single-mode analysis of the RMI when a rarefaction wave is reflected. Two fluids with different densities, $\rho_{1}$ and $\rho_{2}(<\rho_{1})$, are separated by a contact discontinuity at $y=Y_{\rm cd}$. A planar shock propagating through the heavier fluid “1” ($y>Y_{\rm cd}$) hits the interface at a time $t=0$. Here the $x$- and $y$-axis are perpendicular and parallel to the shock surface. The incident shock velocity is $U_{i}(<0)$ and the fluid velocity behind the shock is $V_{1}(<0)$. The pre-shocked pressure $P_{0}$ is uniform in both fluids. The Mach number of the incident shock is defined as $M=|U_{i}|/c_{s1}$ where $c_{s1}=(\gamma P_{0}/\rho_{1})^{1/2}$ is the sound speed of the fluid “1”. The interface has an initial corrugation of a sinusoidal form, $y=Y_{\rm cd}+\psi_{0}\cos(kx)$, where $\psi_{0}$ is a corrugation amplitude, $k=2\pi/\lambda$ is the perturbation wavenumber, and $\lambda$ is the wavelength.

The initial geometry of a magnetic field is assumed to be uniform with the size of $|\bm{B}|=B_{0}$ in the pre-shocked region ($y<Y_{\rm is}$). As for the field direction, we consider the cases of a parallel MHD shock ($B_{y}=B_{0}$) and a perpendicular shock ($B_{x}=B_{0}$). A uniform magnetic field parallel to the incident shock velocity is assumed in the entire region for the parallel shock case. While, for the perpendicular shock case, the initial field is calculated from the following vector potential to avoid the presence of the normal field at the interface,

using the relation $\bm{B}=\bm{\nabla}\times\bm{A}$. The physical quantities in the post-shocked region behind the incident shock are calculated from the Rankine-Hugoniot conditions for MHD shocks.

Only four non-dimensional parameters characterize the initial configuration depicted by Figure 1(a). The sonic Mach number $M$ parameterizes the incident shock velocity. The contact discontinuity is expressed by the density jump $\rho_{2}/\rho_{1}$ and the ratio of the corrugation amplitude to the wavelength $\psi_{0}/\lambda$. The initial field strength is given by the plasma beta $\beta_{0}=8\pi P_{0}/B_{0}^{2}$, which is the ratio between the gas and magnetic pressures defined at the pre-shocked region. Various situations can be examined only by choosing these four parameters; $M$, $\rho_{2}/\rho_{1}$, $\psi_{0}/\lambda$, and $\beta_{0}$.

The basic equations of ideal MHD are solved to study the interaction between the RMI and a magnetic field;

where $\rho$, $\bm{v}$, and $\bm{B}$ are the mass density, velocity, and magnetic field, respectively, and $e$ is the total energy density,

In our simulations, the equation of state for the ideal gas is adopted with the adiabatic index $\gamma=5/3$.

We solve the MHD equations by a conservative Godunov-type scheme (e.g., Sano et al., 1998). The hydrodynamical part of the equations is solved by a second-order Godunov method, using the exact solutions of a simplified MHD Riemann problem. The one-dimensional Riemann solver is simplified by including only the tangential component of a magnetic field. The characteristic velocity is the magnetosonic wave alone, and then the MHD Riemann problem can be solved in a way similar to the hydrodynamical one (Colella & Woodward, 1984). The remaining terms, the induction equation and the magnetic tension part in the equation of motion, are solved by the consistent MoC-CT method (Stone & Norman, 1992; Clarke, 1996), guaranteeing $\nabla\cdot\mbox{\boldmath$B$}=0$ within round-off error throughout the calculation (Evans & Hawley, 1988). The numerical scheme includes an additional numerical diffusion in the direction tangential to the shock surface to care for the carbuncle instability (Hanawa et al., 2008).

Calculations are carried out in a frame moving with the velocity $v^{\ast}$ which is the interface velocity after the interaction with the incident shock. For convenience, we define the locus of $y=0$ as where the incident shock reaches the contact discontinuity at $t=0$. Then the post-shocked interface stays at $y=0$ if it is not corrugated initially or stable for the RMI. The simulations start before the incident shock hits the corrugated interface ($|Y_{\rm is}-Y_{\rm cd}|\gg\psi_{0}$).

The system of equations is normalized by the initial density and sound speed of the heavier fluid “1”, $\rho_{1}=1$ and $c_{s1}=1$, and the wavelength of the density fluctuation $\lambda=1$. The sound crossing time is also unity in this unit, $t_{s1}=\lambda/c_{s1}=1$. Most of the calculations in this paper use a standard resolution of $\Delta_{x}$ = $\Delta_{y}$ = $\lambda/256$ unless otherwise stated. A periodic boundary condition is used in the $x$-direction, and an outflow boundary condition is adopted in the $y$-direction. The size of the computational box in the $x$-direction is always set to be $L_{x}=\lambda$. The $y$-length, on the other hand, is taken to be sufficiently wider, $L_{y}\geq 10\lambda$, for the transmitted shock not to reach the edge of the computational domain.

The RMI occurs for any Atwood number or any ratio of the interface densities. When the adiabatic index $\gamma$ is assumed to be identical in both fluids, the light-to-heavy (heavy-to-light) configuration corresponds to the shock-reflected (rarefaction-reflected) case. In our previous work (Sano et al., 2012), the role of the magnetic field on the RMI has been examined intensely only for the shock-reflected cases. However, the RMI can grow when the reflected wave is a rarefaction. Thus we examine the rarefaction-reflected cases in this work and attempt to understand the magnetohydrodynamical evolutions of the RMI comprehensively. In this subsection, we will demonstrate the importance of the rarefaction-reflected cases in a simple system.

Let us suppose a situation when a planar shock front passes through a dense layer of the density $\rho_{b}$ ($>\rho_{a}$), which is depicted by Figure 2(a). Both sides of the dense layer are subjected to the RMI if it is corrugated. The Mach number of the incident shock, $M_{i}=U_{i}/c_{sa}$, determines the growth velocity of the RMI at the front side, $v_{\rm lin,F}$, where $c_{sa}$ is the sound speed outside of the layer. If the decay of the transmitted shock in the dense layer is negligible, the Mach number of the transmitted shock $M_{t}=U_{t}/c_{sb}$ characterizes the growth velocity at the rear surface, $v_{\rm lin,R}$. Assuming the same corrugation wavelength and amplitude at both sides, the ratio of the growth velocities can be derived as a function of the incident Mach number $M_{i}$ and the density ratio $\rho_{b}/\rho_{a}$. Figure 2(b) shows the ratio of the growth velocities at the front and rear sides of the dense layer for a given $M_{i}$. For the cases of weaker shock $M_{i}\lesssim 2$, the growth velocity at the rear side exceeds that of the front side. In general, when the density jump is small ($\rho_{b}/\rho_{a}\lesssim 10$), the RMI at the rear side is as fast as that at the front side, which suggests that both of the shock-reflected and rarefaction reflected RMI are important equally. By contrast, the rear side instability becomes relatively unimportant when the density jump is high ($\rho_{b}/\rho_{a}\gtrsim 100$). In the limit of high Mach number, larger density ratio brings slower growth velocity at the rear side.

The unstable motions of the RMI amplify ambient magnetic fields efficiently for the cases when a rarefaction wave is reflected. This feature is almost identical to that seen in the RMI for the shock-reflected cases (Sano et al., 2012).

The interface is stretched with the magnetic field lines through the growth of the RMI. Figures 3(a) shows the density distribution at the nonlinear stage of the RMI, which is taken at $kv_{\rm lin}t=10$. As for the initial conditions of this fiducial run, the density jump and corrugation amplitude of the interface are $\rho_{2}/\rho_{1}=0.1$ and $\psi_{0}/\lambda=0.1$, respectively. The adiabatic index is $\gamma=5/3$ in both sides of the fluids. For this case, the injected shock wave evolves into a transmitted shock and a reflected rarefaction after the interaction with the contact discontinuity. The Mach number of the incident shock is assumed to be quite large, $M=100$. The phase reversal of the surface modulation is a characteristic behavior for the rarefaction-reflected cases of the RMI. Then, the growth of the density spike in this run is towards the opposite direction to the initial pattern.

A uniform magnetic field is applied in the $y$-direction, and then this is a parallel shock case with $\bm{B}_{0}\parallel\bm{U}_{i}$. The initial plasma beta is set to be $\beta_{0}=10^{8}$. The magnetic field develops passively to the RMI motions because it is not strong enough to alter the plasma flow. The perpendicular shock case ($\bm{B}_{0}\perp\bm{U}_{i}$) with the same field strength is also performed, but the density structure formed by the RMI motions is quite similar except for tiny fluctuations associated with the mushroom-shaped spike.

The magnetic field lines for the parallel and perpendicular shock cases are shown in Figures 3(b) and 3(c). The regions where the magnetic field lines are dense indicate the amplified magnetic field. The strong magnetic fields appear near the interface on the side of the transmitted shock. Small-amplitude fluctuations of the density appear in the lighter fluid, which is associated with reverse and cumulative jets left after the rippled transmitted shock (Stanic et al., 2012; Dell et al., 2015). These motions generate many filamentary structures of locally amplified magnetic fields between the interface and transmitted shock front. This feature can be recognized in both cases, suggesting it is independent of the field direction.

The magnetic field is amplified significantly in the timescale of the RMI, $t_{\rm rm}$. Figures 4(a) and 4(b) show the time evolution of the maximum strength of the magnetic field for the parallel-shock cases ($\bm{B}_{0}\parallel\bm{U}_{i}$) and for the perpendicular-shock cases ($\bm{B}_{0}\perp\bm{U}_{i}$), respectively. The initial parameters except for the magnetic field strength are the same as in the fiducial run shown in Figure 3. The maximum strength is normalized by the initial strength of $B_{0}$ for the parallel-shock case, while it is divided by the post-shocked field strength $B^{\ast}$ in the fluid “2”, which stands for the maximum field strength when the interface is flat. The characteristic behavior of the maximum field strength shows little dependence on the direction of the initial magnetic field. If the ambient magnetic field is weak ($\beta_{0}\gtrsim 10^{4}$), the amplification factor reaches two orders of magnitude. However, as the initial field strength increases, the maximum field strength decreases gradually. When the initial beta is $\beta_{0}=0.01$, no amplification is observed for both cases. These features on the magnetic field amplification are qualitatively similar to the shock-reflected cases (Sano et al., 2012). The strong magnetic fields are localized predominantly in thin filamentary structures. Thus the average magnetic energy over the entire region is at most a few times larger than the initial value.

The magnetic field amplification is tightly linked to the growth of the RMI. The interface is developed into a mushroom shape with smaller-scale vortical structures on the surface. The width of the mixing layer $d_{sb}$, which is defined by the distance from the spike top to the bubble bottom, increases with time. If the external magnetic field is large enough, then the growth of the RMI is suppressed significantly. It is known that the critical field strength depends on the Mach number of the incident shock (Sano et al., 2013). Figure 5 depicts the time history of the mixing width $d_{sb}$ to demonstrate the dependence on the external magnetic field. When the field strength is relatively weaker $\beta_{0}\gtrsim 1$, the growth of the RMI is unaffected by the magnetic field. However, if the plasma beta becomes $\beta_{0}\sim 0.01$, the RMI is stabilized completely. Then, the mixing layer shrinks into much smaller than that in the hydrodynamic case. Although this figure includes only the parallel shock cases, the role of the magnetic field on the RMI suppression is independent of the initial field direction. Note that the critical plasma beta for the RMI is not given by unity.

The critical value of the field strength varies depending on the model parameters, especially on the Mach number (Sano et al., 2013; Matsuoka et al., 2017). However, it is found that a single parameter characterizes magnetohydrodynamic evolutions of the RMI, that is, the Alfvén number for the RMI. The Alfvén number is defined here as the ratio of the growth velocity to the Alfvén speed,

where $v_{A}^{\ast}=\min\{v_{A1}^{\ast},v^{\ast}_{A2}\}$ is the minimum of the Alfvén speed at the interface after the shock passage. The slower Alfvén speed determines the criterion for the entire suppression of the fluctuations on both sides. It is turned out that the RMI growth is diminished when the Alfvén number is less than unity, or when the Alfvén speed exceeds the growth velocity in the unmagnetized limit. The suppression condition $R_{A}\lesssim 1$ is almost the same meaning as that derived by Sano et al. (2013), where the critical magnetic field is obtained by

This equation can be simplifed as $v_{A,{\rm crit}}=\alpha v_{\rm lin}$ where $\alpha$ is a factor of order unity, so that it is strictly related to the Alfvén number. Through a systematic survey of the broad parameter space in this paper, we find this condition is quite robust and valid for any parameters of the RMI, including the perpendicular shock cases and the parallel shock cases.

Figure 6 shows the maximum magnetic field amplified by the RMI as a function of the Alfvén number $R_{A}=v_{\rm lin}/v_{A}^{\ast}$ for various runs. The initial field direction is perpendicular to the interface for all the cases in Figure 6(a). The maximum field strength is measured at the end of calculations, where the magnetic field evolution is saturated sufficiently. The amplification factor is normalized by $B_{0}$. Even though the data include both the shock-reflected cases (filled marks) and rarefaction-reflected cases (open marks), the amplified magnetic field exhibits the same trend. When the Alfvén number is less than unity, the amplification factor is almost unity. No amplification denotes that the magnetic field stabilizes the RMI in this regime. However, as the Alfvén number increases over unity, the maximum strength increases in proportion to $R_{A}$. Then the maximum field energy of $B_{\max}$ becomes comparable to the turbulent kinetic energy of $v_{\rm lin}$ at the saturated state. The amplification factor follows a relation of $1+R_{A}$ at $R_{A}\lesssim 100$. It should be emphasized that the $R_{A}$-dependence is valid for a wide range of parameters; The Mach number is from 1.2 to 10 and the density jump is from 0.1 to 10.

When $R_{A}\gtrsim 100$, the increase of the maximum strength is weakened obviously. In this parameter range, the amplification factor is broadly scattered depending on the parameters, but it looks like a nearly flat function of $R_{A}$. The amplified magnetic field seems to be higher as the incident Mach number increases. However, the maximum field strength in this regime is found to depend on the numerical resolutions. Then, it might be inappropriate to discuss quantitative features of the magnetic field only by our simulations performed in this paper.

The Alfvén number is an excellent indicator of MHD RMI for the parallel shock cases too. The amplification factors of the magnetic field are shown in Figure 6(b) for various cases when the initial ambient field is in the $x$-direction. The maximum field strength is made dimensionless by $B^{\ast}$ for this case. The same trend as in Figure 6(a) can be seen clearly in a range of $R_{A}\lesssim 100$. To summarize, the RMI is suppressed by the tension force of the magnetic field when $R_{A}\lesssim 1$. The surface modulations oscillate stably, and the mixing layer does not grow at all. On the other hand, if the Alfvén number exceeds unity, the magnetic field is enhanced through the interface stretching and tends to be toward the equipartition between $v_{A,\max}\approx v_{\rm lin}$. The saturated level at $R_{A}\gtrsim 100$ indicates that the larger density jump causes the larger amplification in the parallel-shock runs.

The linear growth of the RMI is reduced due to the presence of a strong magnetic field. The other similar instabilities, the Rayleigh-Taylor instability (RTI) and Kelvin-Helmholtz instability (KHI), are also stabilized by the magnetic field. Here we compare the dependence of the field strength on the critical wavelength to suppress these interfacial instabilities. We focus on the simplest situation where two uniform fluids are separated by a horizontal boundary and examine the effect of a uniform horizontal magnetic field $B$.

If the gravity is working toward the lighter fluid ($\rho_{2}$) from the heavier one ($\rho_{1}$), it becomes unstable to the RTI. The linear growth rate of the RTI including the horizontal field is derived as

(Chandrasekhar, 1961) where $g$ is the gravitational acceleration and $k$ is the horizontal wave number. The first and second terms on the right side are the growth rate for the unmagnetized case and the suppression rate originated from the Lorentz force. The RTI is stabilized by the magnetic field when

where $A=(\rho_{1}-\rho_{2})/(\rho_{1}+\rho_{2})$ is the Atwood number, $\bar{v}_{A}=B/(4\pi\bar{\rho})^{1/2}$ is the Alfvén speed, and $\bar{\rho}=(\rho_{1}+\rho_{2})/2$ is the average density. Thus, the fluctuations in shorter wavelengths are more easily stabilized in terms of the RTI.

The interface is subject to the KHI, if the two fluids have relative horizontal motions ($v_{1}\neq v_{2}$). The linear growth rate is derived analytically as

(Chandrasekhar, 1961). Here a uniform magnetic field parallel to the interface is assumed. The stabilization by the magnetic field is determined from the balance between the two terms on the right side of $\sigma_{\rm kh}$. Unlike the case of the RTI, the criterion is independent of the wavelength. The critical field strength is calculated as

for all wavelengths of perturbations.

The suppression condition for the RMI is obtained empirically by the numerical simulations. The growth of the RMI is dramatically reduced when the Alfvén number is less than unity. The linear growth velocity can be expressed by $v_{\rm lin}=\zeta k\psi_{0}U_{i}$ using the incident shock velocity $U_{i}$. The factor $\zeta$ is calculated from the pre-shocked quantities and takes roughly of the order of 0.1; $\zeta=\zeta(M,A,\gamma)\sim O(0.1)$. Then, the RMI is suppressed by the magnetic field, when $v_{\rm lin}<v_{A}^{\ast}$. This condition is rewritten as

Therefore, for the RMI, the surface fluctuations with longer wavelengths are preferentially stabilized as oppose to the RTI.

The characteristics of the stabilized wavelength are distinctly different for each instability. To be precise, though, this condition should depend on the direction of the magnetic field. Especially, interchange modes are unaffected at all by the restoring force of magnetic fields. However, the relations between the critical wavelength and magnetic field [Equations (11), (13), and (14)] are beneficial to intuitively understand the role of the magnetic field for each instability.

Single-mode analysis has the advantage of understanding the nonlinear features of the RMI in a clear picture. However, the actual situations are not always as simple as such a system. The growth process of more complex multi-mode disturbances is briefly touched upon in the last section.

The interface modulation for our multi-mode analysis is assumed to be

where $\lambda_{n}={\lambda}/{n}$ is the multi-mode wavelength, $a_{n}$ and $\phi_{n}$ are random coefficients between 0 and 1, $\sigma$ is the standard deviation of the summation part to normalize the fluctuation to be $1/\sqrt{2}$, and the initial amplitude $\psi_{0}$ determines the size of the perturbation. The single-mode analysis corresponds to the case of $N=1$. The multi-mode corrugation with $N=10$, for example, consists of ten different wavelengths of $n=1$–10.

The magnetic field amplification for the case of multi-mode fluctuations is examined in Figure 7. The model parameters are mostly the same as the fiducial run of the single-mode case; $M=100$ and $\rho_{2}/\rho_{1}=0.1$. The position of the initial interface is given by Equation (15) with $\psi_{0}/\lambda=0.1$ and $N=10$. As a result of the RMI growth, the interface shows a complicated shape where various-scale mushroom-structures are growing randomly in all directions. However, the thickness of the mixing layer is comparable to that in the corresponding single-mode case shown by Figure 3.

The initial magnetic field is relatively weak as $\beta_{0}=10^{8}$ so that the RMI motions enhance the seed field in this run. The amplified field structure exhibits that a large number of curled filaments fill the entire area between the interface and the transmitted shock front [see Figure 7(b)]. The filamentary feature is caused by the complicated vortex structure left behind the transmitted shock. Figure 8 indicates the probability distribution function of the magnetic field strength for the cases of $N=10$, 3, and 1. The parameters for the interface modulation are listed in Table 1. Here we assume $a_{n}=1$ for a white noise amplitude. The maximum field strength in the multi-mode runs is nearly the same as that in the single-mode. Furthermore, the power-law feature seen in the amplified component is shared by the single-mode and multi-mode runs. The index fitted in a range of $2<|\bm{B}|/B_{0}<20$ is about $-1.7$, which may be the characteristic quantity of the RMI turbulence driven by a shock wave. The power-law features are unaffected by choice of the random numbers for the initial interface shape in Equation (15). It should be noticed that a similar power index is obtained in interstellar turbulence simulations by Inoue et al. (2009), indicating the RMI motions could contribute to the evolutions of interstellar magnetic fields. However, three-dimensional simulations will be essential for further quantitative comparisons with observations.