© ESO, 2017

1. Introduction

The Rayleigh-Taylor Instability (RTI) is a dynamic mixing process which occurs when a lower density fluid pushes into a higher density fluid (Rayleigh 1882; Taylor 1950). This is usually realised by the lighter fluid supporting the heavier against gravity, and manifests as an interpenetration of finger-like plumes. The contact discontinuity between the two fluids is unstable to perturbations that grow by converting potential energy to kinetic energy, causing bubbles of the low-density fluid to rise, and spikes of the high-density fluid to sink. A thorough review of the hydrodynamic RTI is given by Sharp (1984). The mixing region¹ of a simulation of the RTI is shown in Fig. 1.

Fig. 1 The mixing region of Rayleigh-Taylor unstable fluids of density 1 and 10, where the pure fluids above and below the mixing region are not shown.

RTI is observed to play a role in many astrophysical systems, not least of all the Sun: the observed upflows in solar prominences have been confirmed to be RT unstable using 3D numerical magnetohydrodynamic (MHD) simulations (Hillier et al. 2012); the dynamics of back-falling plasma following a failed coronal mass ejection (CME) displayed morphology indicative of RTI, which was confirmed from Alfvén velocity estimates (Innes et al. 2012); further evidence of RTI was found following the same eruption at later times by Carlyle et al. (2014), who used evidence of the RTI to estimate magnetic field strength; RTI has been found to be a plausible mechanism for driving jets in supra arcade downflows (SADs) by Guo et al. (2014); the magnetic RTI has also been investigated in relation to filamentation of emerging flux (Isobe et al. 2005), and in the breaking up of magnetic flux sheets (Cattaneo & Hughes 1988). Further afield in the Universe, the structural formation of the crab nebula has been explained as an occurence of the RTI and a magnetic field strength was calculated by Hester et al. (1996), and more recently Porth et al. (2014) showed that the magnetic field strength in the crab nebula is not sufficient to suppress the instability, using high-resolution MHD simulations. The list goes on, with observations of the RTI reported on many scales, highlighting the importance of a thorough understanding of this process in an astrophysical context.

The magnetic RTI may be thought of as a competition between two forces: the gravitational potential, pulling the higher density fluid through the lower (and vice versa), mixing the two; and the magnetic tension, preventing deformation of the field lines and hence suppressing the mixing. Therefore it is useful to define some parameter relating these two forces so that the simulations may be compared with observations. This parameter was chosen as $\mathit{J}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}}{\mathit{g}\mathrm{\Lambda }}\mathit{,}$(1)where c_A is the Alfvén velocity, g is gravitational acceleration, and Λ is a characteristic length-scale. J is therefore a dimensionless parameter which describes the balance between magnetic and gravitational forces. A system with J ≫ 1 is likely to be dominated by the magnetic forces, and one with J ≪ 1 is likely to be dominated by the gravitational forces.

Using this J parameter, it is possible to compare astrophysical plasmas observed to undergo RTI with ideal MHD simulations. For erupted filamentary solar plasma, these values are approximately c_A = 5 × 10⁶ cm s^-1 and Λ = 10¹⁰ cm from Innes et al. (2012), and g = 10⁴ cm s^-2 by using the surface gravity of the Sun g_surf = 2.74 × 10⁴ cm s^-2 and taking into consideration that the plasma examined is apparently between 0.5−1 solar radii R_⊙. This gives J = 0.25. Alternately, from Hester et al. (1996) we obtain values of approximately c_A = 7 × 10⁶ cm s^-1, g = 3.5 × 10^-3 cm s^-2, and Λ = 10¹⁷ cm, which gives J = 0.14. The strongest magnetic fields studied in previous 3D numerical experiments appear in Stone & Gardiner (2007a) and are conducted in the J = 0.03 regime, so in order to better investigate the RTI in the context of astrophysical plasmas, a higher value of J should be explored.

Much insight can be gained into the RTI from analytic work carried out by Chandrasekhar (1961). If the fluids involved in the RTI are inviscid, the growth rate of the linear phase of the instability is described by the spatial frequency ω: ${\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathit{Akg,}$(2)where A is the Atwood number, defined as (ρ_u − ρ_l)/(ρ_u + ρ_l), ρ is density (and the subscript denotes upper and lower density fluids), k is wavenumber, and g is gravitational acceleration. The amplitude of disturbance to the boundary η a linear mode at a given time t in this case is defined by $\mathit{\eta }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{\eta }}_{\mathrm{0}}\exp \mathrm{\left[}\mathrm{\right(}\mathit{Akg}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{t}\mathrm{]}\mathit{,}$(3)where η₀ = η(t = 0) is the size of the initial small perturbation.

The addition of a magnetic field B parallel to the contact discontinuity (provided the fluids are sufficiently conductive) modifies this linear growth rate through the addition of magnetic tension along the field which will work to suppress high wavenumber perturbations. If the magnetic field is purely in the x-direction, the growth rate is given by ${\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathit{Akg}\mathrm{+}\frac{{\cos }^{\mathrm{2}}\mathit{\theta }{\mathit{k}}^{\mathrm{2}}{\mathit{B}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }\mathrm{(}{\mathit{\rho }}_{\mathrm{u}}\mathrm{+}{\mathit{\rho }}_{\mathrm{l}}\mathrm{)}}\mathit{,}$(4)(Chandrasekhar 1961) where θ is the angle between k and B. Much analytic work has been performed on this linear regime of the RTI, which has been used so frequently to attempt to explain observed astrophysical processes.

It is useful to estimate a characteristic length scale associated with the RTI, which can be achieved with these linear equations. If Eq. (4) is below 0, then ω is imaginary and the system is stable and any perturbations will produce waves in the interface. If ω is real, the system is unstable to that perturbation and the instability will give rise to the bubbles and spikes described. The most unstable wavelength of the instability is always the interchange mode, where θ = π/ 2, because it does no work against the magnetic field. However, the most unstable wavelength of the instability for a given θ, i.e., the characteristic lengthscale of the instability at a particular angle to the magnetic field, will be at the peak of the distribution of ω(k,θ), that is where ∂ω/∂k = 0. Doing this gives ${\mathit{\lambda }}_{\mathrm{u}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{{\mathit{k}}_{\mathrm{u}}}\mathrm{=}\frac{\mathrm{2}{\cos }^{\mathrm{2}}\mathit{\theta }{\mathit{B}}^{\mathrm{2}}}{\mathit{g}\mathrm{(}{\mathit{\rho }}_{\mathrm{u}}\mathrm{-}{\mathit{\rho }}_{\mathrm{l}}\mathrm{)}}\mathit{,}$(5)(e.g. Carlyle et al. 2014) where k_u is the wavenumber of the most unstable mode, and λ_u is the corresponding wavelength of this mode: the dominant growth scale of the instability.

Once the non-linear saturation of the instability has been reached, these equations will no longer describe the development of the system. There is not a definite distinction between the linear and non-linear regimes of the RTI, in fact a transitional stage between the two exists which is not easily described by simple ordinary differential equations. The argument can be made that the non-linear saturation is reached when the contact discontinuity has been deformed in the vertical direction over a distance comparable with 1 /k (Fermi & Von Neumann 1953). Since k = 2π/λ, qualitatively the non-linear saturation could be described as the point at which the vertical scales are comparable to the horizontal. See Sect. 5 of Hillier (2016) for a thorough description of this.

In the non-linear regime, the growth of the (hydrodynamic) RTI becomes self similar, and is described by $\frac{\mathit{\partial h}}{\mathit{\partial t}}\mathrm{=}\mathrm{2}\mathrm{(}\mathit{\alpha Agh}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(6)(Ristorcelli & Clark 2004; Cook et al. 2004) where h is the distance from the initial interface which the mixing region has penetrated the homogenous fluid (i.e. the height of the bubbles or the depth of the spikes), and t is time. The dimensionless coefficient α is referred to as the non-linear growth rate, and is insensitive to initial conditions. By taking the positive roots (i.e., those which are physically realisable), for constant α, A and g, the solution to Eq. (6) is $\mathit{h}\mathrm{=}\mathit{\alpha Ag}{\mathit{t}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{(}\mathit{\alpha Ag}{\mathit{h}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{t}\mathrm{+}{\mathit{h}}_{\mathrm{0}}\mathit{.}$(7)If t = 0 is chosen as the onset of non-linearity, then h₀ is the thickness of the mixing region at this point (the extent of the mixing region from the initial contact discontinuity; Cabot & Cook 2006). At later times, the first term of Eq. (7) dominates, and the latter two may be neglected, so we have $\mathit{h}\mathrm{\approx }\mathit{\alpha Ag}{\mathit{t}}^{\mathrm{2}}\mathit{.}$(8)For more in-depth analytic work into the magnetic RTI, the authors also recommend the text books by Goedbloed & Poedts (2004) and Goedbloed et al. (2010).

By comparing, it can be seen that the exponential growth of the linear regime will compete with the t² growth of the non-linear regime at smaller k, and so h should be greater than the maximum wavelength of the system in order for Eq. (8) to be applicable. It is also apparent that in a simulation, the characteristic domain size L should be greater than the characteristic separation of the most unstable mode λ_u, and as a good rule of thumb this is assumed to be of the same order as h. Therefore, non-linear analysis of the magnetic RTI should be conducted in a system which satisfies L>h>λ_max.

Non-linear growth rate α has been determined from multiple laboratory experiments to be approximately 0.057, however, studies of simulations of the RTI typically give a value a factor two smaller than this. Glimm et al. (2001) conclude that numerical dissipation effects (such as mass diffusion and viscosity) due to algorithmic differences and differences in simulation duration are the main reasons for the observed spread in non-linear growth rate across studies, and Dimonte et al. (2004) argue that the reduced growth rate in simulations arises from band-limited initial perturbations.

Jun et al. (1995) studied the linear and non-linear regimes of the RTI using 2D MHD simulations, investigating the effect of a magnetic field tangential to the initial interface as well as fields normal to this. They found enhanced growth (relative to the hydrodynamic case) in the normal case, the material collimating along field lines as the instability sets in, but there is an upper limit to the magnetic field strength, beyond which the growth is greatly suppressed. However, as pointed out by Hillier (2016), unlike the hydrodynamic case (or the case of Jun et al. 1995, where the magnetic fields are normal to the simulated plane), the evolution of the instability is no longer isotropic due to the addition of the magnetic field (tangential to the interface), and so a 2D simulation would only be able to capture the growth of a single mode from the whole spectrum of preferred modes; a fundamentally 3D system cannot be truly captured by 2D simulations.

The non-linear phase of the RTI with a magnetic field has been studied in 3D MHD simulations: Stone & Gardiner (2007a) showed how the shape of resulting bubbles is affected by magnetic field configuration, and Stone & Gardiner (2007b) demonstrated that whilst the instability was slowed by the addition of a strong magnetic field during the initial onset of the instability, at later times the non-linear growth rate was actually enhanced relative to the hydrodynamic case. This is attributed to the suppression of mixing between the fluids, which occurs through secondary Kelvin-Helmholtz roll-ups at the edges of the bubbles and fingers, preserving the density discontinuity. This directly refutes the postulation above that unidirectional magnetic fields suppress the modes of the instability in one dimension, reducing the overall growth rate of the system. Beyond these investigations, little numerical work has been carried out on the non-linear saturation of the RTI.

This study aims to conduct idealised 3D numerical MHD simulations in a parameter space approximating astrophysical plasma, and the details of these are given. It should be noted that whilst we approach realistic values for the J parameter (Eq. (1)), the simulations are still highly idealised; for example, magnetic field may not always be aligned with the interface in reality (nor indeed will the magnetic fields in each plasma be similarly oriented, or even in strength, necessarily). This work also neglects other physical processes, such as radiative transfer and ionisation balance, which are not thought to have appreciable effects on the growth of plasma instabilities under varying magnetic field. We then present analysis of the non-linear regime of the instability, particularly the growth rate α and the mixing of the system. Finally, we discuss the implications and validity of these results with respect to application to the observations which inspired this study.

2. Numerical MHD simulations

This work used the Athena code for astrophysical MHD (see Stone et al. 2008, for a complete description of this code), which solves the equations of ideal MHD with a constant gravitational acceleration, g = (0,0,−g): $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\\ & & \frac{\mathit{\partial \rho }{v}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathit{\rho }{v}{v}\mathrm{-}{BB}\mathrm{)}\mathrm{+}\mathrm{\nabla }\mathit{P}\mathrm{=}\mathit{\rho }{g}\\ & & \frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{\times }\mathrm{(}{v}\mathrm{\times }{B}\mathrm{)}\mathrm{=}\mathrm{0}\\ & & \frac{\mathit{\partial E}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left[}\mathrm{\right(}\mathit{E}\mathrm{+}\mathit{P}\mathrm{)}{v}\mathrm{-}{B}\mathrm{(}{B}\mathrm{·}{v}\mathrm{\left)}\mathrm{\right]}\mathrm{=}\mathit{\rho }{v}\mathrm{·}{g}\mathit{,}\end{array}$where we have used total pressure P ≡ P_g + (B·B)/2, gas pressure P_g = (γ−1)ϵ, total energy density E ≡ ϵ + ρ(v·v/ 2) + (B·B)/2, internal energy ϵ, and the adiabatic index γ = 5/3. This is not the value which would necessarily be expected from e.g. solar plasma, however, the simulations are conducted at the incompressible limit by using a large enough sound speed such that all fluid motions are highly subsonic, and so varying the adiabatic index has little effect on the results (Stone & Gardiner 2007b). Note that these equations have been normalised to dimensionless units such that sound speed c_s = 1 (at the interface between the fluids) for B = 1 and ρ = ρ_u = 1, and the characteristic length scale of the system Λ = 1. In this model, g = 0.1 and so $\sqrt{\mathit{g}\mathrm{\Lambda }}\mathit{/}{\mathit{c}}_{\mathrm{s}}\mathrm{\ll }\mathrm{1}$, which indicates that the induced flows are almost incompressible.

The equations are solved using a second-order Godunov scheme. Perhaps the most important element of this scheme is the Riemann solver, which calculates time-averaged fluxes of all conserved quantities at cell interfaces. Here a multi-state Harten-Lax-van Leer Discontinuities (HLLD) approximate Riemann solver is used, since this is as accurate as the well studied Roe approximate Riemann solver and less computationally demanding (Miyoshi & Kusano 2005). This is combined with the constrained transport (CT) technique which evolves the induction equation in a way which ensures zero divergence of the poloidal (constrained) field components to within machine round-off error (Evans & Hawley 1988). Discretization is based on cell-centered volume averages for mass, momentum, and energy, and face-centered area averages for the magnetic field. Athena has been shown to be successful at conducting MHD simulations of the RTI in three dimensions (Stone & Gardiner 2007a) into the non-linear regime (Stone & Gardiner 2007b), and as such it was deemed suitable for conducting this investigation.

The x and y boundaries of the domain are periodic whilst the z boundaries are reflective, and the origin is in the centre of the domain. The regular cartesian grid used has dimensions of 256 × 64 × 1024. Resolution in the z-dimension (height) was doubled relative to x and y – that is, dz = dx/ 2 = dy/ 2 – so as to achieve a high precision and accuracy of measurement of height and therefore growth rate.

The system is initially in hydrostatic equilibrium, and the gas pressure is chosen such that the sound speed (c_s) in the light fluid at the interface is unity, and so $\mathit{P}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{5}}\mathrm{-}\mathit{g\rho z}\mathrm{+}\frac{{\mathit{B}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{·}$(13)The characteristic length scale Λ is roughly an order of magnitude larger than the scales predicted by Eq. (5) for desired magnetic field strengths. The width of the domain in the direction of magnetic field used for the first set of simulations to be L_x = 0.4Λ. This width is chosen to allow L_x ≥ λ_u for all simulations, and resolves the dominant wavelengths λ_u with at least 44 grid points.

RT modes perpendicular to the magnetic field behave as the hydrodynamic case, and so the smallest scales are favoured. Numerical diffusion in the simulations was of the order 0.01Λ for the resolution used, which is measured from the extent over which the contact discontinuity reaches after the two pure fluids are allowed to relax with no perturbations given to the system. Therefore, L_y = 0.1Λ is used as the depth of the domain, allowing sufficient space for interchange structures to develop. A height of L_z = 0.8Λ is used to ensure the L_z ≫ h is not violated, and to prevent the growth of the bubbles being affected by the reflective upper boundary, as it is for this reason that Stone & Gardiner (2007a) discarded 20% of their data.

The lowest J (corresponding to the weakest magnetic field strength, see Eq. (1)) corresponds to the Athena RTI test case (and as such has been rigorously analysed and tested for accuracy), however, higher J simulations have not previously been conducted; the highest J (strongest field) used here is at the limit of L ≃ λ_u, violating the requirement of L>h>λ_max (this was further investigated in the second set of simulations, which are described below). A larger L was not used as the simulations were already computationally demanding; a lower resolution was also avoided as the current setup should lead to approximately 50 pixels per λ_u, and lower resolution is not desirable as it is important that the simulation allows all scales dictated by the physics to develop, and not be inhibited for computational reasons. The magnetic field is initially applied uniformly along the x axis, that is (B_x = const., B_y = 0, B_z = 0). Seven simulations were run in this set, and are described in Table 1.

The mixing layer (that is all fluid with 1.5 ≤ ρ ≤ 9.5, where the initial setup has ρ_l = 1 for z< 0 and ρ_u = 10 for z ≥ 0) of B1, B3, B5 and B7 are shown at three points along the run in Fig. 2. The chosen start time of 0.1 rather than 0 is to show the interface; at t = 0, the lower half of the domain is filled with ρ_l = 1 material and the upper half with ρ_u = 10, so no mixing layer is visible. As the simulations progress, bubbles of scales predicted by Eq. (5) can be seen developing along x, the direction along which B is directed. The scales of these so-called undular modes are seen to increase as J (and hence magnetic field strength) increases, whilst the scales across the magnetic field, the interchange mode, remain apparently constant for all simulations: one bubble along the y-direction can be seen.

Fig. 2 Visualisations of the mixing region 1.5 ≤ ρ ≤ 9.5 of simulations B1, B3, B5 and B7 (cf. Table 1) at t = 2, 4, 6.

The structures created in the linear and non-linear phases of the instability, as shown in Fig. 2, are clearly dependent on the existence and strength of the magnetic field. Generally the instability drives the creation of filamentary structure that is aligned with the direction of the magnetic field. Looking at t = 2 for the B3 simulation, approximately five peaks can be seen across the length of the box, but if these were just undular modes we would expect to see about three (see Table 1). This implies that the formation of mixed modes, with structure across the magnetic field, is playing an important role. As the instability develops in its non-linear phase, larger structures develop both across and along the magnetic field, where the simulations with stronger magnetic fields maintain their larger aspect ratio between the along field and across field scales even into this regime. These impressions from the figure have been confirmed through a Fourier transform of the data, presented in Fig. 3. This shows that the power in higher frequencies is reduced for stronger magnetic field only in the direction of the magnetic field, indicating that smaller scales are suppressed by magnetic field.

Fig. 3 Fourier transform of Bz at t = 1 (a, d), t = 2 (b, e), and t = 6 (c, f), of the central z-slice for selected B simulations, showing the spectral power of different spatial frequencies (a–c) shows the scales aligned with magnetic field (d–f) shows scales across.

Since the only two parameters which vary between B1–B7 are J and L_x/λ_u (the ratio of domain width to dominant linear wavelength of the RTI), it is necessary to ensure the observed change in growth rate is due to the former rather than the latter. In order to achieve this, a second set of simulations is run, where magnetic field is kept constant but the width of the domain is made progressively smaller. If a dependence of α on L/λ_u can be seen in simulations of constant magnetic field similar to this dependence in the previous set of simulations, then the magnetic field can not be said to be causing the postulated effect on growth rate. Another benefit of this is to investigate behaviour which violates the L_x>h>λ_max requirement for a constant magnetic field. Four simulations are conducted with J = 0.03 and other parameters detailed in Table 1. Figure 4 shows the mixing region for W2, 3 and 4 for t = 2, 4, 6. Note that W1 has the same initial conditions as B1.

Fig. 4 Visualisations of the mixing region 1.5 ≤ ρ ≤ 9.5 of simulations W2, W3 and W4 (cf. Table 1) at t = 2, 4, 6.

3. Analysis

The measure of bubble height h is taken as being the highest point at which the average density over the x−y-plane is ⟨ ρ_z ⟩ ≤ 9.5, which should return a position at the average height of all bubbles; Figs. 2 and 4 show the sharp density gradients at the edges of the mixing region. Figure 5 shows the development of the highest ⟨ ρ_z ⟩ = 9.5 (i.e., Eq. (8)) for each simulation in the magnetic field varying set. The gradient of the slopes defines the relative rate at which bubbles grow, and it is apparent that this decreases across the simulations from B1 to B7 from visual inspection, suggesting that increased magnetic field strength will yield a reduced (non-linear) growth rate (as well as agreeing with the analytic prediction that linear growth rate decreases with magnetic field strength; see Eq. (4)).

Fig. 5 Development of bubble height as a function of Agt2 for the B simulations (increasing magnetic field strength). Dotted lines mark t = 2 and t = 4; simulations end at t = 6. Dashed line shows the slope of α = 0.04.

The early linear phase of the RTI can be seen in Fig. 5, characterised by a rapid growth. The rate of growth appears to then suddenly decrease at the same point in all simulations, continuing thereon with a relatively steady dependence on t². Some of the B simulations diverge from this dependence towards later times, but this is likely to be the result of the formation of large, coherent flows developing which weaken the statistics of the averaging process.

The non-linear growth rate is calculated by finding the rate of change of h relative to Agt² for each time-step (by fitting a linear regression to the surrounding 100 datapoints) in the non-linear regime, and taking the mean value. This gives a value for α, as well as the standard deviation, σ, for each simulation. These are given in Table 1, and are plotted against J for all B simulations in Fig. 6, with error bars representing the standard deviation.

Fig. 6 α plotted against J for all B simulations – note the scaling of the axes is logarithmic.

A mixing parameter may also be defined as Θ = 4 ⟨ f_hf_l ⟩ (Stone & Gardiner 2007a), where f denotes the fractional amount of upper and lower density material, respectively. That is f_l = | 10 − ρ | /9 and f_h = | 1 − ρ | /9, where these are averaged over the x−y-plane. These data are plotted for all B simulations at t = 4 in Fig. 7.

Fig. 7 The average degree of mixing over the height of all B simulations at t = 4.

Similar plots of the W simulations are shown in Fig. 8, which also appear to show a change in α across the simulations, listed in Table 1, where L/λ_u is the only parameter which is decreasing between simulations. This suggests that magnetic field is not necessarily having a strong effect on the non-linear growth of the RTI, though the trend in α is much clearer and more correlative in the initial (B) set. α is also found for the W simulations over the non-linear phase; these values are listed in Table 1.

Fig. 8 Development of bubble height as a function of Agt2 for the W simulations (decreasing domain width). Dashed lines mark t = 2 and t = 4; simulations end at t = 6.

Fig. 9 Panel a: contour plot of the horizontal mean of the density for B1. Panel b: same for B8. Panel c: same for W4. The dashed white and black lines signify α = −0.047 and α = 0.047, respectively, for panel a, α = −0.055 and α = 0.0354 for panel b, and α = −0.044 and α = 0.038 for panel c. Panel d: mean density distribution plotted against z/Agt2 for the final 200 snapshots of the simulation for B1. Panel e: same but for B8. Panel f: same but for W4. Note that the values that would be calculated for α from this figure are over estimates due to the initial expansion of the system.

Figure 9 shows the expansion of the mixing layer for simulations B1 and B8. The top three panels are contour plots of the ⟨ ρ ⟩ for B1, B8 and W4 respectively with z as the horizontal axis and Agt² for the vertical. The dashed white and black lines highlight the lower and upper limits of the expansion of the mixing with t² dependence. The values for α associated with this are −0.047 and 0.047 for B1, −0.055 and 0.0354 for B8, and α = −0.044 and α = 0.0378 for W4. The bottom three panels of this figure show the average of the mean density profiles over the final 200 snapshots of the simulation for B1 (panel d), B8 (panel e) and W4 (panel f). Note that the values for α that can be calculated from this plot are overestimates due to the fast initial expansion of the layer; this issue is discussed in Cabot & Cook (2006).

4. Discussion

Simulations of RT-unstable plasma (governed by an ideal equation-of-state) are conducted close to the incompressible limit, with homogenous magnetic fields aligned perpendicular to gravity (parallel to the density-jump interface in the initial conditions), and the impact of the strength of the embedded magnetic fields on the development of the non-linear regime of the RTI is examined. From visual inspection of Figs. 5 and 6, and from the calculated α in Table 1, it is apparent that increased magnetic field strength leads to reduced non-linear growth rate. However, it is not possible to draw a linear regression trend through all data in Fig. 6; that is, the decrease in non-linear growth rate as magnetic field strength increases is not smooth. Figure 6 (as well as Table 1) shows, for example, that the growth rate for the B2 simulation is much lower than that of B3 and even B4, each of which have progressively stronger magnetic fields.

Some values of α returned in this work are not necessarily equal to that of other studies (in fact a single value of α is not agreed upon between very many studies). This may be explained by Glimm et al. (2001), who show that numerical dissipation effects (such as mass diffusion and viscosity) due to algorithmic differences and differences in simulation duration are the main reasons for this discrepancy across studies, whilst within this study these effects are constant across simulations and so the calculated non-linear growth rates are suitable to be compared to one another, though comparisons of absolute values with other work is less trustworthy.

The reduced growth rate for stronger fields could be due to an increased magnetic tension (= (B·∇)B/μ₀) in stronger magnetic fields. Greater magnetic tension reduces the free energy of the system along the direction of the magnetic field line, suppressing the modes of the instability in this dimension, and with the reduced contribution to the whole system from these modes, the overall growth-rate will be reduced. This is not in disagreement with Stone & Gardiner (2007a), who conclude that the RTI non-linear growth rate is faster when a magnetic field is added relative to the hydrodynamic case, explained by the reduced interface-mixing from suppressed Kelvin-Helmholtz roll-ups by magnetic tension.

Whilst magnetic field strength is the only parameter altered between successive runs of the initial B set, this leads to a secondary constraint on the physics of the system: the ratio of the dominant scale of the instability (λ_u, from Eq. (5)) to domain width (L_x) may constrain the bubbles which develop along the direction of magnetic field. If this is below unity, Eq. (5) gives a characteristic scale size of the simulation as larger than the simulated domain. Moreover, as explained in Sect. 1, in order for Eq. (8) to apply to the system (i.e. the equation relating h and α), we must fulfil L>h>λ_max. In order to investigate whether λ_u/L affects the growth rate of the (either linear or non-linear) RTI, as well as to examine the behaviour of a system which violates the L>h>λ_max constraint, a further set of simulations were run with constant magnetic field strength but variable width. These displayed non-linear growth rates which showed similar variation to the first B set, but also with no apparent correlation and greater uncertainty. This is highlighted by comparing Figs. 5 and 8, as well as the standard deviations in Table 1. The height vs. time-squared plots are visually more correlative for the B set of simulations than the W set, and the standard deviation on α is somewhat higher for the W simulations (see Table 1). The lack of correlation in α for the W simulations gives validity to the results of the B simulations which do appear to display a trend, indicating that enhanced magnetic field strength leads to reduced non-linear growth-rates.

Measured non-linear growth rate α for the W simulations also highlighted the importance of the L>h>λ_max requirement. It is apparent that where this is violated (W4) the RTI does not develop in the same manner as all other simulations conducted – from visual examination of Fig. 8, the slope of W4 mixing region height does not follow the same trend as other simulations plotted in this figure, nor in any from Fig. 5. This can also be seen in Table 1, which shows the standard deviation on W4’s calculated α is the highest of all simulations.

There are further implications for the 3D RTI with unidirectional embedded magnetic fields; if the magnetic tension opposes the deformation of the field lines induced by plasma motions, but only in the direction of the field lines, this is essentially reducing the free energy of the system in only this dimension. It would therefore follow that by introducing progressively stronger magnetic fields, the system approaches a quasi-2D domain, and hence the growth rate of the instability would be reduced. This is demonstrated by Kane et al. (2000), who show that α is reduced by roughly 30% in the 2D case relative to the 3D case. This can explain the result obtained here for reduced α with increased magnetic field strength.

In order to measure the non-linear growth rate of the RTI, the gradient of a curve such as those plotted in Fig. 5 is often commonly used. However, a precise value of α for each simulation is difficult to measure, as this is an attempt to describe the average behaviour of the non-linear system; the non-linearity itself implies fluctuations which will change α on small timescales. For this reason, the mean value of the rate of change of the height of the mixing region over the non-linear regime was used to calculate alpha.

Furthermore, the transition between linear and non-linear regimes of the instability is ill-defined, so the α measurement starting point can be difficult to choose. As a good approximation, the eigenfunction for the vertical velocity v_z is given by ${\mathit{v}}_{\mathit{z}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathit{A}{\mathrm{e}}^{\mathrm{-}\mathit{k}\mathrm{\left|}\mathit{z}\mathrm{\right|}}$(14)(Chandrasekhar 1961), which implies that 1 /k can be used as the vertical scale through which the perturbation can travel before it reaches its non-linear saturation. Figures 2 and 4 indicate that this is achieved just after t = 2 for the low-strength magnetic field cases, and at later times as this is increased. It is possible that non-linear saturation has not been reached for the strongest magnetic fields by the end of the simulation; the bottom-right of Fig. 2 suggests that the height of the mixing region has not yet reached the observed wavelength. However, from Fig. 5, B8 indeed displays the shift of behaviour into an apparent t² dependence at roughly the same time as all other simulations, suggesting that a vertical scale which depends on the magnetic field strength (indirectly, i.e., through modification of 1 /k) may not be the best method of identifying non-linear saturation. Note it has been shown that the initial expansion of the system when dominated by linear modes can also display the same t² dependence (Hillier 2016), which is consistent with our findings in this paper.

The results for the mixing parameter in Fig. 7 show that increasing the magnetic field does not seem to drastically increase the density mixing of the system, which is in contrast to the results found for including shear in the magnetic field of Stone & Gardiner (2007a).

The non-linear growth rate is also calculated from the upper and lower limits of the expansion of the mixing layer (see Fig. 9). This finds that for stronger magnetic field cases, there is an asymmetry between the upwards growth of the bubbles and the downwards growth of the spikes. The difference between the two α values is approximately 0.09 for both B1 and B8, so though B1 grows quicker upwards, the overall expansion of the layers in both simulations is approximately the same. W4, however, has a smaller difference of approximately 0.08 which we can attribute to the simulation becoming more 2D like which is know to reduce the value of α. This leads to the conclusion that the development of asymmetry from B1 to B8 is to some extent a magnetic effect and not purely as a result of the stronger magnetic field rendering the B8 simulation more 2D like. A note of caution must be presented with this, because W4 by necessity has poorer statistics for the averaging (due to less structure being averaged across resulting in individual plumes carrying greater weight in the averaging process), and so it is harder to draw strong conclusions.

From panels (d) and (e) in Fig. 9, it can be seen that the asymmetry in the expansion of the mixing layer of B8 also manifests itself in the distribution of the mean density of the mixing layer. For B1, the density profile is approximately linear, but for B8 a long tail exists in the distribution at the lower level. The cause of this asymmetry is unclear, however, it does suggest that the reduced growth rate for enhanced magnetic field strengths may not be an accurate interpretation of these results. This asymmetry reflects a reduction of the erosion of the upper dense layer (as evidenced by the lower α value). As the erosion of this layer reflects the release of the gravitational potential energy that is driving the mixing, it is no surprise that on calculating the energy release, given by ${\mathit{E}}_{\mathrm{GP}}\mathrm{=}{\mathrm{\int }}_{-0.1}^{\mathrm{0.1}}\mathrm{(}\mathrm{⟨}\mathit{\rho }\mathrm{(}{\mathit{z}}^{\mathrm{\prime }}\mathrm{)}\mathrm{⟩}\mathrm{-}{\mathit{\rho }}_{\mathrm{init}}\mathrm{(}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right)}\mathit{g}{\mathit{z}}^{\mathrm{\prime }}\mathrm{d}\mathit{z}\prime$(15)where z′ = z/Agt² and ρ_init is the initial density distribution, the energy released in simulation B8 is only 76% of that of B1. As energy is a conserved quantity, a reduction in the amount of gravitational energy released must necessarily result in reductions in the amount of energy found in the turbulent components of the velocity and magnetic field, because there is less energy released to drive the growth of these components. Therefore, the fact that the increase in magnetic field strength reduces the gravitational energy release leads to a reduction in the turbulence driven by this instability.

The effect of J (see Eq. (1)) is explored over a parameter space ranging from the highest value conducted in previous work (Stone & Gardiner 2007a), and a maximum value which is an accurate representation of plasma observed to be RT unstable in the crab nebula (Hester et al. 1996), and is factor two lower than that predicted for erupted solar filament plasma such as the study by Innes et al. (2012). The simulations are conducted on plasma at the incompressible limit, and in both of these examples mentioned, the sound speed is believed to be greater than velocities of any fluid motions. Therefore the results of this study are applicable to the observations of astrophysical plasmas which initially prompted this investigation. In the initial conditions of the simulations, a homogenous magnetic field parallel to the interface is used. Such a configuration seems unlikely to exist in nature, however, it can be very difficult to observe the magnetic structure within astrophysical plasmas. Since no previous studies have investigated the effect of the strength of the magnetic field on the non-linear RTI, this is the most straightforward case to begin investigating. The results demonstrate that the strength of magnetic fields embedded in plasmas will affect the development of the non-linear mixing by the RTI.

5. Conclusion

Simulations of the RTI were conducted in order to better understand how magnetic field strength may affect the growth rate of this instability. It has been found from previous work that non-linear growth rate is enhanced when a strong magnetic field is present (cf. the hydrodynamic case), however, this study has found that increasing the strength of the magnetic field leads to a decrease in non-linear growth rate of the rising bubbles of the instability. This is speculated to be due to higher magnetic tension requiring greater energy in order for the frozen-in plasma to move. Non-linear growth rates were found to converge on ~0.039 for the strongest magnetic fields studied, however, the decrease in non-linear growth rate in the results presented in this work is by no means a smooth one; for instance, the growth rate of the second-weakest field is almost equal to that of the second-strongest field.

Since altering magnetic field strength in a fixed domain has the indirect effect of changing the ratio between domain size and dominant mode (i.e. characteristic length-scale) of the RTI, this was investigated for a set of simulations with constant magnetic field strength. Whilst the calculated growth rates for these simulations also showed some difference, the results were not correlative at all, as the first set was, implying that the interpretation of stronger magnetic fields reducing the non-linear development of the RTI is accurate.

Finally, the non-linear growth rate was estimated in a different way from both the rising bubbles and the falling spikes; these results corroborated the earlier conclusion that stronger fields yield lower growth rates in the rising bubbles, however, this also indicated that the growth rate of the falling bubbles is enhanced by approximately the same order. This suggests that stronger fields do not enhance the development of the RTI, but in fact create an asymmetry. This leaves the α value as a measure of the rate at which the dense layer is eroded, and its energy is released, where we find that stronger magnetic fields slow down the release of this energy.

Whilst the simulations may not perfectly describe observed RT-unstable plasma throughout the universe (e.g., the simulations are conducted approximately in the incompressible limit, with homogenous magnetic fields, under the ideal gas equation), they are conducted in a parameter space relevant to the relative effects of magnetic tension and gravity on such astrophysical plasmas. Hence, we have shown that increased magnetic field strength reduces the non-linear development of the RTI in situations where the instability is observed throughout the universe. This arises from the magnetic tension which acts against the deformation of field lines, and hence acting against the instability in the dimension aligned with the field. By reducing the growth in this way (i.e., in one dimension), the development of the instability overall in the system is reduced.

Acknowledgments

The authors wish to thank Davina Innes for inspiring this project, as well as invaluable discussions and advice; and thanks to LiJia Guo for assistance with performing the numerical experiments. Jack Carlyle conducted this work during a Ph.D. joint-funded by University College London and The Max Planck Institute for Solar System Research, as well as whilst receiving funding as a Research Fellow from the European Space Agency. Andrew Hillier is supported by his STFC Emest Rutherford Fellowship grant number ST/L00397X/2.

References

All Tables

All Figures

	Fig. 1 The mixing region of Rayleigh-Taylor unstable fluids of density 1 and 10, where the pure fluids above and below the mixing region are not shown.
In the text

	Fig. 2 Visualisations of the mixing region 1.5 ≤ ρ ≤ 9.5 of simulations B1, B3, B5 and B7 (cf. Table 1) at t = 2, 4, 6.
In the text

	Fig. 3 Fourier transform of Bz at t = 1 (a, d), t = 2 (b, e), and t = 6 (c, f), of the central z-slice for selected B simulations, showing the spectral power of different spatial frequencies (a–c) shows the scales aligned with magnetic field (d–f) shows scales across.
In the text

	Fig. 4 Visualisations of the mixing region 1.5 ≤ ρ ≤ 9.5 of simulations W2, W3 and W4 (cf. Table 1) at t = 2, 4, 6.
In the text

	Fig. 5 Development of bubble height as a function of Agt2 for the B simulations (increasing magnetic field strength). Dotted lines mark t = 2 and t = 4; simulations end at t = 6. Dashed line shows the slope of α = 0.04.
In the text

	Fig. 6 α plotted against J for all B simulations – note the scaling of the axes is logarithmic.
In the text

	Fig. 7 The average degree of mixing over the height of all B simulations at t = 4.
In the text

	Fig. 8 Development of bubble height as a function of Agt2 for the W simulations (decreasing domain width). Dashed lines mark t = 2 and t = 4; simulations end at t = 6.
In the text

Fig. 9 Panel a: contour plot of the horizontal mean of the density for B1. Panel b: same for B8. Panel c: same for W4. The dashed white and black lines signify α = −0.047 and α = 0.047, respectively, for panel a, α = −0.055 and α = 0.0354 for panel b, and α = −0.044 and α = 0.038 for panel c. Panel d: mean density distribution plotted against z/Agt2 for the final 200 snapshots of the simulation for B1. Panel e: same but for B8. Panel f: same but for W4. Note that the values that would be calculated for α from this figure are over estimates due to the initial expansion of the system.

In the text