The nonlinear growth of the magnetic RayleighTaylor instability
^{1} European Space Agency, ESTEC, 2201, Noordwijk, The Netherlands
email: jack.carlyle@esa.int
^{2} College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter EX4, UK
Received: 17 March 2017
Accepted: 24 July 2017
This work examines the effect of the embedded magnetic field strength on the nonlinear development of the magnetic RayleighTaylor instability (RTI) (with a fieldaligned interface) in an ideal gas close to the incompressible limit in three dimensions. Numerical experiments are conducted in a domain sufficiently large so as to allow the predicted critical modes to develop in a physically realistic manner. The ratio between gravity, which drives the instability in this case (as well as in several of the corresponding observations), and magnetic field strength is taken up to a ratio which accurately reflects that of observed astrophysical plasma, in order to allow comparison between the results of the simulations and the observational data which served as inspiration for this work. This study finds reduced nonlinear growth of the rising bubbles of the RTI for stronger magnetic fields, and that this is directly due to the change in magnetic field strength, rather than the indirect effect of altering characteristic length scales with respect to domain size. By examining the growth of the falling spikes, the growth rate appears to be enhanced for the strongest magnetic field strengths, suggesting that rather than affecting the development of the system as a whole, increased magnetic field strengths in fact introduce an asymmetry to the system. Further investigation of this effect also revealed that the greater this asymmetry, the less efficiently the gravitational energy is released. By better understanding the understudied regime of such a major phenomenon in astrophysics, deeper explanations for observations may be sought, and this work illustrates that the strength of magnetic fields in astrophysical plasmas influences observed RTI in subtle and complex ways.
Key words: instabilities / magnetohydrodynamics (MHD) / plasmas / magnetic fields
© ESO, 2017
1. Introduction
The RayleighTaylor 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 fingerlike 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 lowdensity fluid to rise, and spikes of the highdensity fluid to sink. A thorough review of the hydrodynamic RTI is given by Sharp (1984). The mixing region^{1} of a simulation of the RTI is shown in Fig. 1.
Fig. 1 The mixing region of RayleighTaylor unstable fluids of density 1 and 10, where the pure fluids above and below the mixing region are not shown. 

Open with DEXTER 
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 backfalling 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 highresolution 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 (1)where c_{A} is the Alfvén velocity, g is gravitational acceleration, and Λ is a characteristic lengthscale. 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^{6} cm s^{1} and Λ = 10^{10} cm from Innes et al. (2012), and g = 10^{4} cm s^{2} by using the surface gravity of the Sun g_{surf} = 2.74 × 10^{4} 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^{6} cm s^{1}, g = 3.5 × 10^{3} cm s^{2}, and Λ = 10^{17} 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 ω: (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 (3)where η_{0} = η(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 xdirection, the growth rate is given by (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 (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 nonlinear 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 nonlinear 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 nonlinear 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 nonlinear 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 nonlinear regime, the growth of the (hydrodynamic) RTI becomes self similar, and is described by (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 nonlinear 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 (7)If t = 0 is chosen as the onset of nonlinearity, then h_{0} 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 (8)For more indepth 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^{2} growth of the nonlinear 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, nonlinear analysis of the magnetic RTI should be conducted in a system which satisfies L>h>λ_{max}.
Nonlinear 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 nonlinear growth rate across studies, and Dimonte et al. (2004) argue that the reduced growth rate in simulations arises from bandlimited initial perturbations.
Jun et al. (1995) studied the linear and nonlinear 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 nonlinear 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 nonlinear 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 KelvinHelmholtz rollups 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 nonlinear 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 nonlinear 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): 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 , which indicates that the induced flows are almost incompressible.
The equations are solved using a secondorder Godunov scheme. Perhaps the most important element of this scheme is the Riemann solver, which calculates timeaveraged fluxes of all conserved quantities at cell interfaces. Here a multistate HartenLaxvan 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 roundoff error (Evans & Hawley 1988). Discretization is based on cellcentered volume averages for mass, momentum, and energy, and facecentered 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 nonlinear 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 zdimension (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 (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.
Relevant parameter space explored of all simulations conducted, alongside measured nonlinear growth rate α and standard deviation σ.
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 socalled 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 ydirection 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. 

Open with DEXTER 
The structures created in the linear and nonlinear 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 nonlinear 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 B_{z} at t = 1 (a, d), t = 2 (b, e), and t = 6 (c, f), of the central zslice 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. 

Open with DEXTER 
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. 

Open with DEXTER 
3. Analysis
The measure of bubble height h is taken as being the highest point at which the average density over the x−yplane 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 (nonlinear) 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 Agt^{2} 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. 

Open with DEXTER 
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^{2}. 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 nonlinear growth rate is calculated by finding the rate of change of h relative to Agt^{2} for each timestep (by fitting a linear regression to the surrounding 100 datapoints) in the nonlinear 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. 

Open with DEXTER 
A mixing parameter may also be defined as Θ = 4 ⟨ f_{h}f_{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−yplane. 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. 

Open with DEXTER 
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 nonlinear 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 nonlinear phase; these values are listed in Table 1.
Fig. 8 Development of bubble height as a function of Agt^{2} for the W simulations (decreasing domain width). Dashed lines mark t = 2 and t = 4; simulations end at t = 6. 

Open with DEXTER 
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/Agt^{2} 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. 

Open with DEXTER 
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^{2} for the vertical. The dashed white and black lines highlight the lower and upper limits of the expansion of the mixing with t^{2} 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 RTunstable plasma (governed by an ideal equationofstate) are conducted close to the incompressible limit, with homogenous magnetic fields aligned perpendicular to gravity (parallel to the densityjump interface in the initial conditions), and the impact of the strength of the embedded magnetic fields on the development of the nonlinear 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 nonlinear growth rate. However, it is not possible to draw a linear regression trend through all data in Fig. 6; that is, the decrease in nonlinear 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 nonlinear 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/μ_{0}) 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 growthrate will be reduced. This is not in disagreement with Stone & Gardiner (2007a), who conclude that the RTI nonlinear growth rate is faster when a magnetic field is added relative to the hydrodynamic case, explained by the reduced interfacemixing from suppressed KelvinHelmholtz rollups 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 nonlinear) 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 nonlinear 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. timesquared 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 nonlinear growthrates.
Measured nonlinear 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 quasi2D 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 nonlinear 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 nonlinear system; the nonlinearity 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 nonlinear regime was used to calculate alpha.
Furthermore, the transition between linear and nonlinear regimes of the instability is illdefined, 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 (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 nonlinear saturation. Figures 2 and 4 indicate that this is achieved just after t = 2 for the lowstrength magnetic field cases, and at later times as this is increased. It is possible that nonlinear saturation has not been reached for the strongest magnetic fields by the end of the simulation; the bottomright 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^{2} 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 nonlinear saturation. Note it has been shown that the initial expansion of the system when dominated by linear modes can also display the same t^{2} 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 nonlinear 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 (15)where z′ = z/Agt^{2} 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 nonlinear 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 nonlinear 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 nonlinear 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 nonlinear 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 frozenin plasma to move. Nonlinear growth rates were found to converge on ~0.039 for the strongest magnetic fields studied, however, the decrease in nonlinear growth rate in the results presented in this work is by no means a smooth one; for instance, the growth rate of the secondweakest field is almost equal to that of the secondstrongest 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 lengthscale) 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 nonlinear development of the RTI is accurate.
Finally, the nonlinear 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 RTunstable 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 nonlinear 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.
In this paper, the word “mixing” is used to refer to the region in height of RTunstable material that contains a combination of both the higher and lowerdensity fluids. Whilst this is not necessarily a new, mixed phase fluid (i.e., another common definition of “mixing”), the initial regions of pure high and lowdensity plasma do indeed mix together in the RTI and so in the interests of eloquence, this massredistributionlayer is henceforth referred to as the mixing region.
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. jointfunded 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
 Cabot, W. H., & Cook, A. W. 2006, Nat. Phys., 2, 562 [CrossRef] [Google Scholar]
 Carlyle, J., Williams, D. R., van DrielGesztelyi, L., et al. 2014, ApJ, 782, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Cattaneo, F., & Hughes, D. W. 1988, J. Fluid Mech., 196, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford: Clarendon) [Google Scholar]
 Cook, A. W., Cabot, W., & Miller, P. L. 2004, J. Fluid Mech., 511, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Dimonte, G., Youngs, D. L., Dimits, A., et al. 2004, Phys. Fluids, 16, 1668 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [NASA ADS] [CrossRef] [Google Scholar]
 Fermi, E., & Von Neumann, J. 1953, Taylor instability of incompressible liquids, AECU; 2979 (Oak Ridge, Tenn.: United States Atomic Energy Commission, Technical Information Service) [Google Scholar]
 Glimm, J., Grove, J. W., Li, X. L., Oh, W., & Sharp, D. H. 2001, J. Comput. Phys., 169, 652 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Goedbloed, J. P. H., & Poedts, S. 2004, Principles of Magnetohydrodynamics (UK: Cambridge University Press) [Google Scholar]
 Goedbloed, J. P., Keppens, R., & Poedts, S. 2010, Advanced Magnetohydrodynamics (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Guo, L.J., Huang, Y.M., Bhattacharjee, A., & Innes, D. E. 2014, ApJ, 796, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Hester, J. J., Stone, J. M., Scowen, P. A., et al. 1996, ApJ, 456, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, A. S. 2016, MNRAS, 462, 2256 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, A., Berger, T., Isobe, H., & Shibata, K. 2012, ApJ, 746, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Innes, D. E., Cameron, R. H., Fletcher, L., Inhester, B., & Solanki, S. K. 2012, A&A, 540, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Isobe, H., Miyagoshi, T., Shibata, K., & Yokoyama, T. 2005, Nature, 434, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Jun, B.I., Norman, M. L., & Stone, J. M. 1995, ApJ, 453, 332 [NASA ADS] [CrossRef] [Google Scholar]
 Kane, J., Arnett, D., Remington, B. A., et al. 2000, ApJ, 528 [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, O., Komissarov, S. S., & Keppens, R. 2014, MNRAS, 443, 547 [NASA ADS] [CrossRef] [Google Scholar]
 Rayleigh, L. 1882, Proc. London Math. Soc., s1, 170 [CrossRef] [Google Scholar]
 Ristorcelli, J. R., & Clark, T. T. 2004, J. Fluid Mech., 507 [Google Scholar]
 Sharp, D. H. 1984, Physica D: Nonlinear Phenomena, 12, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Gardiner, T. 2007a, ApJ, 671, 1726 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Gardiner, T. A. 2007b, Phys. Fluids, 19, 094104 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJSS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, G. 1950, Proc. Roy. Soc. Lond. A, 201, 192 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Relevant parameter space explored of all simulations conducted, alongside measured nonlinear growth rate α and standard deviation σ.
All Figures
Fig. 1 The mixing region of RayleighTaylor unstable fluids of density 1 and 10, where the pure fluids above and below the mixing region are not shown. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 3 Fourier transform of B_{z} at t = 1 (a, d), t = 2 (b, e), and t = 6 (c, f), of the central zslice 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. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 5 Development of bubble height as a function of Agt^{2} 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. 

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

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

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

Open with DEXTER  
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/Agt^{2} 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. 

Open with DEXTER  
In the text 