How can largescale twisted magnetic structures naturally emerge from buoyancy instabilities?
Abstract
We consider the threedimensional instability of a layer of horizontal magnetic field in a polytropic atmosphere where, contrary to previous studies, the field lines in the initial state are not unidirectional. We show that if the twist is initially concentrated inside the unstable layer, the modifications of the instability reported by several authors (see e.g. cattaneo1990) are only observed when the calculation is restricted to two dimensions. In three dimensions, the usual interchange instability occurs, in the direction fixed by the field lines at the interface between the layer and the fieldfree region. We therefore introduce a new configuration: the instability now develops in a weakly magnetised atmosphere where the direction of the field can vary with respect to the direction of the strong unstable field below, the twist being now concentrated at the upper interface. Both linear stability analysis and nonlinear direct numerical simulations are used to study this configuration. We show that from the smallscale interchange instability, largescale twisted coherent magnetic structures are spontaneously formed, with possible implications to the formation of active regions from a deepseated solar magnetic field.
keywords:
Sun: magnetic fields, MHD1 Introduction
Active regions at the surface of the Sun are believed to be the visible manifestation of deepseated intense magnetic fields. The generation of these predominantly toroidal fields is associated with the differential rotation in the tachocline. However, the transport of these strong fields from the lower part of the convective zone up to the photosphere remains one of the major unknowns of the solar cycle. It is clear that superequipartition fields are buoyant and could therefore rise up to the surface, but the observed size, coherence and twist of the magnetic structures forming active regions remain largely unexplained.
In order to obtain some understanding of the transport of magnetic structures by buoyancy, highly idealised models have been considered. An early model (parker1955) considers the buoyancy properties of the magnetic field when it is concentrated in slender pressureconfined flux tubes. A simple model describing the dynamical evolution of such structures was proposed by spruit1981 and was subsequently used by several authors (see for example moreno1983). This model has successfully reproduced Joy’s law for the tilt of bipolar regions at the solar photosphere (choudhuri1987; dsilva1993; caligari1995) as well as the relationship between the tilt angle and the total flux of active regions (fan1994). Many authors have sought to move on from this simple model and have considered the full set of compressible magnetohydrodynamic (MHD) equations (or the anelastic approximation) to study the evolution of a simplified magnetic structure mimicking a flux tube. These models have shown that a certain amount of twist must be present inside the tube in order for it to remain coherent as it rises; such a twist corresponds to a nonzero axial current inside the flux tube. Without twist, the flux tube is shredded apart by the vorticity generated by its own motion (schussler1979; longcope1996; emonet1998; fan1998; jouve2007). Even more problematic than the issue of tube coherence is the interaction between the rising structures and the smallscale convective motions. Some studies have considered the rise of gradually twisted magnetic flux tubes through the convective zone both in Cartesian (fan2003; abbett2004) and spherical (jouve2009) geometries. Again, a given amount of twist is necessary for the flux tubes to remain coherent and rise through the convective zone. There is observational evidence of nonpotential magnetic fields in active regions (seehafer1990; leka1996; pevtsov2001), with the implication that the flux tubes are twisted before they emerge as active regions. While the twist is no more than one turn in general (chae2005), the corresponding Lorentz forces are implicated in chromospheric flux eruptions (torok2005). The origin of the twist is still unclear as there are various different explanations. The effect of the Coriolis force (fan2000) and differential rotation (devore2000) on the rising structures can only contribute a small fraction of the observed twist values (holder2004).
Other mechanisms have been suggested to explain the origin of nonpotential magnetic fields. One possibility is that the tube becomes twisted due to the effect of smallscale helical motions in the convective zone acting on the rising magnetic field, the socalled effect (longcope1998). Another suggestion, which is closely related to the content of this paper, is that the twist arises by accretion of poloidal fields during the rise of a flux tube (choudhuri2003; chatterjee2006). However indirect observational evidence suggests the presence of an initial twist in flux tubes at the base of the solar convective zone (holder2004), so that the twist may be produced during the formation process of the flux tube itself.
Another issue, which we don’t yet fully understand, is the mere existence of localised magnetic structures similar to flux tubes below the convective zone. What is known is that the differential rotation profile produces a strong toroidal field at the base of the solar convective zone. We can attempt to model the appearance of isolated flux structures by looking at a uniform layer of toroidal field where the vertical structure of the magnetic layer is initially imposed; then the horizontal length scale will be set naturally by the buoyancy instability. As a result of the magnetic pressure, the density of polytropic atmosphere in the magnetised region is reduced, assuming the temperature is not modified by the presence of the magnetic field. Consequently, the upper interface of the magnetic slab is susceptible to RayleighTaylortype instabilities. This paper is principally concerned with the threedimensional evolution of these instabilities when the initial magnetic field has nonzero twist. The twodimensional evolution of states with no initial twist was investigated by cattaneo1988; since no variation was permitted along the original field direction this treatment could only encompass the interchange instability, with the field lines remaining straight. This configuration has also been considered numerically in three dimensions by several authors (matthews1995; wissink2000), whereas alternative initial conditions (in which the temperature profile is altered by the magnetic field instead of the density) were considered by fan2001. It has proved very difficult to produce flux structures of significant size from this type of instability. Whatever the scale of the original instability secondary vortex instabilities due to the down flows between the plumes disrupt the magnetic structures leaving a rather diffuse and dynamically inactive field . Note that an alternative model has been suggested by kersale2007, where they considered the nonlinear evolution of magnetic buoyancy instabilities resulting from a smoothly stratified horizontal magnetic field. cattaneo1990 conducted twodimensional simulations of an initial twisted field structure. They found that in their constrained geometry large structures could appear as result of the twist. We will show in this paper that this effect disappears if we allow the motion to be fully threedimensional.
The earlier simulations to examine the evolution of buoyant magnetic structures assumed that the region above the initial magnetic slab was field free. In this case the initial potential energy contained in the initial condition is eventually transferred to kinetic energy and is finally dissipated. Any possible effects of the magnetic field existing above the unstable slab are neglected. Is this a reasonable simplification?
As shown by many simulations of overshooting convection in the presence of magnetic field (see for example tobias2001), the strong toroidal field initially injected in the system or generated by a shear flow mimicking the tachocline (guerrero2011) is predominantly contained in the stable radiative zone. However, despite the effects of turbulent pumping, a nonnegligible part of the field is still present in the convection zone, both due to the redistribution of largescale magnetic fields and to local smallscale dynamo action. If the strong toroidal field is buoyantly unstable, it will rise and interact with the small scale magnetic perturbations present in the convective zone.
The purpose of this paper is to investigate the simplest model that allows us to consider the buoyancy instability of a toroidal magnetic field in a weakly magnetised atmosphere above. The layer of strong toroidal field models the field induced by shear in the tachocline whereas the magnetised atmosphere above represents the magnetic perturbations in the convective zone. For this first, illustrative study, the field in the convective zone is assumed to be unidirectional but with a different direction as in the layer of strong field below. This is of course highly idealised, but allows us to derive a model involving few free parameters. The instability of a sheared magnetic field has already been studied (cattaneo1990; cattaneo1990b; kusano98; nozawa2005). In this case, the twist is uniformly distributed inside the unstable slab and the atmosphere above is fieldfree. To our knowledge, this type of instability has only been considered in published papers using twodimensional numerical simulations. The configuration we consider later in the paper is different as the field lines are unidirectional except in a thin region where the twist and current are concentrated. A similar setup was considered by stone2007a; stone2007b where they look at the effect of a transverse magnetic field on the RayleighTaylor instability.
The paper will proceed as follows: in section 2, we describe the governing equations, the model and physical parameters and the numerical methods. We then extend the results by cattaneo1990 and kusano98 to the threedimensional case in section 3. Finally, the interchange instability in a weakly magnetised atmosphere is considered in section 4, using both linear stability analysis and nonlinear numerical simulations in three dimensions.
2 Model and method
2.1 Model and governing equations
We consider the evolution of a planeparallel layer of compressible fluid, bounded above and below by two impenetrable, stressfree boundaries, a distance apart. The upper boundary is held at fixed temperature, whereas a vertical temperature gradient is fixed at the lower boundary. The geometry of this layer is defined by a Cartesian grid, with and corresponding to the horizontal coordinates. The axis points vertically downwards, parallel to the constant gravitational acceleration . The horizontal size of the fluid domain is defined by the aspect ratios and , so that the fluid occupies the domain , and . The physical properties of the fluid, namely the specific heats and , the shear viscosity , the thermal conductivity , the magnetic permeability and the magnetic diffusivity , are assumed to be constant. The model is identical to that used by, for example, matthews95b, silvers2009 and favier2012 .
It is convenient to introduce dimensionless variables, so we adopt the scalings described in matthews95b and bushby08. Lengths are scaled with the depth of the layer, . The temperature, , and the density, , are scaled with their values at the upper surface, and respectively. The velocity, , is scaled with the isothermal sound speed, , at the top of the layer, where is the gas constant. We adopt the same scaling for the Alfvén speed, which implies that the magnetic field, , is scaled with . Finally, we scale time by an acoustic time scale .
We now express the governing equations in terms of these dimensionless variables. The equation for conservation of mass is given by
(1) 
Similarly, the dimensionless momentum equation can be written in the following form,
(2) 
where is the pressure given by the equation of state and is the rate of strain tensor defined by
(3) 
Several nondimensional parameters appear in equation (2) including , the vertical temperature gradient at the lower boundary. Without magnetic field, it is initially equal to where is the temperature difference between the upper and lower boundaries. corresponds to the polytropic index. The dimensionless thermal diffusivity is given by and is the Prandtl number. The dimensionless strength of the magnetic field is defined by , where is the amplitude of the imposed magnetic field. It is related to the plasma by . The induction equation for the magnetic field is
(4) 
where is the ratio of magnetic to thermal diffusivity at the top of the layer (inverse Roberts number). The magnetic field is solenoidal so that
(5) 
Finally, the heat equation is
(6) 
where .
To complete the specification of the model, some boundary conditions must be imposed. In the horizontal directions, all variables are assumed to be periodic. As has already been described, the upper and lower boundaries are assumed to be impermeable and stressfree, which implies that at (the upper boundary) and (the lower boundary). Having nondimensionalised the system, the thermal boundary conditions at these surfaces correspond to fixing at and at . For the magnetic field boundary conditions, we choose appropriate conditions for perfectlyconducting boundaries, which implies that at and .
2.2 Initial conditions and parameters
In this paper, we investigate the instability of a layer of purely horizontal magnetic field. The vertical component of the magnetic field is therefore initially zero and we define a vertical profile for the initial horizontal field . In this paper, we consider two main initial magnetic field configurations. In section 3, the horizontal field is nonzero only in a prescribed layer , where is the depth of the unstable interface and is the lower boundary of the numerical domain. In section 4, we include a weak horizontal field above the strong horizontal field in the magnetic layer to give a magnetised atmosphere throughout the domain. This new configuration allows the magnetic tension to alter the development of the instability. For these prescribed initial magnetic fields an equilibrium static state () can be found by solving equations (2) and (6) for the density and temperature profiles respectively. The temperature gradient is slightly increased from the polytropic atmosphere due to the presence of ohmic heating, whereas the density profile is modified in order to keep the interface in pressure equilibrium (as seen later in figure 2(b)), making the layer of strong field buoyantly unstable. This magnetised polytropic layer, coupled with a small random thermal perturbation, is an appropriate initial condition for our numerical simulations.
Variable  Parameter  Value 

Thermal diffusivity  
Prandtl number  
Inverse Roberts number  
Thermal stratification  
Polytropic index  
Dimensionless magnetic field strength  
Horizontal aspect ratio 
There are many dimensionless parameters in this system and it is not viable to complete a systematic survey of the whole parameter space. Throughout this paper, the polytropic index is fixed at whereas the ratio of specific heats is given by (the appropriate value for a monatomic gas). This ensures that the initial polytropic atmosphere is stable to convective motions (even with the increase in the temperature gradient due to ohmic heating) and any resulting instability is due to magnetically induced buoyancy effects. It is also an appropriate value for the transition between the nearlyadiabatic solar convective zone and the stablystratified radiative zone below, which is the main focus of the present model.
We fix the kinetic Prandtl number to be and the inverse Roberts number to be for the work presented in this paper. Note that tachocline values are much lower than those that can be realised numerically. These two values are slightly larger than the ones observed in the literature. This is mostly justified by numerical reasons (larger values of and allow for larger time steps). However, we have run additional simulations using smaller values for and (down to ), along with cases with (as it should be for the tachocline), and the results are qualitatively unchanged. Note we did not explore the possibility of the double diffusive instability (silvers2009; silvers2010), which is much more demanding numerically and will be followed up at a later point. The thermal stratification is fixed to be , which corresponds to a moderately stratified layer (the layer contains approximately three pressure scale heights).
As will be explained later, increasing the stratification is not useful in the present model. We do not vary the dimensionless strength of the magnetic field in this paper and so is fixed to be . However we note here that we also tried to vary this parameter without any qualitative changes in the results. The value of in the solar tachocline is expected to be much smaller than the one considered here, but this regime is challenging numerically and is therefore left for future studies.






The numerical domain is chosen to be elongated in the direction in order to accommodate for the long wavelength modulation of the secondary instability in this direction, as reported by matthews1995. In the transverse direction, the numerical domain is large enough to accommodate approximately most unstable wavelengths in the case without transverse magnetic field. A summary of our choice of parameters is given in Table 1.
2.3 Numerical method
The given set of equations is solved using a modified version of the mixed pseudospectral/finite difference code originally described by matthews95b. Due to periodicity in the horizontal direction, horizontal derivatives are computed in Fourier space using fast Fourier transforms. In the vertical direction, a fourthorder finite differences scheme is used and an upwind stencil is applied for the advective terms. The timestepping is performed by an explicit adaptive thirdorder AdamsBashforth technique, with a variable timestep. The standard resolution is gridpoints in each horizontal direction and gridpoints in the vertical direction. We also consider quasi twodimensional simulations for which the variations along the direction (along the field lines of the unstable layer) are neglected. A poloidaltoroidal decomposition is used for the magnetic field in order to ensure that the field remains solenoidal.
A linear stability analysis of the equilibrium state is also performed in section 4.1 to determine how the initial buoyancy instability of the strong toroidal field is affected by the addition of the transverse field . The system of equations (1)(6) is perturbed, linearised and the ideal gas law and the magnetic solenoidality condition are used to eliminate the pressure and one component of the fluctuating magnetic field. All perturbations are of the form , where is the complex growth rate. This gives rise to the a seventh order system of equations for the fluctuating variables which may be numerically approximated by
(7) 
where is the vector of discretized eigenvectors and is the matrix of discretized differential operators, using a fourthorder finite difference approximation. The boundary conditions are applied through the relevant matrix coefficients using finite difference schemes. The resulting matrix equation is then solved for the eigenvalues and the corresponding eigenfunctions using the Linear Algebra PACKage (LAPACK).
3 Sheared magnetic field
In this section, we extend the configuration already considered by cattaneo1990 and kusano98 to the threedimensional case. We consider two different cases: in the first one, the initial field is unidirectional, say , whereas in the second case, a horizontal transverse field is added. For these two different cases, the horizontal magnetic field is the same and we compare both twodimensional (2D) and threedimensional (3D) simulations. By twodimensional simulation, we mean that all the variables are invariant in one of the horizontal direction (namely in the direction). The basic state is defined as follows: the initial horizontal field only depends on and is the same in both cases,
(8) 
where is the depth corresponding the top of the unstable layer whereas corresponds to the width of the transition between the layer and the nonmagnetised atmosphere above. This choice is compatible with the boundary conditions defined in section 2. In the following, we choose and . The first reference case is obtained by taking so that the field is initially unidirectional, in the direction. In the sheared case, and following cattaneo1990, the is chosen to be linearly dependent on the depth according to
(9) 
(this profile is actually smoothed using error function to avoid spurious numerical effects due to discontinuities). The depth at which the transverse field goes to zero is called the resonant layer by cattaneo1990. We choose here . The initial amplitude of the transverse field is derived from (the maximum value of initially occurs at and is ). Note that we tried different configurations with different resonant layers and different amplitudes for the field, but the results are qualitatively similar.


From this initial condition, equations (1)(6) are solved using the numerical method described in section 2.3. The simulations are run until the upper and lower boundary conditions have a significant effect on the generated flow. We note here that throughout all the results, the time is presented in terms of sound crossing time units. We present in figure 1 a comparison of twodimensional and threedimensional simulations for the two cases with and without transverse field. The cases with a unidirectional field are shown in the upper part of the figure, whereas the cases where initially are shown in the lower part. It can be seen that in the unidirectional case, both 2D and 3D simulations are very similar. The main difference is due to the lack of vortex stretching in the 2D simulation, leading to different mixing properties. The striking similarity between 2D and 3D simulations is a confirmation of the interchange nature of the initial instability, which is purely twodimensional. On the other hand, in the sheared case (i.e. when initially, lower line of figure 1), a drastic difference is observed between 2D and 3D simulations. In 2D, we recover the largescale instability already observed by several authors (see for example cattaneo1990 and kusano98). The smallscale interchange instability observed in the unidirectional case is altered by the tension due to the transverse field, leading to a long wavelength instability. This strong modification is not observed in 3D, where the usual interchange instability occurs. The similarity between the 3D simulations with and without field is clear. Figure 1 also shows a threedimensional visualisation comparing the unidirectional and sheared cases using the software VAPOR^{1}^{1}1http://www.vapor.ucar.edu/ (vapor1). It is clear that the nature of the instability is similar in both cases, apart from the change in direction of the field lines. In 2D, due to geometrical constraints, it is not possible for the field to undergo interchange instability, whereas it is possible in 3D. The results are similar for many different configurations (e.g. changing , , and the shape of the field from linear to quadratic). Thus we conclude from what we have observed that imposing a transverse field inside the unstable layer is not enough to alter interchange instability in three dimensions, as was previously thought using 2D simulations. We cannot rule out the possibility that this conclusion depends on the parameters and initial conditions considered, but we simply conclude that it might not be as easy as in the 2D case.
In the next section, we examine another initial configuration, which leads to a significant modification of the instability, while working both in two and three dimensions.
4 Instability in a magnetised atmosphere
In this section, unlike in section 3 and previous studies, the atmosphere above the layer of strong field is now also magnetised. Initially, we fix the vertical field to be zero everywhere whereas the horizontal field is only varying in the vertical direction according to
(10) 
where is the depth corresponding the top of the unstable layer and corresponds the width of the transition between the layer and the weakly magnetised atmosphere. As in section 3, we choose and . In the following, the region is called the atmosphere whereas the region is called the layer. The component of the magnetic field is derived from
(11) 
The profile of horizontal field can be seen on figure 2(a), for .
We define the pitch angle as being so that varies between and , depending on the value of . In the present setup, the horizontal field is six times weaker in the atmosphere than in the layer. The actual amplitude of the field in the atmosphere is rather unimportant and we conducted additional simulations varying this parameter. The results were qualitatively similar apart from the cases where the vertical gradients of the component of the field become strong enough to undergo interchange instability, which makes the analysis of the results difficult. We therefore choose to focus on the simpler case where is large enough to have a dynamical effect on the development of the instability but small enough in order to remain stable to interchange modes during the simulations. We would like to stress that this setup is different from the case presented in section 3, where the twist was concentrated inside the unstable layer. Here, we consider the usual interchange instability altered by the presence of an oblique field in the atmosphere above.
4.1 Linear analysis
From the linear stability analysis we find that the most readily destabilised modes are interchange modes (). The instability is localised at the interface between the magnetic layer and the magnetised atmosphere. Figure 3 shows how the growth rate of the instability varies as a function of the horizontal wave number for different pitch angles . As the pitch angle is increased the growth rate of the most unstable mode decreases due to the increase in magnetic tension at the interface. This lower growth rate means that the instability will take longer to develop for larger values of the twist parameter . The mode of maximum growth also decreases as the twist is increased, however the wavelength of the most unstable mode remains smaller than the depth of the domain for all values of the twist parameter , this will be discussed further in the following sections. Although our configuration differs, our conclusion is qualitatively similar to the one drawn by cattaneo1990 and kusano98, i.e. as the intensity of the transverse field is increased, both the growth rate of the instability and the transverse wave number decrease.

4.2 Nonlinear 3D numerical simulations
In this section, we discuss the result from nonlinear simulations in 3D. The parameter that we vary here is the initial pitch angle of the field in the atmosphere contained between and . We first run a reference simulation where is initially zero everywhere. This corresponds to the classical interchangetype instability already considered in many papers (see, for example, matthews1995 and wissink2000), apart from the initial presence of a weak magnetic field in the atmosphere above the layer. Using this as a reference case, we examine three different simulations that only differ by the parameter as defined in equation (11). The parameter varies from (unidirectional field, ) to (). The reason why we don’t consider angles larger than will be explained later. For each of these initial configurations, we solve equations (1)(6) using the numerical method described previously.
We present in figure 4 vertical slices in the plane where contours of the toroidal field are plotted at different times. Bright colours correspond to large values whereas dark tones correspond to low values of . The upper part of the figure correspond to the unidirectional field case (i.e. initially). As already observed in section 3, the initial instability consists of alternation of upward and downward flows. As they increase in amplitude, a secondary instability, of KelvinHelmholtz type, develops in the shear region generated at the interface between the upward and downward motions. This tends to generate horizontal vorticity, which further mixes field lines with different amplitudes of . The end product of this type of simulations is qualitatively the same for different parameters: the initial coherent layer is disrupted by the vortical motion generated by the secondary instability. The resulting toroidal field is diffuse and no intense coherent magnetic structures are observed. It is worth noting that the presence of a weak field in the upper part of the domain doesn’t modify the nature of the instability, as can be seen comparing the top middle panel of figure 4 with the top middle panel of figure 1. Note however that, as noted by vasil2009, one way to damp the secondary instability is to consider very large values of the magnetic Prandtl number. By doing so, the secondary instability is damped by strong viscosity whereas the magnetic field remains coherent. It is however difficult to justify such a parameter regime as the magnetic Prandtl number based on molecular diffusivities is expected to be very small in the tachocline (gough2007).






Now consider a second simulation where the component of the magnetic field is initially non zero in the atmosphere above the unstable layer. We first consider the case where the pitch angle is (as depicted in figure 2(a)) in order to illustrate the dynamical effect of the transverse field . Initially, as seen on the lower left panel of figure 4, the nature of the instability is unchanged, apart from a larger characteristic horizontal length scale, as predicted by the previous linear analysis of section 4.1. From the simulation, the most unstable wave number in the direction is roughly twice as small in the case than in the reference case (see also the magnetic energy spectra in figure 9 below). At later times (, lower middle panel of figure 4), it is apparent that the nonlinear development of the instability is also modified due the magnetic tension existing at the interface between the unstable layer and the atmosphere above. Figure 4 shows an alternation of strong and weak down flows due to the merging of small magnetic structures into larger ones. The reduction of the growth rate as predicted by the linear analysis of section 4.1 is also apparent, as the magnetic fluctuations are much more spread across the layer in the unidirectional case. At the final time , the difference between the two cases is very clear: the presence of the field in the atmosphere tends to generate large coherent magnetic structures in the direction containing fluctuations on smaller scales generated during the initial stages of the instability. The toroidal field is also much less diffuse than in the unidirectional case, and a clear interface between field and fieldfree regions is observed. At later times, the magnetic and velocity fluctuations are too close to the upper and lower boundaries, and we stop the simulations to avoid spurious confinement effects.




It is informative to look in details at one of the merging event, where two small magnetic structures create a larger one in order to overcome the magnetic tension existing at the interface. A time series of a subpart of the plane is presented on figure 5. We plot both the current density as contours along with arrows representing the velocity field in this particular vertical plane. Initially, a strong down flow in the middle of the panel is detected, dragging down the magnetic field . This vertical down flow must also advect the weak field existing in the atmosphere, therefore working against magnetic tension. At time (third panel), it is apparent that the downward flow is now weaker and disconnected from the upper atmosphere. The secondary instability still occurs as a vortex has been created inside the layer. However, the vorticity is not continuously fed by the down flow. The final stage consists of a unique magnetic structure containing small scale perturbations slowly diffusing away. This reconnection process is clearly visible in the current density, where two rising structures are pushed together due to magnetic tension, leading to the formation of one large current sheet at the interface.
Finally, and before comparing quantitatively the different simulations, we present in figure 6 a volume rendering of the magnetic energy . Again, we compare the unidirectional case with the case , at the same time . Without any initial transverse field , the results are similar to what was already observed by matthews1995 and wissink2000. The initial linear instability is purely twodimensional and corresponds to the interchange of horizontal straight field lines. Eventually the interaction of counterrotating vortices generated by the secondary instability leads to the arching of the magnetic structures. Note that, although these magnetic structures are clearly observed on figure 6, they are very weak and lack twist. The use of isosurfaces, as in wissink2000, can be misleading as it gives the impression the magnetic structures are coherent with a sharp interface, which is not the case. We therefore expect that the presence of overshooting convective motions at the interface between the convective and radiative zones will inevitably disrupt these smallscale magnetic features. In the case however, we observe larger magnetic structures, containing a significant amount of twist. Note that these larger structures are far from being welldefined localised flux tubes (cattaneo2006). We however stress that the fact that our polytropic atmosphere is stablystratified and the magnetic tension resulting from the presence of the oblique field tends to act against the rise of these structures. These limiting effects are the results of our oversimplified model, but more refined approaches could eventually lead to the formation of localised flux tubes from the breakup of a magnetic layer. The last panel on the right in figure 6 shows a closeup with several magnetic field lines in blue. As the strong field rises, twist accumulates at the apex of the magnetic structure, damping its disruption by vortical motions and increasing its longevity against eventual convective motions. This is reminiscent of the scenario suggested by choudhuri2003, where the twist is generated by the rise of strong toroidal fields inside the convective zone containing a poloidal field. Note however that choudhuri2003 considered already formed flux tubes and not a uniform layer of magnetic field.






We now discuss more quantitatively the effect of the transverse field on the development of the instability. Figure 7 shows the evolution with time of several global quantities. In the following, the brackets denote a spatial average over all coordinates. In particular, we define the total kinetic energy and the total magnetic energy . As previously, we compare the results from simulations without initial field () and with various pitch angles. After the decay of the initial perturbation, the kinetic energy grows exponentially for a few sound crossing times, as seen on the left panel of figure 7. Note that, as predicted by the linear analysis of section 4.1, the growth rate decreases when increases.
After approximately ten sound crossing times, the instability saturates. The evolution after this stage strongly depends on . For the case , the kinetic energy decreases with time. In this case, most of the initial potential energy contained in the initial condition is released during the linear phase, leading to a decaying nonlinear regime. For , much less kinetic energy is initially released due to the magnetic tension at the interface, leading to a growth of the kinetic energy with time up to the end of the simulation. Of course, as no external source of energy is included in the model, we expect the kinetic energy to ultimately decay in all cases. One of the important steps of the magnetic buoyancy instability is the generation of strong horizontal vortices by the secondary KelvinHelmholtz instability. It is clear that the field, if strong enough, will tend to damp this generation of vorticity. In order to create strong down and up flows, it is now necessary to work against the tension of the weak magnetic field . The evolution of the total enstrophy in the numerical domain is plotted versus time on figure 7, where . In all cases, we observe a rapid increase of the enstrophy with time. As the instability saturates, the growth rate of the enstrophy decreases. As for the kinetic energy, we expect the enstrophy to ultimately decay. The main effect of the transverse field is to reduce the production of enstrophy by the secondary instability, so that the enstrophy is always smaller in the cases where than in the case where . This leads to a reduction of the viscous dissipation, therefore reducing the rate at which the kinetic energy is dissipated.
The two others panels in figure 7 show the evolution with time of the total magnetic energy and current density , where . In all cases, the magnetic energy decays whereas the current density increases with time. As increases, the magnetic energy decreases less rapidly. This is due to the damping of the vortical motions that reduces the mixing of the magnetic field lines and the creation of dissipative currents. Note that the current is initially slightly larger when , as expected from the initial conditions.
Following stone2007a, it is useful to quantify the amount of mixing due to the instability. Instead of defining the mixing in terms of density, as for the RayleighTaylor instability, we define the fraction of strong horizontal field as
(12) 
where . corresponds to the initial amplitude of the horizontal field at the upper boundary, and can be considered as the weak field. is the initial horizontal field at the lower boundary and is therefore the strong field. The fraction of weak field is simply . The amount of mixing in the system can be estimated by the following quantity
(13) 
which is equal to when the averaged horizontal field is equal to its initial value at the vertical boundaries, and if the field is perfectly mixed and is equal to the average between these two values.

Figure 8 shows the mixing parameter as a function of for the four cases considered at . The thin black line corresponds to the value of derived from the initial condition, which is nearly the same in all four cases. As time increases, the two fronts of mixing propagate vertically. It is clear that due to the magnetic tension at the interface, the overall mixing of the upper front is reduced. The vertical extent of the mixing region is reduced as increases. This result can be explained by the reduced growth rate on the initial linear instability (see figure 3) in the presence of a transverse field and by the damping of the secondary instability leading to weaker vortical motions (see figure 7). In the unidirectional case, the potential energy is released to counter the stable stratification, leading to the rise of the toroidal magnetic field. When , the magnetic tension also acts against the instability. This is a clear limitation of our simple model, and will be further discussed in section 5. Note finally that figure 8 also shows that the lower front is roughly unchanged by the presence of transverse field.


The generation of largescale coherent magnetic structures is clearly visible in the magnetic energy spectra defined by
(14) 
where is the onedimensional Fourier transform in the direction of , whereas the star denotes the complex conjugation. We focus here on the magnetic energy spectra in the direction perpendicular to the initial strong field (the results in the direction are rather similar between the different cases).
Figure 9 shows the magnetic energy spectrum for and , at four different times spanning from the initial linear instability () to the late stage of the nonlinear phase (). Initially, the most unstable horizontal wave number corresponds to the linear results already discussed in section 4.1. The most unstable wave number is roughly for and for , which can be compared to the results presented in figure 3. The time evolution of the magnetic energy spectra can be described in two phases. First, the smallscale instability grows exponentially from to . The effect of the transverse field is to reduce the initial linear growth rate, as already discussed in figures 3 and 7. From , the evolution of the instability is constrained by the mixing due to the vortical motions. Here, the effect of the transverse field is to reduce the direct cascade of magnetic energy toward small scales, as we observe steeper energy spectra for than for , especially at early times. Finally, during the late stage of the instability, we observe a decay of the magnetic energy at small scales whereas the magnetic energy at large scales remains constant for . In contrast, the magnetic energy is still growing at large scales for . The effect of the transverse field is indeed clearly visible on the small wave numbers region . For , a clear peak is observed at large scales, which corresponds to the coherent magnetic structures already observed on figures 4 and 6. This is also consistent with the larger values of the global magnetic energy observed on figure 7 when . Without transverse field, the magnetic energy peaks at larger horizontal wave numbers. Note that this difference cannot be attributed to linear effects, as the most unstable wave numbers as predicted by the linear analysis are much larger (see figure 3). During the nonlinear phase, the smallscale features (i.e. ) are roughly independent of , and we observe a continuity of magnetic scales up to the resistive scale. There is no clear scaling in our case since the diffusivities do not allow for a fullydeveloped turbulent state.

Finally, the amount of twist in the magnetic field lines, as observed in figure 6, can be estimated looking at the current helicity, defined by . Since both magnetic energy and current density are timedependent (see figure 7), we consider the relative current helicity defined by
(15) 
Note that given our initial configuration, the amount of current helicity depends on . We show on figure 10 the evolution with time of for different values of . The effect of the initial condition is clearly apparent, as when , whereas initially increases with . Without initial transverse field, the current helicity remains statistically negligible, although we observe a very weak positive correlation, which might be due to the interaction between the vorticity and the rising magnetic field. However, as soon as , we observe an increase of the relative current helicity from its small initial value. This is particularly visible in the case. After a transient stage where increases and then decreases, we observe a monotonic increase of the relative current helicity up to the end of the simulation. Note that continues to increase even after the saturation of the initial linear phase, which roughly corresponds to (see figure 7). Again, since the atmosphere is stably stratified, we expect the magnetic field to ultimately stop rising and the current helicity to ultimately decay. If the entropy gradient was nearly negligible, as expected in the convective zone, the magnetic structures observed in our simple model could continue to rise unimpeded, eventually becoming kinkunstable (fan1999). The increase in the relative current helicity observed here is an indication that the twist observed in active regions could be due to the interaction between a rising toroidal structure and the weaker field inside the convective zone, as suggested by choudhuri2003. Of course, due to numerical constraints, the large diffusivities considered in this paper do not allow for the formation of isolated structures.
5 Discussion
In this work, we sought to examine the instability of a layer of horizontal magnetic field in a polytropic atmosphere, where the direction of the field lines is depthdependent. In section 3 we examined the instability of a layer of horizontal magnetic field in a polytropic atmosphere, where the direction of the field lines in the layer varied with depth. We showed that the initial idea of building largescale coherent magnetic structures from the buoyancy instability of a sheared magnetic layer, as studied by cattaneo1990 and kusano98, seemed limited to twodimensional geometry. The same configuration drastically changed in three dimensions, where the interchange instability was able to occur in an oblique direction. We stress that it could be possible to find a configuration of this type that strongly modifies the nature of the magnetic buoyancy instability, however this was not the aim of the paper.
In section 4 we introduced another initial configuration, where the twist was concentrated at the interface between a strong layer of horizontal magnetic field and a weakly magnetised atmosphere above. The weakly magnetised atmosphere was justified by the presence of convection, which is expected to contain nonnegligible magnetic fluctuations, both due to the redistribution of largescale magnetic fields, and due to local smallscale dynamo action. The presence of a strong unidirectional field was justified by the presence of radial differential rotation in the tachocline. We presented a linear stability analysis and numerical simulations of the magnetic buoyancy instability in sections 4.1 and 4.2. We have shown that the initial interchange instability is only weakly modified by the presence of the transverse magnetic field. During the nonlinear phase of the instability, the smallscale magnetic structures merged together to overcome magnetic tension, leading to the formation of magnetic structures on scales larger than those of the initial linear instability. These structures were shown to be more coherent than those in the unidirectional case, since the production of vorticity by the secondary instability was reduced. In addition we found a significant amount of twist was generated due to the rise of the toroidal field in the weakly magnetised atmosphere. This mechanism could provide the initial twist necessary for the magnetic structures to rise coherently through the convective zone.
It is however important to stress the limitations of the present model. In order to limit the number of parameters, we have considered two unidirectional fields with a current sheet at the interface. One could argue that the magnetic field in the convective zone is interacting with smallscale motions and is certainly not unidirectional as assumed here. Although it remains to be verified using a more refined model, we expect that locally the toroidal field will behave similarly in the presence of smallscale magnetic fluctuations.
There are two stabilising effects in the present model, the stably stratified atmosphere (which is necessary to suppress convective motions) and the magnetic tension generated by the twisted atmosphere. The tension becomes important as the atmosphere is compressed by the rising structures, and the toroidal magnetic field will eventually stop rising. We have reached this regime in some of our simulations where the pitch angle is equal or larger that . The stabilising effect of the tension is clearly overestimated by the present model, as the tension exists everywhere at the interface between the rising layer and the atmosphere. In a more realistic situation, the convective zone would be filled with smallscales magnetic structures. Locally these structures could have a similar effect to the global tension considered in our simulations, although we didn’t check it numerically here. We plan in a forthcoming paper to look at a configuration where the global magnetic tension is removed at some depth, leading the rise and formation of local coherent twisted magnetic structures. Another possibility would be to extend the work of fan2003 and jouve2009, assuming that the convective zone is filled with magnetic perturbations. This would confirm whether the generation of twist, via the interaction of a rising flux tube and smallscale magnetic perturbations, is sufficient for the flux tube to remain coherent as it rises through the convective zone. Such simulations are underway in spherical geometry using the ASH code (pinto2012).
We conclude by arguing that the merging mechanism illustrated by this simple model could have implications on the formation of largescale magnetic structures from a deepseated toroidal field. The twist observed in active regions could therefore be attributed to a progressive accumulation of local twist as the toroidal field rises through the magnetised convective zone.
Acknowledgements This work has been financially supported by STFC. CPU time was provided by the HPC resources of CALMIP under the allocation 2012P1118 and by the UKMHD supercomputing facility located in Warwick.