BIPOLAR MAGNETIC STRUCTURES DRIVEN BY STRATIFIED TURBULENCE WITH A CORONAL ENVELOPE

Sunspots are visible surface manifestations of magnetic fields inside the Sun. The magnetic fields in sunspots are sufficiently strong to suppress the transport of heat via convective motions and therefore they appear as dark on the solar disk. We are now able to study them with high-resolution telescopes (Scharmer et al. 2011) as well as with realistic numerical simulations (Rempel et al. 2009) to discover new properties of sunspots. However, their formation mechanism is still a subject of active discussion.

It is broadly believed that sunspots are caused by emerging flux tubes (Parker 1955) that are generated at the bottom of the convection zone (Parker 1975). The magnetic flux tubes become unstable and rise until they pierce the solar surface to form bipolar regions and pairs of sunspots (Caligari et al. 1995). The rise of flux tubes and the decay of the resulting active regions is also an ingredient in some dynamo theories to explain the 22 yr cyclic behavior of the solar magnetic field (e.g., Nandy & Choudhuri 2001; Choudhuri 2003). However, these theories assume the existence of thin, isolated magnetic flux tubes at the bottom of the convection zone with fields of the order of 10⁵ G. In direct numerical simulations (DNS) of solar dynamos, thin flux tubes of comparable strength have not yet been found (e.g., Guerrero & Käpylä 2011; Nelson et al. 2013; Käpylä et al. 2013). There is no conclusive evidence of rising magnetic flux tubes from helioseismology either; see Birch et al. (2010, 2013), who also comment on a detection reported by Ilonidis et al. (2011). These difficulties have led to an alternative approach in which sunspots are formed as shallow phenomena near the solar surface (Brandenburg 2005).

The formation of bipolar regions has recently been found in realistic radiation-magnetohydrodynamics (MHD) simulations of near-surface convection (Stein & Nordlund 2012), were kG magnetic fields were inserted at the bottom of a 20 Mm deep box. The authors argue that deep downdrafts associated with the supergranulation concentrate the magnetic field and determine thereby the scale of their separation. They also emphasize the importance of radiative cooling at the surface.

Two alternative mechanisms for the spontaneous formation of flux concentrations have been discussed in the literature. One is based on a turbulent thermo-magnetic instability in small-scale turbulence with radiative boundaries, where the magnetic quenching of the turbulent convective heat transport is used in mean-field simulations (MFSs) to produce flux concentrations (Kitchatinov & Mazur 2000). The other mechanism is based on the suppression of the total turbulent pressure (the sum of hydrodynamic and magnetic contributions) by the magnetic field. This leads to the so-called negative effective magnetic pressure instability (NEMPI), which is another mechanism that forms flux concentrations in a stratified turbulent medium (see, e.g., Brandenburg et al. 2011; Kemel et al. 2012). The original idea of NEMPI goes back to early work by Kleeorin et al. (1989, 1990) and has led to numerous studies both analytically (Kleeorin & Rogachevskii 1994; Kleeorin et al. 1996; Rogachevskii & Kleeorin 2007) and through DNS and MFS (Brandenburg et al. 2010, 2012; Losada et al. 2012, 2013a, 2013b; Jabbari et al. 2013).

More recently, Brandenburg et al. (2013b) found super-equipartition magnetic field concentrations in DNS by imposing a relatively weak vertical magnetic field on forced strongly stratified turbulence. Their results may be related to earlier simulations of turbulent convection in which a vertical magnetic field was found to segregate into magnetized and unmagnetized regions (Tao et al. 1998), but it is still unclear whether those results are also caused by NEMPI or by some other mechanism that still remains to be understood. However, neither the turbulent thermo-magnetic instability nor NEMPI have yet been able to produce bipolar magnetic regions. Except for the work of Stein & Nordlund (2012), bipolar regions have previously been studied by advecting a semi-torus shaped twisted flux tube of 9 kG through the bottom boundary at a 7.5 Mm depth of thermally relaxed convection (Cheung et al. 2010).

In the current study we demonstrate that bipolar regions can be obtained naturally in DNS in the presence of an outer coronal layer. Indeed, all previous simulations of NEMPI and the turbulent thermo-magnetic instability have a rigid upper boundary, which is clearly unrealistic. In this Letter, we alleviate the constraints from a rigid upper boundary by combining NEMPI in a forced turbulent layer with a coronal envelope at the top. This approach has been successful in generating plasmoids reminiscent of coronal mass ejections (Warnecke et al. 2011). In that work, a simplified corona is combined with a large-scale dynamo that is either generated by forced turbulence in a Cartesian domain (Warnecke & Brandenburg 2010) or in spherical coordinates (Warnecke et al. 2011), or through self-consistent turbulent convection (Warnecke et al. 2012b). This approach leads to a more realistic "boundary condition" at the interface between the dynamo and corona without using, for example, radiative transfer. Not only has it led to coronal ejections, but also to the discovery of a change of sign of magnetic helicity some distance away from the star (Warnecke et al. 2012a), and the realization that both the dynamo and the differential rotation may be affected by the presence of a free boundary as modeled by the coronal envelope (Warnecke et al. 2013).

For a more realistic model, several other factors need to be accounted for. Except for radiation and ionization that would potentially lead to new effects (see, e.g., Barekat & Brandenburg 2013), which might distort our interpretation, there would be the issue of the origin of the magnetic field itself. It should really come from a self-consistent large-scale dynamo beneath the surface. Again, this leads to new complications, through its interaction with NEMPI, which have only recently been addressed by Jabbari et al. (2013) and Losada et al. (2013a). To avoid those complications, we take a horizontal magnetic field as the initial condition. We mainly consider the case of a uniform magnetic field that is then also present in the corona at all times (case A), and compare with a case in which it is initially only in the turbulence layer, but slowly decaying (case B). These fields are weak, but they might have interesting unexpected and unexplored effects that will be the topic of the current Letter. At this point we cannot be certain that it has direct applications to the Sun, but in view of the fact that research into NEMPI has been growing rapidly, it is important to explore all aspects of it before we can make definitive statements about its applicability to sunspots.

Our model is similar to that of Brandenburg et al. (2013b), where they solve the isothermal hydromagnetic equations in the presence of vertical gravity g = (0, 0, − g), but here the random nonhelical forcing f acts just in the lower part (−π < z < 0) of a Cartesian domain of size L_x × L_y × L_z, where L_x = L_y = 2π and L_z = 3π. Also, instead of a vertical field, we now take a weak uniform horizontal magnetic field, either as an imposed one throughout the domain, B_imp = (0, B₀, 0) with A = 0 (case A), or as an initial one in the lower part (z < 0) by putting A_x = −max (0, −z) B₀ and B_imp = 0 (case B). Similar to Warnecke & Brandenburg (2010), we employ a forcing function f(x, t) that is modulated by

$$\begin{equation} \theta _w(z)={\textstyle {1\over 2}}\left(1-{\rm erf}{z\over w}\right), \end{equation} \tag{ 1 }$$

where w is the width of the transition. We thus solve the equations for the velocity u, the magnetic vector potential A, and the density ρ:

$$\begin{equation} {{D} {}{\boldsymbol {u}} {}\over {D} {}t}={\boldsymbol {g}} {}+\theta _w(z){\boldsymbol {f}} {}+{1\over \rho } \left[-c_{\rm s}^2{\boldsymbol {\nabla }} {}\rho +{\boldsymbol {J}} {}\times {\boldsymbol {B}} {}+{\boldsymbol {\nabla }} {}\cdot (2\nu \rho {\boldsymbol {{\sf S}}} {})\right],\quad \end{equation} \tag{ 2 }$$

$$\begin{equation} {\partial {\boldsymbol {A}} {}\over \partial t}={\boldsymbol {u}} {}\times {\boldsymbol {B}} {}+\eta \nabla ^2{\boldsymbol {A}} {}, \end{equation} \tag{ 3 }$$

$$\begin{equation} {{D} {}\ln \rho \over {D} {}t}=-{\boldsymbol {\nabla }} {}\cdot {\boldsymbol {u}} {}, \end{equation} \tag{ 4 }$$

where ν is the kinematic viscosity, η is the magnetic diffusivity, ${\boldsymbol {B}} {}={\boldsymbol {B}} {}_{\rm imp}+{\boldsymbol {\nabla }} {}\times {\boldsymbol {A}} {}$ is the magnetic field, ${\boldsymbol {J}} {}={\boldsymbol {\nabla }} {}\times {\boldsymbol {B}} {}/\mu _0$ is the current density, μ₀ is the vacuum permeability, ${\sf S}_{ij}={\textstyle ({1/2})} (u_{i,j}+u_{j,i})-{\textstyle ({1/3})} \delta _{ij}{\boldsymbol {\nabla }} {}\cdot {\boldsymbol {u}} {}$ is the traceless rate of strain tensor, and commas denote partial differentiation. The forcing function f consists of random, white-in-time, plane non-polarized waves (Haugen et al. 2004) with an average wavenumber k_f = 30 k₁, where k₁ = 2π/L_x is the lowest wavenumber in the domain in the horizontal direction. The forcing strength is such that the turbulent rms velocity is approximately independent of z with $u_{\rm rms}=\langle {\boldsymbol {u}} {}^2\rangle ^{1/2}\approx 0.1\,c_{\rm s}$. The value of the gravitational acceleration g is related to the density scale height $H_\rho =c_{\rm s}^2/g$ and is chosen such that k₁H_ρ = 1, which leads to a density contrast between bottom and top of exp (3π) ≈ 1.2 × 10⁴.

Except for the profile of the forcing function and the extent of the domain, our setup is also similar to Kemel et al. (2012), who used the same values of the fluid Reynolds number  Re ≡ u_rms/νk_f = 38, the magnetic Prandtl number  Pr_M ≡ ν/η ≡ 1/2, and thus of the magnetic Reynolds number  Re_M ≡  Re  Pr_M = 19, which is known to lead to negative effective magnetic pressure for the mean magnetic fields, $\overline{\boldsymbol {B}}{}$, in the range $0<|\overline{\boldsymbol {B}}{}|/B_{\rm eq}<0.4$. Here, $\overline{\boldsymbol {B}}{}$ is obtained by averaging over the scale of several turbulent eddies. The magnetic field is expressed in units of the local equipartition field strength, $B_{\rm eq}=\sqrt{\mu _0\rho } \, u_{\rm rms}$, while the initial magnetic field B₀ is specified in units of the value at z = 0, namely, $B_{\rm eq0}=\sqrt{\mu _0\rho _0} \, u_{\rm rms}$, where ρ₀ = ρ(z = 0). We use B₀/B_eq0 = 0.02, which is also the field strength used in the main run of Brandenburg et al. (2013b). Time is expressed in turbulent-diffusive times, $\tau _{\rm td}=(\eta _{\rm t0}k_1^2)^{-1}$, where η_t0 = u_rms/3k_f is the estimated turbulent magnetic diffusivity.

The simulations are performed with the Pencil Code,⁵ which uses sixth-order explicit finite differences in space and a third-order accurate time-stepping method. We use a resolution of 256 × 256 × 512 mesh points in the x-, y-, and z-directions. We adopt periodic boundary conditions in the xy plane and present our results by shifting our coordinate system such that the regions of interest lie at the center around x = y = 0. On z = −π we apply a stress-free perfect conductor condition and on z = 2π a stress-free vertical field condition.

We report the spontaneous formation and decay of bipolar magnetic regions at the surface (z = 0), which is the boundary between regions with and without forcing. These two parts resemble the upper convection zone and a simplified corona of the Sun. In Figure 1, we show for case A the bipolar region as the normalized vertical magnetic field B_z/B_eq (left panel) and the normalized magnetic energy ${\boldsymbol {B}} {}^2/B_{\rm eq}^2$ (right panel) at the moment of maximum strength, t/τ_td = 2. Note that the y-direction points to the right, while the positive x-direction points downward, so the coordinate system has been rotated by 90° to allow a view that is more similar to how bipolar regions are oriented on the solar disk. The shapes of both magnetic spots are nearly circular, but are still disturbed by the turbulent motion acting on the magnetic field. The field outside this bipolar region is weak: almost all the magnetic field is concentrated inside the bipolar region. At and slightly above z = 0, we find field strengths significantly above the equipartition value. This is seen more clearly in Figure 2(a), where we show for case A profiles of B_z(y)/B_eq through x = 0 at three heights. We normalize the magnetic field by its local equipartition value, B_eq(z) (Figure 2a) or by the strength of the initial field B₀ (Figure 2b). The field B_eq(z) is also shown as a thick line in Figure 2(b). At the surface we have B_eq(0)/B₀ ≈ 50, but for z > 0 it drops sharply to values below 20. At each height, we have computed maximum and minimum field strengths of the mean field, $\overline{\boldsymbol {B}}{}$, defined here by Fourier filtering to include only wavenumbers below k_f/2, as functions of z, i.e., $\overline{B}_z^{\max }(z)$ and $-\overline{B}_z^{\min }(z)$, respectively. It turns out that $\overline{B}_z^{\min }(z)\approx -\overline{B}_z^{\max }(z)$, and that both decline more slowly with height than B_eq, so the field reaches super-equipartition strength in the outer parts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We recall that in case A, we have B₀ = B_imp over the whole domain. However, the field quickly becomes tangled by the random velocity field in z < 0 where the forcing is acting to produce small-scale magnetic fields on the scale of the turbulence. In the upper part, however, this horizontal field stays roughly unchanged up to the instant when it becomes affected by a large-scale instability (NEMPI). As shown in Figure 2(c), the rms values of the three components of the magnetic field at the surface (z = 0) grow rapidly until they saturate at t/τ_td ≈ 0.2. At t/τ_td ≈ 1, the magnetic field has attained a strong vertical component while the horizontal one declines. By the time t/τ_td = 2, the vertical field is stronger than the horizontal field until all three components decay rapidly to a lower value and saturate there.

To see whether the increase of vertical magnetic field and structure formation is related to NEMPI, we show in Figure 2(d) that the effective magnetic pressure ${\mathcal P}_{\rm eff}$ is indeed negative below the surface in the turbulent region for case A. It can be calculated using an approach applied by Brandenburg et al. (2012):

$$\begin{equation} {\mathcal P}_{\rm eff}={\overline{\boldsymbol {B}}{}^2\over 2B_{\rm eq}^2}+{\left(\Pi _{xx}^{(B)}-\Pi _{xx}^{(0)} +\Pi _{zz}^{(B)}-\Pi _{zz}^{(0)}\right) \over B_{\rm eq}^2}, \end{equation} \tag{ 5 }$$

where $\Pi _{ij}^{(B)}\equiv \overline{\rho u_iu_j} +{\textstyle ({1/2})}\delta _{ij}\mu _0^{-1}\overline{{\boldsymbol {b}} {}^2}-\mu _0^{-1}\overline{b_ib_j}$ and $\Pi _{ij}^{(0)}$ are the components of the total (Reynolds and Maxwell) turbulent stress tensor with and without mean magnetic influence, respectively, and δ_ij is the unit Kronecker tensor.

Horizontal cross-sections of B_z(x, y)/B_eq through z = 0 are shown in Figure 3 for case A at different times. At t/τ_td = 1, structures begin to form that become more coherent and more circular while decreasing their distance to a minimum until they lie directly next to each other (t/τ_td = 2). This is also the time of maximum field strength and maximum coherence, and agrees with the peak of $B_z^{\rm rms}$ in Figure 2(c). After this time, the distance and field strength of the two polarities decrease until no large-scale structures are visible (t/τ_td = 3.5).

Zoom In Zoom Out Reset image size

In Figure 4 we compare yz slices through the bipolar region for cases A and B. The vertical magnetic field is color coded and arrows indicate magnetic field vectors in the plane. The field in the region of z > 3 is not strongly disturbed by the structure formation and represents B_imp in case A and is ≈0 in case B. The weaker field concentration in case B might be related to the decay of the initial field. In both cases, the bipolar spot orientation is peculiar because an emerging flux tube with similar field direction as the imposed field would cause an inverted vertical flux configuration.

Zoom In Zoom Out Reset image size

The formation of magnetic structures can be caused by the negative contribution of turbulence to the effective large-scale magnetic pressure (the sum of turbulent and non-turbulent contributions). For large magnetic Reynolds numbers, the turbulent contributions are larger than the non-turbulent ones, and the effective magnetic pressure becomes negative; see Figure 2(d). This results in a large-scale instability, which causes a redistribution of mass so that a large-scale flow is generated. This flow may also drive the magnetic field patches together. Since turbulence has produced similar strengths for all three components of the magnetic field, and since NEMPI allows for stronger vertical fields than horizontal ones (Brandenburg et al. 2013b), the result is the formation of strong vertical field structures.

It is important to realize that our setup corresponds to an initial value problem in the sense that the magnetic field affects the effective magnetic pressure. It changes the horizontally symmetric background state, which is unstable with respect to NEMPI. This leads to the formation of magnetic structures that tend to stabilize the system. This is the reason why, with our present setup, a bipolar magnetic region occurs only once. Of course, if we apply this mechanism to the Sun, the imposed magnetic field would be provided by a dynamo acting in the convection zone, which would certainly show cycles and fluctuations.

In the coronal region, magnetic loops form that connect the far ends of spots of opposite polarity, as can be seen by plotting 1.5 times the full periodic length in the y-direction. In Figure 5 we show for case A the temporal evolution of the normalized magnetic energy density ${\boldsymbol {B}} {}^2 /B_{\rm eq}^2$ in such a yz slice. The large-scale loop is strongest at t/τ_td ≈ 2 when the two magnetic spots are most concentrated. An exploratory run with a larger extent in the y-direction resulted in similarly oriented bipolar structures with somewhat larger separation.

Zoom In Zoom Out Reset image size

The current study has shown for the first time that a bipolar magnetic region can emerge and later decay as a natural consequence of stratified turbulence with a coronal envelope and a sufficiently large domain size. For the formation process of the bipolar region, radiative cooling and a particular large-scale flow topology as in Stein & Nordlund (2012) seem to be unimportant. However, they will certainly play a role in determining the detailed structure of sunspots (Rempel et al. 2009; Cheung et al. 2010). The interaction with a weak pre-existing coronal field is not unrealistic. A similar interaction has previously been employed in connection with the production of X-ray jets and flares (Isobe et al. 2005). It supports the formation of a bipolar region (case A), but is otherwise not critical (case B). On the other hand, preliminary calculations with a potential field boundary condition and no corona have not yet resulted in bipolar regions.

We find the two polarities approaching each other in the forming process and separating in the decaying phase. This is an important aspect that was not previously expected from turbulent processes (Spruit 2012). The observed formation of the bipolar regions in cases A and B is consistent with NEMPI, which predicts that, even though the optimal field strength for the excitation of NEMPI is significantly below the equipartition value, the magnetic field strengths inside the bipolar region can be significantly above the equipartition value of the turbulence at and slightly above the surface (Brandenburg et al. 2013a).

Useful extensions of the model include more realistic stratifications such as a polytropic one in the lower turbulent part (Losada et al. 2013b) combined with a sharp drop in temperature in the upper layer caused by ionization and radiation. Furthermore, above the chromosphere the temperature should increase again in the solar corona (cf. Warnecke et al. 2013).

We thank the referee for useful comments and acknowledge the computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and the High Performance Computing Center North in Umeå. This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952, by the Swedish Research Council under the project grants 621-2011-5076 and 2012-5797, by EU COST Action MP0806, by the European Research Council under the Atmospheric Research Project No. 227915, and by a grant from the Government of the Russian Federation under contract No. 11.G34.31.0048.