# How can large-scale twisted magnetic structures naturally emerge from buoyancy instabilities?

Benjamin Favier, Laurène Jouve, William Edmunds, Lara J. Silvers, Michael

R. E. Proctor

► To cite this version:

Benjamin Favier, Laurène Jouve, William Edmunds, Lara J. Silvers, Michael R. E. Proctor. How can large-scale twisted magnetic structures naturally emerge from buoyancy instabilities?. Monthly Notices of the Royal Astronomical Society, 2012, 10.1111/j.1365-2966.2012.21920.x. [https://hal.science/hal-01287839](https://hal.science/hal-01287839)

HAL Id: hal-01287839

https://hal.science/hal-01287839

Submitted on 12 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# How can large-scale twisted magnetic structures naturally emerge from buoyancy instabilities?

B. Favier,<fnref n="1" />* L. Jouve,<fnref n="2" />,<fnref n="3" /> W. Edmunds,<fnref n="1" /> L. J. Silvers<fnref n="4" /> and M. R. E. Proctor<fnref n="1" />

Accepted 2012 August 13. Received 2012 June 7

## ABSTRACT

We consider the 3D 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 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 field-free region. We therefore introduce a new configuration: the instability now develops in a weakly magnetized 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 non-linear direct numerical simulations are used to study this configuration. We show that from the small-scale interchange instability, large-scale twisted coherent magnetic structures are spontaneously formed, with possible implications to the formation of active regions from a deep-seated solar magnetic field.

**Key words:** instabilities – magnetic fields – MHD – Sun: magnetic topology.

## 1 INTRODUCTION

Active regions at the surface of the Sun are believed to be the visible manifestation of deep-seated 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 super-equipartition 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 idealized models have been considered. An early model (Parker 1955) considers the buoyancy properties of the magnetic field when it is concentrated in slender pressure-confined flux tubes. A simple model describing the dynamical evolution of such structures was proposed by Spruit (1981) and was subsequently used by several authors (see for example Moreno-Insertis 1983). This model has successfully reproduced Joy's law for the tilt of bipolar regions at the solar photosphere

(Choudhuri 1987; D'Silva & Choudhuri 1993; Caligari, Moreno-Insertis & Schüssler 1995) as well as the relationship between the tilt angle and the total flux of active regions (Fan, Fisher & McClymont 1994). Many authors have sought to move on from this simple model and have considered the full set of compressible magnetohydrodynamic 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 non-zero axial current inside the flux tube. Without twist, the flux tube is shredded apart by the vorticity generated by its own motion (Schüssler 1979; Emonet & Moreno-Insertis 1996; Longcope, Fisher & Arendt 1996; Fan, Zweibel & Lantz 1998; Jouve & Brun 2007). Even more problematic than the issue of tube coherence is the interaction between the rising structures and the small-scale convective motions. Some studies have considered the rise of gradually twisted magnetic flux tubes through the convective zone both in Cartesian (Fan, Abbett & Fisher 2003; Abbett et al. 2004) and spherical (Jouve & Brun 2009) 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 non-potential magnetic fields in active regions (Seehafer 1990; Leka et al. 1996; Pevtsov, Canfield & Latushko 2001), with the implication that the flux tubes are twisted

<footnote n="1">¹Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA</footnote>

<footnote n="2">²Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, UPS-OMP, Toulouse, France</footnote>

<footnote n="3">³Institut de Recherche en Astrophysique et Planétologie, CNRS, 14 Avenue Édouard Belin, 31400 Toulouse, France</footnote>

<footnote n="4">⁴Centre for Mathematical Science, City University London, Northampton Square, London EC1V 0HB</footnote>

<footnote>*E-mail: b.favier@damtp.cam.ac.uk, bff21@damtp.cam.ac.uk</footnote>
before they emerge as active regions. While the twist is no more
than one turn in general (Chae & Moon 2005), the correspond-
ing Lorentz forces are implicated in chromospheric flux eruptions
(Török & Kliem 2005). The origin of the twist is still unclear as
there are various different explanations. The effect of the Coriolis
force (Fan & Gong 2000) and differential rotation (DeVore 2000)
on the rising structures can only contribute a small fraction of the
observed twist values (Holder et al. 2004).

Other mechanisms have been suggested to explain the origin
of non-potential magnetic fields. One possibility is that the tube
becomes twisted due to the effect of small-scale helical motions in
the convective zone acting on the rising magnetic field, the so-called
Σ effect (Longcope, Fisher & Pevtsov 1998). 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 (Choudhuri 2003; Chatterjee, Choudhuri & Petrovay 2006).
However, indirect observational evidence suggests the presence of
an initial twist in flux tubes at the base of the solar convective zone
(Holder et al. 2004), so that the twist may be produced during the
formation process of the flux tube itself.

Another issue, which we do not yet fully understand, is the mere
existence of localized magnetic structures similar to flux tubes be-
low 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 den-
sity of a polytropic atmosphere in the magnetized region is reduced,
assuming that the temperature is not modified by the presence of the
magnetic field. Consequently, the upper interface of the magnetic
slab is susceptible to Rayleigh–Taylor-type instabilities. This paper
is principally concerned with the 3D evolution of these instabilities
when the initial magnetic field has non-zero twist. The 2D evolu-
tion of states with no initial twist was investigated by Cattaneo &
Hughes (1988); 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 configura-
tion has also been considered numerically in three dimensions by
several authors (Matthews, Hughes & Proctor 1995a; Wissink et al.
2000), whereas alternative initial conditions (in which the temper-
ature profile is altered by the magnetic field instead of the density)
were considered by Fan (2001). It has been proved very difficult to
produce flux structures of significant size from this type of instabil-
ity. Whatever the scale of the original instability, secondary vortex
instabilities due to the downflows between the plumes disrupt the
magnetic structures leaving a rather diffuse and dynamically inac-
tive field. Note that an alternative model has been suggested by
Kersalé, Hughes & Tobias (2007), where they considered the non-
linear evolution of magnetic buoyancy instabilities resulting from
a smoothly stratified horizontal magnetic field. Cattaneo, Chiueh
& Hughes (1990a) conducted 2D simulations of an initial twisted
field structure. They found that in their constrained geometry, large
structures could appear as a result of the twist. We will show in
this paper that this effect disappears if we allow the motion to be
fully 3D.

The earlier simulations to examine the evolution of buoyant mag-
netic structures assumed that the region above the initial magnetic
slab was field free. In this case the potential energy contained in
the initial condition is eventually transferred to kinetic energy and
is finally dissipated. Any possible effects of the magnetic field ex-

isting above the unstable slab are neglected. Is this a reasonable
simplification?

As shown by many simulations of overshooting convection in
the presence of a magnetic field (see for example Tobias et al.
2001), the strong toroidal field initially injected in the system or
generated by a shear flow mimicking the tachocline (Guerrero &
Käpylä 2011) is predominantly contained in the stable radiative
zone. However, despite the effects of turbulent pumping, a non-
negligible part of the field is still present in the convection zone,
both due to the redistribution of large-scale magnetic fields and
due to local small-scale 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 magnetized atmosphere above. The layer
of strong toroidal field models the field induced by shear in the
tachocline, whereas the magnetized 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 idealized but allows us
to derive a model involving few free parameters. The instability of
a sheared magnetic field has already been studied (Cattaneo et al.
1990a; Cattaneo, Chiueh & Hughes 1990b; Kusano, Moriyama &
Miyoshi 1998; Nozawa 2005). In this case, the twist is uniformly
distributed inside the unstable slab and the atmosphere above is
field free. To our knowledge, this type of instability has only been
considered in published papers using 2D 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
Stone & Gardiner (2007a,b), where they looked at the effect of a
transverse magnetic field on the Rayleigh–Taylor 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 Cattaneo et al. (1990a) and Kusano et al. (1998) to the 3D case in Section 3. Finally, the interchange instability in a weakly magnetized atmosphere is considered in Section 4, using both linear stability analysis and non-linear numerical simulations in three dimensions.

2 MODEL AND METHOD

## 2.1 Model and governing equations

We consider the evolution of a plane-parallel layer of compressible
fluid, bounded above and below by two impenetrable, stress-free
boundaries, a distance *d* apart. The upper boundary is held at fixed
temperature *T*₀, whereas a vertical temperature gradient θ is fixed
at the lower boundary. The geometry of this layer is defined by a
Cartesian grid, with *x* and *y* corresponding to the horizontal co-
ordinates. The *z*-axis points vertically downwards, parallel to the
constant gravitational acceleration *g* = *g*̂. The horizontal size of
the fluid domain is defined by the aspect ratios λₓ and λᵧ, so that
the fluid occupies the domain 0 < *z* < *d*, 0 < *x* < λₓ*d* and 0 <
*y* < λᵧ*d*. The physical properties of the fluid, namely the specific
heats *c*<sub>p</sub> and *c*<sub>v</sub>, the shear viscosity *μ*, the thermal conductivity *K*,
the magnetic permeability *μ*₀ and the magnetic diffusivity η, are
assumed to be constant. The model is identical to that used by, for
example, Matthews, Proctor & Weiss (1995b), Silvers et al. (2009)
and Favier & Bushby (2012).
It is convenient to introduce dimensionless variables, so we adopt the scalings described in Matthews et al. (1995b) and Bushby et al. (2008). Lengths are scaled with the depth of the layer, *d*. The temperature, *T*, and the density, *ρ*, are scaled with their values at the upper surface, $T_0$ and $\rho_0$, respectively. The velocity, $\mathbf{u}$, is scaled with the isothermal sound speed, $\sqrt{R_*T_0}$, at the top of the layer, where $R_*$ is the gas constant. We adopt the same scaling for the Alfvén speed, which implies that the magnetic field, $\mathbf{B}$, is scaled with $\sqrt{\mu_0\rho_0 R_* T_0}$. Finally, we scale time by an acoustic time-scale $d/\sqrt{R_*T_0}$.

We now express the governing equations in terms of these dimensionless variables. The equation for conservation of mass is given by

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0. \qquad (1) $$

Similarly, the dimensionless momentum equation can be written in the form

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\frac{1}{\rho} \nabla P + \frac{F}{\rho} (\nabla \times \mathbf{B}) \times \mathbf{B} + \theta(m+1)\hat{z} + \frac{\kappa \sigma}{\rho} \nabla \cdot \boldsymbol{\tau}, \quad (2) $$

where $P$ is the pressure given by the equation of state $P = \rho T$ and $\boldsymbol{\tau}$ is the rate of strain tensor defined by

$$ \tau_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} - \frac{2}{3} \delta_{ij} \frac{\partial u_k}{\partial x_k}. \qquad (3) $$

Several non-dimensional parameters appear in equation (2) including $\theta$, the vertical temperature gradient at the lower boundary. Without a magnetic field, it is initially equal to $\Delta T/T_0$, where $\Delta T$ is the temperature difference between the upper and lower boundaries. $m = gd/R_*\Delta T - 1$ corresponds to the polytropic index. The dimensionless thermal diffusivity is given by $\kappa = K/d\rho_0 c_p (R_* T_0)^{1/2}$, and $\sigma = \mu c_p/K$ is the Prandtl number. The dimensionless strength of the magnetic field is defined by $F = B_0^2/(\mu_0\rho_0 R_* T_0)$, where $B_0$ is the amplitude of the imposed magnetic field. It is related to the plasma $\beta$ by $F = 2/\beta$. The induction equation for the magnetic field is

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B} - \zeta_0 \kappa \nabla \times \mathbf{B}), \qquad (4) $$

where $\zeta_0 = \eta c_p \rho_0 / K$ is the ratio of magnetic to thermal diffusivity at the top of the layer (inverse Roberts number). The magnetic field is solenoidal so that

$$ \nabla \cdot \mathbf{B} = 0. \qquad (5) $$

Finally, the heat equation is

$$ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = -(\gamma - 1) T \nabla \cdot \mathbf{u} + \frac{\kappa \gamma}{\rho} \nabla^2 T + \frac{\kappa(\gamma - 1)}{\rho} \left( \sigma \tau^2 / 2 + \zeta_0 F |\nabla \times \mathbf{B}|^2 \right), \qquad (6) $$

where $\gamma = c_p/c_v$.

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 stress free, which implies that $u_{x,z} = u_{y,z} = u_z = 0$ at $z = 0$ (the upper boundary) and $z = 1$ (the lower boundary). Having non-dimensionalized the system, the thermal boundary conditions at these surfaces correspond to fixing $T = 1$ at $z = 0$ and $\partial_z T = \theta$ at $z = 1$. For the magnetic field boundary conditions, we choose appropriate conditions for perfectly conducting boundaries, which implies that $B_z = B_{x,z} = B_{y,z} = 0$ at $z = 0$ and $1$.

## 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 $B_z$ is therefore initially zero, and we define a vertical profile for the initial horizontal field $B_{H(z,t=0)} = \sqrt{B_z^2(z) + B_{\theta(z)}^2}$. In this paper, we consider two main initial magnetic field configurations. In Section 3, the horizontal field is non-zero only in a prescribed layer $z_t < z < 1$, where $z_t$ is the depth of the unstable interface and $z = 1$ 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 magnetized 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 ($\mathbf{u} = \mathbf{0}$) 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 Fig. 2b), making the layer of strong field buoyantly unstable. This magnetized polytropic layer, coupled with a small random thermal perturbation, is an appropriate initial condition for our numerical simulations.

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 $m = 1.6$, whereas the ratio of specific heats is given by $\gamma = 5/3$ (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 nearly adiabatic solar convective zone and the stably stratified radiative zone below, which is the main focus of the present model.

We fix the kinetic Prandtl number to be $\sigma = 0.05$ and the inverse Roberts number to be $\zeta_0 = 0.05$ for the work presented in this paper. Note that tachocline values are much lower than those that can be realized numerically. These two values are slightly larger than the ones observed in the literature. This is mostly justified by numerical reasons (larger values of $\sigma$ and $\zeta_0$ allow for larger time steps). However, we have run additional simulations using smaller values for $\sigma$ and $\zeta_0$ (down to $\sigma = \zeta_0 = 10^{-2}$), along with cases with $\sigma < \zeta_0$ (as it should be for the tachocline), and the results are qualitatively unchanged. Note that we did not explore the possibility of the double diffusive instability (Silvers et al. 2009, 2010), which is much more demanding numerically and will be followed up at a later point. The thermal stratification is fixed to be $\theta = 2$, 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 $F$ is fixed to be $F = 0.1$. However, we note here that we also tried to vary this parameter without any qualitative changes in the results. The value of $F$ 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 y direction in order to accommodate the long-wavelength modulation of the secondary instability in this direction, as reported by Matthews et al. (1995a). In the transverse direction, the numerical domain is large enough to accommodate approximately 20 most unstable
<table><caption>**Table 1.** Parameter values considered in this paper.</caption><thead><tr><th>Variable</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>κ</td><td>Thermal diffusivity</td><td>5 × 10<sup>-4</sup></td></tr><tr><td>σ</td><td>Prandtl number</td><td>5 × 10<sup>-2</sup></td></tr><tr><td>ζ<sub>0</sub></td><td>Inverse Roberts number</td><td>5 × 10<sup>-2</sup></td></tr><tr><td>θ</td><td>Thermal stratification</td><td>2</td></tr><tr><td>m</td><td>Polytropic index</td><td>1.6</td></tr><tr><td>F</td><td>Dimensionless magnetic field strength</td><td>0.1</td></tr><tr><td>λ<sub>x</sub> × λ<sub>y</sub></td><td>Horizontal aspect ratio</td><td>2 × 8</td></tr></tbody></table>

wavelengths in the case without a 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 pseudo-spectral/finite difference code originally described by Matthews et al. (1995b). Due to periodicity in the horizontal direction, horizontal derivatives are computed in Fourier space using fast Fourier transforms. In the vertical direction, a fourth-order finite difference scheme is used and an upwind stencil is applied for the advective terms. The time-stepping is performed by an explicit adaptive third-order Adams–Bashforth technique, with a variable time step. The standard resolution is 256 grid points in each horizontal direction and 240 grid points in the vertical direction. We also consider quasi-2D simulations for which the variations along the y direction (along the field lines of the unstable layer) are neglected. A poloidal–toroidal 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 $B_y$ is affected by the addition of the transverse field $B_x$. The system of equations (1)–(6) is perturbed and linearized, and the ideal gas law and the magnetic solenoidivity condition are used to eliminate the pressure and one component of the fluctuating magnetic field. All perturbations are of the form $\exp(ik_x x + ik_y y + st)$, where s 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

$$s\mathbf{x}_n = A_n\mathbf{x}_n, \qquad (7)$$

where $\mathbf{x}_n$ is the vector of discretized eigenvectors and $A_n$ is the matrix of discretized differential operators, using a fourth-order 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 $s$ and the corresponding eigenfunctions $\mathbf{x}_n$ using the Linear Algebra PACKage (LAPACK).

## 3 SHEARED MAGNETIC FIELD

In this section, we extend the configuration already considered by Cattaneo et al. (1990a) and Kusano et al. (1998) to the 3D case. We consider two different cases: in the first one, the initial field is unidirectional, say $B_y$, whereas in the second case, a horizontal transverse field $B_x$ is added. For these two different cases, the horizontal magnetic field $B_H$ is the same and we compare both 2D and 3D simulations. By 2D simulation, we mean that all the variables are invariant in one of the horizontal directions (namely in the y direction). The basic state is defined as follows: the initial horizontal

field $B_H(z, t = 0)$ only depends on $z$ and is the same in both cases,

$$B_H(z) = \sqrt{B_x^2 + B_y^2} = \frac{1}{2} \left[ 1.0 + \text{erf} \left( \frac{z - z_t}{\Delta z} \right) \right], \qquad (8)$$

where $z_t$ is the depth corresponding to the top of the unstable layer, whereas $\Delta z$ corresponds to the width of the transition between the layer and the non-magnetized atmosphere above. This choice is compatible with the boundary conditions defined in Section 2. In the following, we choose $z_t = 0.6$ and $\Delta z = 1/60$. The first reference case is obtained by taking $B_x = 0$ so that the field is initially unidirectional, in the y direction. In the sheared case, and following Cattaneo et al. (1990a), the $B_x$ is chosen to be linearly dependent on the depth $z$ according to

$$B_x = \begin{cases} \epsilon (z - z_r) & \text{if } z_t < z < z_r \\ 0 & \text{otherwise} \end{cases} \qquad (9)$$

(this profile is actually smoothed using the error function to avoid spurious numerical effects due to discontinuities). The depth $z_r$ at which the transverse field goes to zero is called the resonant layer by Cattaneo et al. (1990a). We choose here $z_r = 0.8$. The initial amplitude of the transverse field is derived from $\epsilon = -0.75$ (the maximum value of $B_x$ initially occurs at $z = 0.6$ and is $B_x = 0.15$). Note that we tried different configurations with different resonant layers and different amplitudes for the x 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 Fig. 1 a comparison of 2D and 3D simulations for the two cases with and without a transverse field. The cases with a unidirectional field $B_y$ are shown in the upper part of the figure, whereas the cases where $B_x \neq 0$ 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 2D. On the other hand, in the sheared case (i.e. when $B_x \neq 0$ initially, lower line in Fig. 1), a drastic difference is observed between 2D and 3D simulations. In 2D, we recover the large-scale instability already observed by several authors (see for example Cattaneo et al. 1990a; Kusano et al. 1998). The small-scale 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 the $B_y$ field is clear. Fig. 1 also shows 3D visualization comparing the unidirectional and sheared cases using the software VAPOR<fnref n="1" /> (Clyne et al. 2007). It is clear that the nature of the instability is similar in both cases, apart from the change in the direction of the field lines. In 2D, due to geometrical constraints, it is not possible for the x field to undergo interchange instability, whereas it is possible in 3D. The results are similar for many different configurations (e.g. changing $z_t$, $z_r$, $\epsilon$ and the shape of the $B_x$ 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

<footnote n="1">¹ http://www.vapor.ucar.edu/</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Comparison between 2D and 3D simulations. Top line: the initial field is unidirectional, directed perpendicularly to the plane. Bottom line: the initial magnetic layer is sheared. In the left column, we show contours of $B_y$ from the 2D simulations, whereas the middle column corresponds to one particular $(x, z)$ plane inside the 3D domain. In the right column, we show a volume rendering of the magnetic energy from the 3D simulations at the same time. All visualizations are performed at the same time $t \approx 20$, except for the bottom-left figure, which corresponds to $t \approx 50$ (due to the fact that the growth rate of the instability is strongly reduced in this particular case; see, for example, Cattaneo et al. 1990a).</figcaption></figure>

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 in both two and three dimensions.

## 4 INSTABILITY IN A MAGNETIZED ATMOSPHERE

In this section, unlike in Section 3 and previous studies, the atmosphere above the layer of strong field is now also magnetized. Initially, we fix the vertical field to be zero everywhere whereas the horizontal field is only varying in the vertical direction according to

$$B_H = \sqrt{B_x^2 + B_y^2} = 0.2 + 0.5 \left[ 1.0 + \operatorname{erf} \left( \frac{z - z_t}{\Delta z} \right) \right], \quad (10)$$

where $z_t$ is the depth corresponding to the top of the unstable layer and $\Delta z$ corresponds to the width of the transition between the layer and the weakly magnetized atmosphere. As in Section 3, we choose $z_t = 0.6$ and $\Delta z = 1/60$. In the following, the region $0 < z < z_t$ is called the atmosphere, whereas the region $z_t < z < 1$ is called the layer. The $x$ component of the magnetic field is derived from

$$B_x = \frac{\epsilon}{2} \left[ 1.0 - \operatorname{erf} \left( \frac{z - z_t}{\Delta z} \right) \right]. \quad (11)$$

The profile of the horizontal field can be seen in Fig. 2(a), for $\epsilon = 0.15$.

We define the pitch angle as being $\alpha = \operatorname{atan}(B_y/B_x)$ so that $\alpha$ varies between 0 and $\pi/2$, depending on the value of $\epsilon$. 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 $x$ 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 $B_x$ 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 destabilized modes are interchange modes ($k_y = 0$). The instability is localized at the interface between the magnetic layer and the magnetized atmosphere. Fig. 3 shows how the growth rate of the instability varies as a function of the horizontal wavenumber $k_x$ for different pitch angles $\alpha$. 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 $\alpha$. The wavenumber of the most unstable mode 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 $\alpha$. This will be discussed further in the following sections. Although our configuration differs, our conclusion is qualitatively similar to the one drawn by Cattaneo et al. (1990a) and Kusano et al. (1998), i.e. as the intensity of the transverse field is increased, both the growth rate of the instability and the transverse wavenumber decrease.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Initial condition for the horizontal components of the magnetic field, the temperature, *T*, and the density, *ρ*. In the top image, the dashed lines correspond to the vertical profile of $B_x$ and $B_y$ for the particular case $\alpha = \pi/4$. In the bottom image, the dashed lines represent the basic state of temperature and density in the absence of a magnetic field. The parameter *F* has been increased from *F* = 0.1 (the value used in this paper) to *F* = 1 in order to enhance the modifications of the basic field-free polytropic state.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Growth rate *s* of the instability versus horizontal wavenumber $k_x$ for the initial configuration presented in Fig. 2. The results are presented for different angles $\alpha$ from 0 (unidirectional field) to $\pi/2$ (the field lines are perpendicular). We stress the case $\alpha = \pi/4$ as it is the largest value considered numerically in the text.</figcaption></figure>

## 4.2 Non-linear 3D numerical simulations

In this section, we discuss the result from non-linear simulations in 3D. The parameter that we vary here is the initial pitch angle of the field in the atmosphere contained between $z = 0$ and 0.6. We first run a reference simulation where $B_x$ is initially zero everywhere. This corresponds to the classical interchange-type instability already considered in many papers (see, for example, Matthews et al. 1995a; Wissink et al. 2000), 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 $\epsilon$ as defined in equation (11). The parameter $\epsilon$ varies from $\epsilon = 0$ (unidirectional field, $\alpha = 0$) to $\epsilon = 0.15$ ($\alpha = \pi/4$). The reason why we do not consider angles larger than $\pi/4$ will be explained later. For each of these initial configurations, we solve equations (1)–(6) using the numerical method described previously.

We present in Fig. 4 vertical slices in the plane $(x, z)$ where contours of the toroidal field $B_y$ are plotted at different times. Bright colours correspond to large values, whereas dark tones correspond to low values of $B_y$. The upper part of the figure corresponds to the unidirectional field case (i.e. $B_x = 0$ initially). As already observed in Section 3, the initial instability consists of an alternation of upward and downward flows. As these flows increase in amplitude, a secondary instability, of Kelvin–Helmholtz 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 $B_y$. 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 $B_y$ field in the upper part of the domain does not modify the nature of the instability, as can be seen comparing the top-middle panel of Fig. 4 with the top-middle panel of Fig. 1. Note however that, as noted by Vasil & Brummell (2009), 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 (Gough 2007).

Now consider a second simulation where the *x* 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 $\pi/4$ (as depicted in Fig. 2a) in order to illustrate the dynamical effect of the transverse field $B_x$. Initially, as seen in the lower-left panel of Fig. 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 wavenumber in the *x* direction is roughly twice as small in the case $\alpha = \pi/4$ as in the reference case $\alpha = 0$ (see also the magnetic energy spectra in Fig. 9). At later times ($t = 16.6$, lower-middle panel of Fig. 4), it is apparent that the non-linear development of the instability is also modified due to the magnetic tension existing at the interface between the unstable layer and the atmosphere above. Fig. 4 shows an alternation of strong and weak downflows 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
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Evolution of the instability in the (x,z) plane perpendicular to the initial field lines in the layer. In the upper row, the magnetic field lines are unidirectional and directed perpendicularly to the plane of visualization. In the lower row, the field lines for *z* < 0.6 are orientated at 45° with respect to the field lines inside the layer. Bright colours correspond to the strong *B*<sub>y</sub> field, whereas dark tones correspond to the weaker field. Time is increasing from left to right: *t* = 8.3, 16.6 and 37.4.</figcaption></figure>

the unidirectional case. At the final time $t = 37.4$, the difference between the two cases is very clear: the presence of the *x* field in the atmosphere tends to generate large coherent magnetic structures in the *x* direction containing fluctuations at smaller scales generated during the initial stages of the instability. The toroidal field $B_y$ is also much less diffuse than in the unidirectional case, and a clear interface between field and field-free 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 detail at one of the merging events, 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 (x, z) plane is presented in Fig. 5. We

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Details of the merging mechanism. The contours show the current density $j^2$ in a small part of the (x, z) plane for different times. Dark tones correspond to weak current, whereas light colours correspond to strong current. The same colour scale is used for the different times. The velocity field in the plane (x, z) is also shown as arrows whose length depends on the velocity amplitude.</figcaption></figure>

plot the current density $j^2 = (\nabla \times B)^2$ as contours along with arrows representing the velocity field in this particular vertical plane. Initially, a strong downflow in the middle of the panel is detected, dragging down the magnetic field $B_y$. This vertical downflow must also advect the weak $B_x$ field existing in the atmosphere, therefore working against magnetic tension. At time $t = 20.8$ (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 downflow. 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 Fig. 6 a volume rendering of the magnetic energy $B^2/2$. Again, we compare the unidirectional case $\alpha = 0$ with the case $\alpha = \pi/4$ at the same time $t = 24$. Without any initial transverse field $B_x$, the results are similar to what was already observed by Matthews et al. (1995a) and Wissink et al. (2000). The initial linear instability is purely 2D and corresponds to the interchange of horizontal straight field lines. Eventually, the interaction of counter-rotating vortices generated by the secondary instability leads to the arching of the magnetic structures. Note that, although these magnetic structures are clearly observed in Fig. 6, they are very weak and lack twist. The use of isosurfaces, as in Wissink et al. (2000), can be misleading as it gives the impression that 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 small-scale magnetic features. In the $\alpha = \pi/4$ case, however, we observe larger magnetic structures, containing a significant amount of twist. Note that these larger structures are far from being well-defined localized flux tubes (Cattaneo, Brummell & Cline 2006). We however stress the fact that
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Volume rendering of the magnetic energy. The case $\alpha = 0$ is presented on the left, whereas the case $\alpha = \pi/4$ is shown in the middle. The visualizations are realized at the same time $t = 24$. The figure on the right presents a subvolume of the case $\alpha = \pi/4$ where magnetic field lines have been added in order to enhance the twisted nature of the field.</figcaption></figure>

our polytropic atmosphere is stably stratified and the magnetic tension resulting from the presence of the oblique $B_x$ 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 localized flux tubes from the break-up of a magnetic layer. The last panel on the right in Fig. 6 shows a close-up with several magnetic field lines in blue. As the strong $B_y$ 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 Choudhuri (2003), where the twist is generated by the rise of strong toroidal fields inside the convective zone containing a poloidal field. Note however that Choudhuri (2003) considered already formed flux tubes and not a uniform layer of magnetic field.

We now discuss more quantitatively the effect of the transverse field $B_x$ on the development of the instability. Fig. 7 shows the evolution with time of several global quantities. In the following, the brackets $\langle \cdot \rangle$ denote a spatial average over all coordinates. In particular, we define the total kinetic energy $\langle \rho u^2 \rangle / 2$ and the total magnetic energy $\langle B^2 \rangle / 2$. As previously, we compare the results from simulations without the initial $B_x$ field ($\alpha = 0$) 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 in the left-hand panel of Fig. 7. Note that, as predicted by the linear analysis of Section 4.1, the growth rate decreases when $\alpha$ increases.

After approximately 10 sound crossing times, the instability saturates. The evolution after this stage strongly depends on $\alpha$. For the case $\alpha = 0$, 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 non-linear regime. For $\alpha = \pi/4$, 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 Kelvin-Helmholtz instability. It is clear that the $x$ field, if strong enough, will tend to damp this generation of vorticity. In order to create strong down- and upflows, it is now necessary to work against the tension of the weak magnetic field $B_x$. The evolution of the total enstrophy $\langle \omega^2 \rangle$ in the numerical domain is plotted versus time in Fig. 7, where $\omega = \nabla \times u$. In all cases, we observe a rapid increase in 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 $\alpha \neq 0$ than in the case where $\alpha = 0$. This leads to a reduction of the viscous dissipation, therefore reducing the rate at which the kinetic energy is dissipated.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Evolution of the total kinetic energy, enstrophy, magnetic energy and current density versus time for four different angles from $\alpha = 0$ to $\pi/4$.</figcaption></figure>
The two other panels in Fig. 7 show the evolution with time of the total magnetic energy and current density $\langle j^2 \rangle$, where $j = \nabla \times B$. In all cases, the magnetic energy decays whereas the current density increases with time. As $\alpha$ 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 $\alpha \neq 0$, as expected from the initial conditions.

Following Stone & Gardiner (2007b), 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 Rayleigh–Taylor instability, we define the fraction of strong horizontal field as

$$f_s = \frac{B_H - B_H(z = 0)}{B_H(z = 1) - B_H(z = 0)}, \qquad (12)$$

where $B_H = \sqrt{B_x^2 + B_y^2}$. $B_H(z = 0)$ corresponds to the initial amplitude of the horizontal field at the upper boundary and can be considered as the weak field. $B_H(z = 1)$ is the initial horizontal field at the lower boundary and is therefore the strong field. The fraction of weak field is simply $f_w = 1 - f_s$. The amount of mixing in the system can be estimated by the quantity

$$\Theta(z) = 4\langle f_s f_w \rangle, \qquad (13)$$

which is equal to 0 when the averaged horizontal field is equal to its initial value at the vertical boundaries and 1 if the field is perfectly mixed and is equal to the average between these two values.

Fig. 8 shows the mixing parameter $\Theta$ as a function of $z$ for the four cases considered at $t = 30$. The thin black line corresponds to the value of $\Theta$ 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 $\alpha$ increases. This result can be explained by the reduced growth rate of the initial linear instability (see Fig. 3) in the presence of a transverse field and by the damping of the secondary instability leading to weaker vortical motions (see Fig. 7). In the unidirectional case, the potential energy is released to counter the stable stratification, leading to the rise in the toroidal magnetic field. When $\alpha \neq 0$, 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 Fig. 8 also shows that the lower front is roughly unchanged by the presence of a transverse field.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Mixing parameter $\Theta(z)$ as defined by equation (13) versus depth at time $t = 30$. The thin black line corresponds to the initial profile at $t = 0$.</figcaption></figure>

The generation of large-scale coherent magnetic structures is clearly visible in the magnetic energy spectra defined by

$$E_M(k_x) = \frac{1}{2} \sum_z \sum_y \hat{B}(k_x, y, z) \cdot \hat{B}^*(k_x, y, z), \qquad (14)$$

where $\hat{B}(k_x, y, z)$ is the 1D Fourier transform in the $x$ direction of $B(x, y, z)$, whereas the asterisk 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 $y$ direction are rather similar between the different cases).

Fig. 9 shows the magnetic energy spectrum for $\alpha = 0$ and $\pi/4$, at four different times spanning from the initial linear instability ($t = 2$) to the late stage of the non-linear phase ($t = 32$). Initially, the most unstable horizontal wavenumber corresponds to the linear results already discussed in Section 4.1. The most unstable wavenumber is roughly $k_x/2\pi \approx 14$ for $\alpha = 0$ and $k_x/2\pi \approx 8$ for $\alpha = \pi/4$, which can be compared to the results presented in Fig. 3. The time evolution of the magnetic energy spectra can be described in two phases. First, the small-scale instability grows exponentially from $t \approx 2$ to $t \approx 12$. The effect of the transverse field is to reduce the initial linear growth rate, as already discussed in Figs 3 and 7. From $t \approx 12$, 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 towards small

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Magnetic energy spectrum $E_M(k_x)$ as defined by equation (14) at four different times $t = 2$ (light tones), $t = 12.5$, $t = 24$ and $t = 32$ (dark tones). (a) $\alpha = 0$ and (b) $\alpha = \pi/4$. The results are obtained by spatially averaging over $y$ and $z$.</figcaption></figure>
scales, as we observe steeper energy spectra for $\alpha = \pi/4$ than for $\alpha = 0$, 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 $\alpha = 0$. In contrast, the magnetic energy is still growing at large scales for $\alpha = \pi/4$. The effect of the transverse field is indeed clearly visible in the small wavenumber region $k_x/2\pi < 3$. For $\alpha = \pi/4$, a clear peak is observed at large scales, which corresponds to the coherent magnetic structures already observed in Figs 4 and 6. This is also consistent with the larger values of the global magnetic energy observed in Fig. 7 when $\alpha \neq 0$. Without a transverse field, the magnetic energy peaks at larger horizontal wavenumbers. Note that this difference cannot be attributed to linear effects, as the most unstable wavenumbers predicted by the linear analysis are much larger (see Fig. 3). During the non-linear phase, the small-scale features (i.e. $k_x/2\pi > 4$) are roughly independent of $\alpha$, 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 fully developed turbulent state.

Finally, the amount of twist in the magnetic field lines, as observed in Fig. 6, can be estimated looking at the current helicity, defined by $\langle j \cdot B \rangle$. Since both magnetic energy and current density are time dependent (see Fig. 7), we consider the relative current helicity defined by

$$ \mathcal{H}_B = \frac{\langle j \cdot B \rangle}{\sqrt{\langle j^2 \rangle} \sqrt{\langle B^2 \rangle}}. \qquad (15) $$

Note that given our initial configuration, the amount of current helicity depends on $\alpha$. We show in Fig. 10 the evolution with time of $\mathcal{H}_B$ for different values of $\alpha$. The effect of the initial condition is clearly apparent as $\mathcal{H}_B = 0$ when $\alpha = 0$, whereas $\mathcal{H}_B$ initially increases with $\alpha$. Without the 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 $\alpha \neq 0$, we observe an increase in the relative current helicity from its small initial value. This is particularly visible in the $\alpha = \pi/4$ case. After a transient stage where $\mathcal{H}_B$ increases and then decreases, we observe a monotonic increase in the relative current helicity up to the end of the simulation. Note that $\mathcal{H}_B$ continues to increase even after the saturation of the initial linear phase, which roughly corresponds to $t \approx 20$ (see Fig. 7). Again, since the atmosphere is stably stratified, we expect the magnetic field to ultimately stop

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Relative current helicity, $\mathcal{H}_B$, as defined by equation (15) versus time.</figcaption></figure>

rising and the current helicity to ultimately decay. If the entropy gradient were nearly negligible, as expected in the convective zone, the magnetic structures observed in our simple model could continue to rise unimpeded, eventually becoming kink unstable (Fan et al. 1999). 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 Choudhuri (2003). 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 depth dependent. 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 large-scale coherent magnetic structures from the buoyancy instability of a sheared magnetic layer, as studied by Cattaneo et al. (1990a) and Kusano et al. (1998), seemed to be limited to 2D 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 magnetized atmosphere above. The weakly magnetized atmosphere was justified by the presence of convection, which is expected to contain non-negligible magnetic fluctuations, both due to the redistribution of large-scale magnetic fields and due to local small-scale 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 non-linear phase of the instability, the small-scale magnetic structures merged together to overcome magnetic tension, leading to the formation of magnetic structures at 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 that a significant amount of twist was generated due to the rise of the toroidal field in the weakly magnetized 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 small-scale 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 small-scale magnetic fluctuations.

There are two stabilizing effects in the present model, the stability 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 to or larger than $\pi/4$. The stabilizing 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 small-scale magnetic structures. Locally these structures could have a similar effect to the global tension considered in our simulations, although we did not 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 to the rise and formation of local coherent twisted magnetic structures. Another possibility would be to extend the work of Fan et al. (2003) and Jouve & Brun (2009), 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 small-scale 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 (Pinto & Brun, private communication).

We conclude by arguing that the merging mechanism illustrated by this simple model could have implications on the formation of large-scale magnetic structures from a deep-seated 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 magnetized convective zone.

## ACKNOWLEDGMENTS

This work has been financially supported by STFC. CPU time was provided by the HPC resources of CALMIP under the allocation 2012-P1118 and by the UKMHD supercomputing facility located in Warwick.

## REFERENCES

* Abbett W. P., Fisher G. H., Fan Y., Bercik D. J., 2004, ApJ, 612, 557

* Bushby P. J., Houghton S. M., Proctor M. R. E., Weiss N. O., 2008, MNRAS, 387, 698

* Caligari P., Moreno-Insertis F., Schüssler M., 1995, ApJ, 441, 886

* Cattaneo F., Hughes D. W., 1988, J. Fluid Mech., 196, 6

* Cattaneo F., Chiueh T., Hughes D. W., 1990a, J. Fluid Mech., 219, 1

* Cattaneo F., Chiueh T., Hughes D. W., 1990b, MNRAS, 247, 1

* Cattaneo F., Brummell N. H., Cline K. S., 2006, MNRAS, 365, 727

* Chae J., Moon Y. J., 2005, ApJ, 701, 1110

* Chatterjee P., Choudhuri A. R., Petrovay K., 2006, A&A, 449, 781

* Choudhuri A. R., 1987, ApJ, 316, 788

* Choudhuri A. R., 2003, Sol. Phys., 215, 31

* Clyne J., Mininni P., Norton A., Rast M., 2007, New J. Phys., 9, 301

* D'Silva S., Choudhuri A. R., 1993, A&A, 272, 621

* DeVore C. R., 2000, ApJ, 539, 944

* Emonet T., Moreno-Insertis F., 1996, ApJ, 492, 804

* Fan Y., 2001, ApJ, 546, 509

* Fan Y., Gong D., 2000, Sol. Phys., 192, 141

* Fan Y., Fisher G. H., McClymont A. N., 1994, ApJ, 436, 907

* Fan Y., Zweibel E. G., Lantz S. R., 1998, ApJ, 493, 480

* Fan Y., Zweibel E. G., Linton M. G., Fisher G. H., 1999, ApJ, 521, 460

* Fan Y., Abbett W. P., Fisher G. H., 2003, ApJ, 582, 1206

* Favier B., Bushby P. J., 2012, J. Fluid. Mech., 690, 262

* Gough D. O., 2007, in Hughes D. W., Rosner R., Weiss N. O., eds, The Solar Tachocline. Cambridge Univ. Press, Cambridge, p. 3

* Guerrero G., Käpylä P. J., 2011, A&A, 533, A40

* Holder Z. A., Canfield R. C., McMullen R. A., Nandy D., 2004, ApJ, 611, 1149

* Jouve L., Brun A. S., 2007, Astron. Nachr., 328, 1104

* Jouve L., Brun A. S., 2009, ApJ, 629, 1300

* Kersalé E., Hughes D. W., Tobias S. M., 2007, ApJ, 663, L113

* Kusano K., Moriyama K., Miyoshi T., 1998, Phys. Plasmas, 5, 2582

* Leka K. D., Canfield R. C., McClymont A. N., van Driel Gesztelyi L., 1996, ApJ, 462, 547

* Longcope D. W., Fisher G. H., Arendt S., 1996, ApJ, 464, 999

* Longcope D. W., Fisher G. H., Pevtsov A. A., 1998, ApJ, 507, 417

* Matthews P. C., Hughes D. W., Proctor M. R. E., 1995a, ApJ, 448, 938

* Matthews P. C., Proctor M. R. E., Weiss N. O., 1995b, J. Fluid Mech., 305, 281

* Moreno-Insertis F., 1983, A&A, 122, 241

* Nozawa S., 2005, PASJ, 57, 995

* Parker E. N., 1955, ApJ, 121, 491

* Pevtsov A. A., Canfield R. C., Latushko S. M., 2001, ApJ, 549, L261

* Schüssler M., 1979, A&A, 71, 79

* Seehafer N., 1990, Sol. Phys., 125, 219

* Silvers L. J., Vasil G. M., Brummell N. H., Proctor M. R. E., 2009, ApJ, 702, L14

* Silvers L. J., Vasil G. M., Brummell N. H., Proctor M. R. E., 2010, in Brummell N. H., Brun A. S., eds, Proc. IAU Symp. 271, Astrophysical Dynamics: From Stars to Galaxies. Cambridge Univ. Press, Cambridge, p. 111, doi:10.1017/S1743921311017649

* Spruit H. C., 1981, A&A, 102, 129

* Stone J. M., Gardiner T., 2007a, ApJ, 671, 1726

* Stone J. M., Gardiner T., 2007b, Phys. Fluids, 19, 094104

* Tobias S. M., Brummell N. H., Clune T. L., Toomre J. T., 2001, ApJ, 549, 1183

* Török T., Kliem B., 2005, ApJ, 630, L97

* Vasil G. M., Brummell N. H., 2009, ApJ, 690, 783

* Wissink J. G., Hughes D. W., Matthews P. C., Proctor M. R. E., 2000, MNRAS, 318, 501

This paper has been typeset from a TeX/Latex file prepared by the author.