The Rayleigh-Taylor instability (RTI) Lord (1900); Taylor (1950); Chandrasekhar (1961), an important hydrodynamic instability, occurs at the unstable interface when a high density fluid is supported by a lower density fluid under the influence of gravity, or when the interface between two fluids with different densities is accelerated. This instability is ubiquitous in nature and plays an important role in diverse areas of science and technology, including inertial confinement fusion (ICF) Tabak et al. (1990); Zhou et al. (2019); Remington et al. (2006); Betti et al. (1998); Srinivasan & Tang (2012); Srinivasan et al. (2012Srinivasan et al. ( , 2017Srinivasan et al. ( , 2019)); Srinivasan & Tang (2014b,a); Stone & Gardiner (2007); Wang et al. (2017), astrophysics Gamezo et al. (2003); Kifonidis, K. et al. (2003); Hwang et al. (2004); Hester (2008); Loll et al. (2013), geophysics Kaus & Becker (2007), and engineering processes Lyubimova et al. (2019). For instance, the RTI is known to act as an inhibitor in achieving an ignition grade hot spot in ICF targets Srinivasan & Tang (2012); Srinivasan et al. (2012Srinivasan et al. ( , 2019)); Srinivasan & Tang (2014b,a); Zhou (2017a,b). RTI occurs in ICF targets during both the acceleration and deceleration phase of the implosion, leading to undesirable mixing of hot and cold plasmas. The RTI is also observed in various astrophysical phenomena such as supernova explosions and their remnants (Crab Nebula) Gamezo et al. (2003); Kifonidis, K. et al. (2003); Hwang et al. (2004); Hester (2008); Loll et al. (2013). Therefore, a detailed understanding of such instabilities in high-energy-density (HED) plasmas has implications for ignition-grade hot-spots, understanding supernova explosions, and revealing Mega-Gauss (MG) scale magnetic field generation and their turbulence in astrophysical settings. The RTI and their mitigation mechanism in HED plasmas has been thoroughly studied by several authors experimentally as well as theoretically and numerically Remington et al. (2006); Srinivasan & Tang (2012); Srinivasan et al. (2012Srinivasan et al. ( , 2019)); Srinivasan & Tang (2014b,a); Atzeni & Meyer-ter Vehn (2004); Sun et al. (2021); Silveira & Orlandi (2017). However, there exists a substantial disagreement between computer simulation results and highenergy density laboratory experiments or astrophysical observations of the RTI Kuranz et al. (2010); Modica et al. (2013). Most of the experiments or astrophysical observations have noted unusual morphological structure of RTI which are significantly different from the computer simulation results, exhibiting strongly suppressed growth of small scale structures and mass extensions of RT spikes. This is due to the fact that many theoretical and numerical studies use conventional hydrodynamic and magnetohydrodynamic (MHD) depiction where either the self-consistent effect of magnetic fields, viscosity, and resistivity have been ignored or they have been considered in isolation. First observations of the the magneto-Rayleigh-Taylor instability evolution in the presence of magnetic and viscous effects have been made in recent experiments Adams et al. (2015). The impact of magnetic fields on RTI in the presence of self-consistent viscosity and resistivity for experimentally and observationally relevant parameter regimes in HED plasmas remains an open area of research.
The primary purpose of this paper is, therefore, to understand the role of the viscous and resistive effects on RTI in magnetized HED plasmas applicable to astrophysical plasmas as well as ICF-based laboratory experiments. Specifically, this work aims to understand how RTI dynamics is impacted by varying plasma beta (ratio of thermal energy to magnetic energy) and magnetic Prandtl number (ratio of magnetic Reynolds number to Reynolds number). This study focuses on a high Atwood number and high-β regime, where the energy density in the magnetic field is small compared to the thermal energy in the fluid. The Atwood number (A t ) is a dimensionless number defined as, A t = (ρ H -ρ L )/(ρ H + ρ L ); where ρ H and ρ L represent the mass density of the heavy and light fluid, respectively. This distinguishes the current work from previous works that have examined the role of viscosity and resistivity in isolation for ICF applications Srinivasan & Tang (2014a); Song & Srinivasan (2020). In addition, this work also presents the evolution of RTI considering fully varying self-consistent viscosity and resistivity profiles. To study the RTI dynamics in HED plasmas, the magnetohydrodynamic (MHD) equations with the inclusion of viscosity and resistivity are solved in this work. These visco-resistive MHD equations are solved in conservation form in 2D (two dimensions) using the fluid modeling tool PHORCE (Package of High ORder simulations of Convection diffusion Equations) based on the unstructured discontinuous Galerkin finite element method Song (2020); Song & Srinivasan (2021); Hesthaven & Warburton (2007). Under this configuration, simulations have been performed over a wide range of magnetic Prandtl numbers with the presence of a longitudinal external magnetic field to reveal the effect of viscosity and resistivity on the evolution of RTI and magneto-RTI in HED plasmas. It is observed that the inclusion of viscosity and resistivity dramatically changes the growth as well as the structures/morphology of the instability on different length scales. It is seen here that the presence of viscosity stabilizes the growth of the RTI and modifies the morphology of the tip of RTI fingers, inhibiting the traditional mushroom cap structures. On the other hand, the morphology of the RTI spikes is found to be independent of resistivity. The presence of resistivity assists in the development of small scale structures by reducing the magnetic field stabilization. When considering spatially-varying viscosity and resistivity with highly disparate profiles, there is a significant impact on the RTI evolution in the high Atwood number regime studied in this work. In this paper, the numerical growth rates of RTI obtained from the simulations are compared with their corresponding analytical values obtained from linear theory. Furthermore, it is also seen here that the dynamics of magnetic field is independent of viscosity and likewise the resistivity does not affect the dynamics of enstrophy and kinetic energy. In addition, this work presents the power law scaling of enstrophy, kinetic energy, and magnetic field energy in both the injection range and inertial sub-range of power spectra for different viscosity and resistivity cases, which could be useful for understanding the RTI induced turbulent mixing in HED plasmas.
The manuscript has been organized as follows. In Section 2, a brief description of the governing equations is presented to study the RTI process in magnetized HED plasmas. Section 3 discusses the simulation techniques and problem setup for the study. Section 4 presents the simulation results, comparison with theory, and discussions. Section 5 presents the summary and conclusion.
In this section, the basic governing equations are presented for the study of RTI in magnetized HED plasmas in the presence of an applied horizontal magnetic field, viscosity, and resistivity. Thermal conduction is neglected in this study to focus on the impact of viscosity, resistivity, and magnetic fields. The generalized Lagrange multipliermagnetohydrodynamic (GLM-MHD) equations Munz et al. (2001); Dedner et al. (2002) with the inclusion of viscosity and resistivity are solved. The compressible MHD equations are given by,
where ρ, u, p, g and B represent the mass density, fluid velocity, pressure, gravitational field, and magnetic field, respectively. Here = p/(γ -1) + ρu 2 /2 + B 2 /2µ 0 defines the total energy; where γ is the ratio of specific heats, and is normally taken as 5/3 for monatomic gases assuming an ideal gas law. Here p is the pressure. For the equation of state, an ideal gas law p = (Z i + 1)ρk B T i /m i is assumed; where Z i , m i , k b , and T i represent the charge state of ion, mass of the ion, Boltzmann constant, and temperature of the ion, respectively. Here ψ, C h , and C p represent the divergence cleaning potential, hyperbolic cleaning speed, and parabolic cleaning speed, respectively. A user-specified parameter C r = C 2 p /C 2 h is defined to determine the ratio between hyperbolic and parabolic divergence cleaning. If C r is very large, the divergence error will only be transported through the hyperbolic term. C h is calculated based on the grid sizes and CFL number (Dedner et al. 2002). In the simulations presented here, C r = 99999 is set to be very large so that only hyperbolic cleaning dominates. In the above equations π and η represent the viscous stress tensor and electrical resistivity coefficient, respectively. In this study, the Braginskii formulation Braginskii (1965) for calculating viscosity and resistivity co-efficient is used, µ = 0.96n i k B T i τ i and η = m e /1.96n e q 2 e τ e , where τ e and τ i are the collision times for electron and ion, respectively. Note that the viscosity and resistivity can also be presented in terms of Reynolds (Re) and magnetic Reynolds number (Re m ) defined as, Re = ρV L/µ and Re m = µ 0 V L/η; where V and L represent some reference velocity and length, respectively.
This section presents the simulation techniques and problem setup used for studying the role of viscosity and resistivity on RTI in magnetized HED plasmas. The simulations presented here are in planar geometry and in 2D. A significant amount of insight can be gained from 2D studies particularly where observations may be dominated by 2D evolution of perturbation growth. In other words, this is true when the wavelength of perturbation for RTI growth in the considered directions is much smaller than the wavelength of perturbation in the third direction. This approximation would be particularly well-suited for cases where magnetic fields influence RTI growth leading to regimes where the perturbation growth are more "2D-like". Most of the past literature on 2D MHD turbulence, not specific to RTI, has focused on incompressible MHD models Orszag & Tang (1979); Biskamp & Schwarz (2001) whereas this work uses a compressible MHD model with a focus on the evolution of the RTI. However, a fully 3-D RTI turbulence study would be important to understand the RTI induced turbulence accounting for 3-D perturbations and this would constitute future studies. In this paper, the code PHORCE (Package of High ORder simulations of Convection diffusion Equations) Song (2020); Song & Srinivasan (2021) developed at Virginia Tech is used for the 2D RTI study. PHORCE is based on the nodal unstructured discontinuous Galerkin method Hesthaven & Warburton (2007) and solves fluid equations (2.1-2.5) in conservation form. To advance the simulation in time, an explicit fourth-order five-stage strong stability-preserving Runge-Kutta (SSP-RK) Song (2020) scheme has been implemented. Several limiters and filters are applied in PHORCE to preserve the positivity of density and pressure and to diffuse the numerical oscillations that typically occur due to strong discontinuities. The code uses an affine reconstructed discontinuous Galerkin (aRDG) scheme Song (2020); Song & Srinivasan (2021) to solve the diffusion terms and to self-consistently capture the effect of spatially varying Reynolds numbers (viscous effects) and magnetic Reynolds numbers (resistive effects).
The RTI simulations have been performed in a rectangular domain with x ∈ [-L x /2, L x /2], y ∈ [-L y /2, L y /2]; where L x and L y represent the width and height of the simulation domain, respectively. The simulations are performed with 2000 × 1000 cells. The gravitational field g = -g ŷ. The simulations are performed using "conducting wall" boundary conditions along the y-direction and "periodic" boundary condition along the x-direction. In equilibrium, the simulation is initialized using the standard hyperbolic tangent density profile given by,
In the above equation, α defines the width of the hyperbolic tangent function. In the simulations presented here α is taken to be 0.01 in order to provide a sharp gradient at the interface. The pressure profile is initialized as,
where p 0 represents the background pressure of the system. To excite the multimode RTI in the simulation, the y-component of velocity at the interface (y = 0 plane) is perturbed as, v = Σ 40 m=1 0.01R(m)cos(2π(mx/Lx + R(m)))exp(-ξy 2 ) at t = 0; where R(m) and ξ represent the random number generator function of m random numbers and the spatial width along the y direction over which the perturbation falls at the interface, with ξ = 1000.
In this work, all the simulation results are presented in normalized units. The following normalization factors have been used,
Here γ RT = √ A t gk, represents the growth rate of the RTI associated with the wave number k = 2π/λ; where λ is the wavelength of the mode Chandrasekhar (1961). As the simulations have been conducted with multimode perturbations having mode number m = 1 -40, note that the value of γ RT would be different for different modes (or wavelengths). The growth rate becomes maximum for smallest wavelength and minimum for longest wavelength modes. To calculate the the value of γ RT for the normalization of time, the smallest mode of perturbation (m = 40) having wavelength λ = L x /40 has been selected.
In some flows in HED plasmas, such as in ICF and supernovae explosions Sauppe et al. (2019); Burton (2011); Dimonte et al. (2005); Cabot & Cook (2006); Srinivasan & Tang (2014a); Srinivasan et al. (2012), the Atwood number can reach a very high value (A t 0.85) and the temperature can have a large variation in the domain. As a result, a large variation in Reynolds and magnetic Reynolds numbers may exist in the domain. In this work, the plasma parameters are selected to access highly varying density and temperature regimes in laboratory and astrophysical plasmas where the viscosity and resistivity may be important. The parameters are summarized in Table 1 in normalized form. The simulations use an initial plasma beta β ini = 2µ 0 p 0 /B ext x 2 = 5000 whenever an external horizontal magnetic field (B ext x ) exists in the system. Using the parameters given in the Table 1 and using the expressions for isotropic . Note that the resistivity profile is modified in the heavy fluid to increase the minimum value of the magnetic Reynolds number from 20 to 285. For the simulations presented here, the modified resistivity profile has been used to capture the essential physics of RTI in the presence of resistivity. The magnetic Prandtl number, P r m = Re m /Re = ν/η (where ν = µ/ρ is the kinematic viscosity), is a dimensionless quantity that estimates the ratio of momentum and magnetic diffusivity. In Fig. 1b, P r m varies from 2 for y/L x < 0 to 4 × 10 -6 for y/L x > 0 producing a significant variation across the domain.
The simulations have been performed for different values of magnetic Prandtl numbers to elucidate the role of viscosity and resistivity on the Rayleigh-Taylor and magneto-Rayleigh-Taylor instability. Table 2 summarizes all simulation cases performed here for different values of plasma beta (external magnetic field) and magnetic Prandtl numbers (Reynolds numbers and magnetic Reynolds numbers). This section discusses the results and findings of each case that is presented.
Simulations for inviscid (µ = 0) and irresitive (η = 0) cases are performed (see run-1 and run-2 in Table 2). Fig. 2 presents plots of mass density (ρ/ρ L ) at different times for β ini → ∞ (no initial external horizontal magnetic field) and for β ini = 5000 (in the presence of initial external horizontal magnetic field). As expected, the height of the RTI mixing region or the height of the RTI fingers reduces in the presence of an Table 2. Summary of numerical simulations performed here. applied horizontal magnetic field. Note ths suppression of small-scale structures due to the presence of the magnetic field. To calculate the growth rate, the peak bubble-to-spike distance (h/L x ) over the normalized times (tγ RT ) for both β ini → ∞ and β ini = 5000 is presented in Fig. 3. In the simulations, the height has been measured by tracking the difference of the upper and lower boundary of the RTI mixing region. As shown in the subplot of Fig. 3, the numerical growth rates are calculated from the slope of the plot log(h/L x ) vs. tγ RT . The numerical growth rate obtained from the simulations for both β ini → ∞ and β ini = 5000 are 0.75γ RT and 0.5γ RT , respectively. The growth of RTI significantly decreases in the presence of applied horizontal magnetic field as expected. The analytical expression of growth rate (γ RT ) of RTI for purely hydrodynamic flows (no viscosity, no resistivity, and no magnetic field) is given by Chandrasekhar (1961) as, γ RT = A t gk.
(4.1)
Using the parameters given in Table 1 and k = 80π/L x (for λ = L x /40), the analytical values of the growth rate γ RT can be estimated as 2.69 × 10 9 s -1 for a single mode that is estimated to be the fastest growing early in time. The numerical growth rate is 0.75γ RT = 2 × 10 9 s -1 but this is for a multimode growth rate which explains the difference between the analytical and numerical values. As time evolves, the nonlinear interactions between modes significantly changes the dominant wave number. When an applied magnetic field B ext exists, the RTI growth rate becomes Chandrasekhar (1961); Jun & Norman (1996),
Note that the RTI is affected by the horizontal magnetic field (B k) and is not directly impacted by magnetic fields that are normal to the interface when using an MHD model. In Fig. 2 for β ini = 5000, the height of the mixing region is decreased along with suppression of the small structures. In this case, one can approximately calculate the wavelength of RTI fingers by calculating the number of RTI fingers in the domain. This technique suggests approximately 30 RTI spikes growing at this time. Therefore, the effective smallest wavelength is approximately ≈ L x /30. When an appropriately aligned magnetic field is initialized, the value of the peak magnetic field in the system increases with time as RTI grows. For example, the plasma β becomes 226 from an initial value of 5000 at time tγ RT = 13.5. Using the parameters given in Table 1, β = 226, and k x L x = 60π, the analytical values of the growth rate γ B RT can be estimated as γ B RT = 0.63γ RT . The numerical growth rate obtained from the simulation shows good agreement with the analytical value for β ini = 5000 considering that these are estimates for multimode simulations.
The enstrophy (Z), kinetic energy (E), and magnetic field energy (B 2 ) averaged over the vertical direction (y) of the system is defined as,
where ω = ∇×u represents the fluid vorticity. In 2D mixing and turbulence, the enstrophy (Z), kinetic energy (E), and magnetic field energy (B 2 ) are important quantities as they appear to be the only quadratic constants of motion. In Fig. 4, the evolution of enstrophy (Z) and kinetic energy (E) spectra are presented as a function of normalized wave number k x L x at different times for β ini → ∞. Note that there will be no magnetic field for β ini → ∞.
The spectra can be separated into three regions based on the range of k x L x . The first region with k x L x 80π is known as the injection range where all the external perturbation modes exist. All the energy has been injected to the system within these wavelengths. The second region 80π k x L x 600π or the middle range is the inertial sub-range. This is the regime which basically connects the injection range to the dissipation range. The third region where k x L x 600π is the dissipation range which accounts for grid scales as L x = 2000∆x; where ∆x is the grid size along the x-direction. As physical dissipation (viscosity and resistivity) is absent in the system for the simulations in this section, the only dissipation mechanism is, therefore, governed by the numerical dissipation. All the energy for modes smaller than or equal to the grid size is dissipated by numerical dissipation. For β ini → ∞ (see Fig. 4), note that the enstrophy (Z) and kinetic energy (E) increase equally in all the available modes in the system with time as long as tγ RT 13.5. At tγ RT = 17.5, the transfer of kinetic energy as well as enstrophy is seen from short wavelength modes to long wavelength modes. This happens due to the nonlinear interactions of the modes leading to the formation of longer wavelength modes with time. As a result, the small scale structures get modified changing the growth rate in the nonlinear regime for β ini → ∞. The numerically obtained power law scalings for the enstrophy, kinetic energy, and magnetic field energy spectra at both the injection and inertial sub-range are included in Fig. 4. In this case, the spectra of kinetic energy and enstrophy obey the following power scaling laws in the injection range (k x L x 80π),
x .
(4.5)
In the inertial sub-range (80π k x L x 600π), the spectra are found to obey different power laws,
For β ini = 5000, the evolution of enstrophy (Z), kinetic energy (E), and magnetic field energy (B 2 ) spectra as a function of wave number k x L x at different times is shown in Fig. 5. The enstrophy (Z), kinetic energy (E), and magnetic field energy (B 2 ) increase equally in all modes in the system until tγ RT = 17.5 for β ini = 5000. There is no transfer of kinetic energy, enstrophy, and magnetic field energy over the modes. This is because the spectrum still lies in the linear regime due to the presence of a horizontal magnetic field. The magnetic field opposes the growth of the RTI and decreases the vertical velocity of the fluid. In this case, the spectra of kinetic energy, enstrophy, and magnetic field energy, obtained from the numerical simulations, obey the following power laws in the injection range (k x L x 80π)
(4.10)
Similarly, the power law in the inertial sub-range (80π k x L x 600π) for β ini = 5000 is found to be
x (4.12)
The slope of the spectra in the inertial sub-range decreases with the presence of a horizontal magnetic field. The slope of the inertial sub-range measures the rate at which the energy is transferred from large scale to small scales or vice versa. In other words, it defines the rate at which the larger scales get fragmented into smaller scales and vice versa due to mixing. Therefore, this shows that the rate of small scale formation due to
x Table 3. Summary of power laws for the numerical simulations in injection range.
RTI mixing decreases with the application of horizontal magnetic field. The scaling of these power laws in both injection and inertial sub-range for these cases (run-1-2) have been summarized in Tables 3 and4. Note that the numerical dissipation is active in the range of k x L x > 600π. As a result, all the energy is also seen to grow proportionally with time in this regime.
4.2. Simulation results for constant viscosity, irresistive cases (P r m = ∞): run-3-6
Constant viscosity is introduced throughout the domain in the simulation. The simulations are performed for two different values of constant Reynolds numbers, Re = 2 × 10 3 and Re = 2×10 6 , but with no resistivity (Re m = ∞). As η = 0 for these simulations, this study corresponds to the cases of very large magnetic Prandtl number (P r m → ∞). In this study, the case without magnetic field β ini → ∞ and with magnetic field β ini = 5000 at tγ RT = 0 are considered. The relevant simulation parameters are shown in Table 2 under run-3-6. In Fig. 6, the mass density (ρ/ρ L ) is shown at different times for Re = 2 × 10 3 and Re = 2 × 10 6 for β ini → ∞. It is seen that the growth of the RTI decreases with decreasing Re or increasing viscosity (µ). Fig. 7 shows mass density (ρ/ρ L ) at different times for the two Reynolds numbers Re = 2 × 10 3 and Re = 2 × 10 6 , but for β ini = 5000. Here, the size of the RTI fingers decreases further when applying a horizontal magnetic field compared with the inviscid case presented in Section 4.1. The magnetic field has a stabilizing effect in addition to viscous stabilization on the growth of RTI. To further illustrate the complementary role of viscous and magnetic field stabilization, the runs inertial sub-range power law run-1
x Table 4. Summary of power laws for the numerical simulations in inertial sub-range. the simulations is 0.55γ RT and 0.64γ RT for Re = 2 × 10 3 and Re = 2 × 10 6 , respectively. The analytical expression for the growth rate of RTI in a compressible viscous fluid is given by Menikoff et al. (1977),
where ω = ν2 k 3 /A t g and ν = (µ l + µ h )/(ρ l + ρ h ) is the density averaged kinematic viscosity. In Fig 9, the analytical form of γ vis RT /γ RT is shown as a function of wave number k x L x for Re = 2 × 10 and Re = 2 × 10 6 . For Re = 2 × 10 3 , it is seen that the analytical growth rate is maximum for k x L x ≈ 60π which corresponds to a wavelength of approximately L x /30. Similarly, for Re = 2 × 10 6 , the analytical growth rate becomes maximum for k x L x ≈ 76π or a wavelength of approximately L x /38. This is consistent with the simulation results from Fig. 7. The theoretical growth rate of the mode having wavelength L x /30 and for the mode having wavelength L x /38 are approximately 0.56γ RT and 0.65γ RT , respectively. The growth rates obtained from simulations show good agreement with the analytical results.
Note that, when viscosity increases, the morphology of the RTI spikes appear to be smooth and exhibit different characteristics as seen in Fig. 7. Due to the presence of viscosity, the traditional mushroom cap structure on the tip of the RTI fingers gets inhibited and forms smooth structures. The presence of viscosity also strongly suppresses the growth of the small scale structures and short-wavelength modes.
The plasma β as a function of peak bubble-to-spike distance (h/L x ) for different Re for β ini = 5000 is presented in Fig. 10. Note that plasma β is independent of Re if presented as a function of the peak bubble-to-spike amplitude instead of as a function of time. This shows that the dynamics of magnetic field is not affected by the viscosity for the same amplitude of the RTI growth but the actual RTI growth as a function of time is impacted by the different Re as noted from Fig. 8. Also note that plasma β decreases with time or height as RTI grows for all Re considered. This is because the value of magnetic field increases as RTI grows in the system. Figure 11 presents enstrophy(Z), kinetic energy (E) and magnetic field energy (B 2 (k)) spectra at time tγ RT t = 17.5 as a function of wave number k x L x for different values of Re. The scaling of these power laws in both injection and inertial sub-range for these cases (run-3-6) have been summarized in Tables 3 and4. Note that the spectral power of the magnetic energy does not change with Re but the spectral power of enstrophy and kinetic energy increases with increasing the value of Re for all available modes. This shows that the dynamics of magnetic field energy is independent of Re or viscosity. It is shown by Kulsrud et al. (1997) that the dynamics of the magnetic field can be completely described by ion fluid vorticity in the absence of viscosity and resistivity but in the presence of a Biermann battery, which is not considered in this work. Including the viscosity and resistivity into the MHD equations considered here, a theoretical treatment is included to illustrate the dynamics of magnetic field and vorticity in presence of viscosity and resistivity. Following the same method as shown by Kulsrud et al. Kulsrud et al. the momentum equation (Eq. 2.2) can be written in terms of vorticity (ω) as, where J represents the net current density. Similarly, Eq. ( 2.4) can be modified in terms of ion cyclotron frequency (ω ci = Z i eB/m i ) as,
The last term on right hand side of equations (4.15) and (4.16) are responsible for the dissipation of the vorticity and magnetic field, respectively. The dynamics of vorticity and kinetic energy depend on the viscous stress tensor π and the corresponding Re. This is consistent with the numerical results presented here. On the other hand, the dynamics of vorticity is independent of resistivity η or magnetic Reynolds number Re m , but the dynamics of magnetic field depends on the Re m . To illustrate this, simulations are performed for different constant values of Re m discussed in the next section. In this section, simulation results are presented for different constant magnetic Reynolds numbers (Re m ) but with no viscosity (µ = 0) (see run 7-8 in Table 2). In this study, P r m = 0. In all these simulations, an initial horizontal magnetic field with β ini = 5000 is applied. In Fig. 12, the mass density (ρ/ρ L ) profile at different times is presented for Re m = 285 and Re m = 1105. It is seen that the growth of the RTI increases with a decrease in magnetic Reynolds numbers (Re m ) or increase of resistivity (η). This is because the resistivity diffuses the magnetic fields and reduces the magnetic stabilization. As a result, the RTI growth increases due the reduction of effective magnetic field tension. In this figure, it is to be noted that the morphology of the RTI spikes in terms of mushroom cap structures on the tip of the fingers are seen to be independent of Re m . Also of note is the appearance of additional small scale structures for higher resistivity cases. This is also expected as the magnetic field opposes development of the small scale structures. In Fig. 13 the effect of magnetic Reynolds number on the growth rate of RTI in HED plasmas. It is found that the growth rate increases with increase in resistivity. The numerical growth rates are obtained from the simulations for Re m = 285 and Re m = 1105 as 0.68γ RT and 0.53γ RT , respectively. Including a finite constant resistivity η, Jukes (1963) has shown that the analytical growth rate of RTI changes with resistivity η as,
(4.17)
The growth rates obtained from the simulations also obey the analytical scaling. The plasma β is plotted as a function of peak bubble-to-spike distance (h/L x ) for different Re m in Fig. 14. Note that plasma β decreases with peak bubble-to-spike distance for all values of Re m but at different rates depending on the value of Re m . The rate at which the plasma beta decreases is larger for high Re m . This shows that the dynamics of the magnetic field is not independent of resistivity. This is due to the fact that the magnetic field gets diffused more for low Re m leading to a higher plasma β.
In Fig. 15, the plot of enstrophy(Z), kinetic energy (E), and magnetic field energy (B 2 (k)) spectra at time (tγ RT t = 17.5) as a function of wave number k x L x has been shown for different values of Re m . The scaling of these power laws in both injection and inertial sub-range for these cases (run-7-8) have been summarized in Tables 3 and4. It is observed that the magnetic field spectra changes significantly by changing the value of Re m , whereas the spectra of enstrophy and kinetic energy does not show any significant dependence on the value of Re m . The spectral power of magnetic field energy increases with increasing the value of Re m for all the available modes. This justifies that the dynamics of magnetic field energy depends on Re m or η. But the dynamics of enstrophy and kinetic energy does not depend on Re m . This is consistent with equations (4.15) and (4.16).
Simulations have also been performed for different values of constant Re with the inclusion of different constant values of Re m (see run-9-12 in Table 2). In this case, all the simulations use an applied horizontal magnetic field corresponding to β ini = 5000. Fig. 16 ) for Re = 2 × 10 6 is presented in Fig. 17. Note that the morphology of the RTI fingers doesn't exhibit a strong dependence on Re m for the values considered here, but shows more significant dependence with Re. The mushroom caps on the tip of the RTI fingers are inhibited for high viscosity. When viscosity is held constant, the growth rate increases with an increase in resistivity. On the other hand, the growth rate decreases with an increase in viscosity when resistivity is held constant.
The power law scaling of the enstrophy (Z), kinetic energy (E), and magnetic field energy (B 2 (k)) spectra as a function of wave numbers k x L x is quantified for these runs (run-9-12) in both the injection range as well as the inertial sub-range. The scalings of these power laws are given in Tables 3 and4 for run-9-12. with Re. Similarly, the enstrophy(Z), kinetic energy (E), and magnetic field energy (B 2 (k)) spectra at time tγ RT t = 17.5 as a function k x L x for different values of Re m are plotted in Fig. 19 holding Re constant at Re = 2 × 10 3 for β ini = 5000. Note that Re m does not affect the spectra of enstrophy and kinetic energy, whereas the magnetic field spectra depends on Re m . These findings are consistent with those in Sections 4.2 and 4.3. 4.5. Simulation results for fully varying viscosity, irresistive cases (P r m = ∞):
run-13-14
Next, the self-consistent fully varying Re profile shown in Fig. 1 is considered without resistivity (see run-13-14 in Table 2). The simulations have been performed using both β ini → ∞ and β ini = 5000. In Fig. 20, the mass density (ρ/ρ L ) profile is presented at different times for β ini → ∞ and β ini = 5000. To further illustrate the effect of a fully varying Re profile on the RTI, the peak bubble-to-spike distance (h/L x ) over time (tγ RT ) is presented for β ini → ∞ and β ini = 5000 for this case in Fig. 21 along with the bubble-to-spike amplitudes for constant Re cases. The growth and nature of the RTI for fully varying viscosity for β ini → ∞ and β ini = 5000 is close to that of the high viscosity case or low Reynolds number (Re = 2 × 10 3 ) case. This is because, the RTI fingers largely grow in the lower density regime (y < 0) at the interface due to the high Atwood number considered here. The mixing is not significant in the high density regime. In the lower density regime, the value of Re = 2 × 10 3 , which has significantly higher viscosity compared to the high density regime (y > 0). Therefore, the evolution of RTI is dominated by the high viscosity regime. Hence, viscosity, even if disparate, plays an important role in the RTI process in such parameter regimes with and without an applied horizontal magnetic field. Similar to the previous cases, the power law scaling of enstrophy (Z), kinetic energy (E), and magnetic field energy (B 2 ) spectra as a function of wave number k x L x in both injection and inertial sub-range are summarized in Tables 3 and4 under run-13-14. In these studies, an applied horizontal magnetic field corresponding to β ini = 5000 is included as before. Fig. 22 shows the mass density (ρ/ρ L ) profile at different times for Re m = 285 and Re m = 1105. It is seen that the growth of the RTI spikes increases with the decrease of Re m as expected. In Fig. 23, the peak bubble-to-spike distance (h/L x ) over time (tγ RT ) for different values of Re m is shown. Note that the growth of RTI is higher for high resistivity (blue solid line) compared to that obtained for low resistivity (red solid line) when also including the fully varying viscosity. The plasma β as a function of peak bubble-to-spike distance (h/L x ) for different values of constant Re m (solid blue and red line) is presented in Fig. 24, where β ini = 5000 and fully varying Re are considered. The magnetic field decreases for the lower value of Re m = 285 which corresponds to higher η. Furthermore, it is observed here that the morphology of the RTI fingers are not significantly affected by the resistivity. The power law scaling of enstrophy (Z), kinetic energy (E), and magnetic field energy (B 2 ) spectra as a function of wave number k x L x in both injection and inertial sub-range for these cases is summarized in Table 3 and 4 in the column under run-15-16. The final set of simulations are performed for a fully varying Re along with a fully varying Re m profile. These correspond to P r m ranging from 2 -4 × 10 -6 . Note that the resistivity profile used for this case is the modified resistivity profile shown in Fig. 1. In this case, an applied horizontal magnetic field corresponding to β ini = 5000 is included as in the previous cases. of the RTI spikes exhibits less Kelvin-Helmholtz formation and shows the suppression of small scale structures more significantly than the higher Re m = 1105, fully varying Re case presented in Fig. 22. In Fig. 23, the peak bubble-to-spike distance (h/L x ) over time(tγ RT ) for fully varying Re m and fully varying Re profiles is presented along with the constant Re m cases (see yellow solid line). The growth rate for the fully varying resistivity case is close to the the growth rate obtained for the constant Re m = 1105 case. This is because the RTI mostly grows in the low density regime where Re m = 1105. Therefore, the dynamics of RTI for the high Atwood number regime can be described by the physical parameter space of the lower fluid, which is governed by the viscosity and resistivity of the lower fluid. The plasma β as a function of peak bubble-to-spike distance (h/L x ) for fully varying Re m and Re is shown in Fig. 24 (see solid yellow line), where β ini = 5000 is considered. The dynamics of the magnetic field and its corresponding growth, as noted by the decreasing plasma β, for fully varying Re m and Re is different from the constant magnetic Re m cases. The field strength obtained lies inbetween the regimes of the upper and lower fluid (with their corresponding resistivities). The power law scaling of enstrophy (Z), kinetic energy (E), and magnetic field energy (B 2 ) spectra as a function of wave number k x L x in both injection and inertial sub-range is summarized in Table 3 and 4 under run 17.
In summary, the role of viscosity and resistivity on Rayleigh-Taylor and magneto-Rayleigh-Taylor instabilities is studied for a high Atwood number and high plasmaβ regime in high-energy-density (HED) plasmas applicable to both laboratory and astrophysical settings. This work describes 2D RTI evolution and resulting turbulence when surveying plasma-β and magnetic Prandtl number, P r m , for these regimes. The simulations have been performed using fluid simulation techniques based on the unstructured discontinuous Galerkin finite element method Song (2020); Song & Srinivasan (2021); Hesthaven & Warburton (2007). Using a visco-resistive-magnetohydrodynamic (MHD) model, a detailed investigation of RTI in 2D planar geometry for experimentally and observationally relevant parameters is presented. It has been shown here that the inclusion of viscosity and resistivity in the system drastically changes the growth of the instability as well as modifies its internal structure on smaller scales. The presence of viscosity inhibits the development of small scale structures and significantly modifies the morphology of the RTI spikes. On the other hand, the morphology of the RTI spikes is found to be independent of resistivity but it assists in the development of small scale structures via the diffusion of the magnetic fields. The reduced magnetic field strength that results in time permits shorter wavelength modes to grow. Considering fully varying viscosity and fully varying resistivity profiles in the simulation due to the strong dependence of viscosity and resistivity on the disparate temperature profile across the interface, the effect of both viscosity and resistivity is shown to be significant on the evolution of RTI in HED plasmas. Furthermore, it is also found that the dynamics of the magnetic field is explicitly independent of viscosity and likewise the resistivity does not affect the dynamics of enstrophy and kinetic energy. Also presented here is the power law scaling of enstrophy, kinetic energy, and magnetic field energy over a wide range of viscosity and resistivity in both injection range and inertial sub-range of spectra. This could provide a useful tool for understanding RTI induced turbulent mixing in high Atwood number HED plasmas and could aid in interpretation of observations of RTI-induced turbulence spectra.
† Email address for correspondence: bratan@vt.edu ‡ Email address for correspondence: srinbhu@vt.edu