A new look at sunspot formation using theory and observations

TL;DR

This paper explores a novel theoretical approach to sunspot formation based on the negative effective magnetic pressure instability, supported by numerical simulations and helioseismic observations, challenging traditional flux tube models.

Contribution

It introduces the negative effective magnetic pressure instability as a new mechanism for sunspot formation, supported by simulations and helioseismic analysis.

Findings

Numerical simulations produce sunspot-like flux concentrations.

Sunspots may form near the surface due to this instability.

Helioseismic methods can help understand sunspot pre-emergence.

Abstract

Sunspots are of basic interest in the study of the Sun. Their relevance ranges from them being an activity indicator of magnetic fields to being the place where coronal mass ejections and flares erupt. They are therefore also an important ingredient of space weather. Their formation, however, is still an unresolved problem in solar physics. Observations utilize just 2D surface information near the spot, but it is debatable how to infer deep structures and properties from local helioseismology. For a long time, it was believed that flux tubes rising from the bottom of the convection zone are the origin of the bipolar sunspot structure seen on the solar surface. However, this theory has been challenged, in particular recently by new surface observation, helioseismic inversions, and numerical models of convective dynamos. In this article we discuss another theoretical approach to the formation of sunspots: the negative effective magnetic pressure instability. This is a large-scale instability, in which the total (kinetic plus magnetic) turbulent pressure can be suppressed in the presence of a weak large-scale magnetic field, leading to a converging downflow, which eventually concentrates the magnetic field within it. Numerical simulations of forced stratified turbulence have been able to produce strong super-equipartition flux concentrations, similar to sunspots at the solar surface. In this framework, sunspots would only form close to the surface due to the instability constraints on stratification and rotation. Additionally, we present some ideas from local helioseismology, where we plan to use the Hankel analysis to study the pre-emergence phase of a sunspot and to constrain its deep structure and formation mechanism.

Tables

Table 1. Table 1: Main characteristics of the presented paradigms on the formation of sunspots

	Traditional approach	Alternative approach
Dynamo type	tachocline flux-transport	distributed
Role of sunspots	part of dynamo	natural outcome
Origin of sunspots	rising flux tubes	turbulence & stratification
Sunspot formation depth	tachocline	near the surface
Turbulent convection	sunspot decay	sunspot formation

Equations

$$\overline{\mbox{\boldmath$u$}{}\times\mbox{\boldmath$b$}{}}=\alpha\overline{\mbox{\boldmath$B$}}{}-\eta_{t}\overline{\mbox{\boldmath$J$}}{}.$$

$${\partial\overline{\mbox{\boldmath$U$}}{}\over\partial t}+\overline{\mbox{\boldmath$U$}}{}\cdot\mbox{\boldmath$\nabla$}{}\overline{\mbox{\boldmath$U$}}{}=-{1\over\overline{\rho}}\mbox{\boldmath$\nabla$}{}\overline{P}_{\rm gas}+\mbox{\boldmath$g$}{}+\overline{\mbox{\boldmath${\cal F}$}}{}_{\rm M}+\overline{\mbox{\boldmath${\cal F}$}}{}_{\rm K}.$$

$$P_{\rm Tot}=P_{\rm gas}+{\overline{\mbox{\boldmath$B$}}{}^{2}\over 2\mu_{0}}+p_{\rm turb}\equiv\tilde{P}_{\rm gas}+{\cal P}_{\rm eff},$$

$$p_{\mathrm{turb}}=\frac{E_{\mathrm{M}}}{3}+\frac{2E_{\mathrm{K}}}{3}=\frac{2}{3}(E_{\mathrm{K}}+E_{\mathrm{M}})-\frac{1}{3}{E}_{\mathrm{M}}.$$

$${\cal P}_{\rm eff}=\left[1-q_{\rm p}(\overline{\mbox{\boldmath$B$}}{})\right]{\overline{\mbox{\boldmath$B$}}{}^{2}\over 2\mu_{0}}.$$

$$P_{\rm Tot}=\tilde{P}_{\rm gas}+\left(1-q_{\rm p}(\overline{\mbox{\boldmath$B$}}{})\right){\overline{\mbox{\boldmath$B$}}{}^{2}\over 2\mu_{0}}.$$

$$q_{\mathrm{p}}(\beta )=\frac{q_{p0}}{1+{\beta }^2/{\beta }_{\mathrm{p}}^2}=\frac{{\beta }_{\star }^2}{{\beta }_{\mathrm{p}}^2+{\beta }^2},$$

$$\frac{\lambda }{{\eta }_{\mathrm{t}}{k}^2}\approx 3{\beta }_{\star }\frac{k_{\mathrm{f}}/k}{k{H}_{\rho }}-1.$$

$$\displaystyle\text{DNS:}\quad{{\rm D}{}\mbox{\boldmath$U$}{}\over{\rm D}{}t}$$

$$\displaystyle\text{MFS:}\quad{{\rm D}{}\overline{\mbox{\boldmath$U$}}{}\over{\rm D}{}t}$$

Full text

††volume: 327††journal: IAUS 327: Fine Structure and Dynamics of the Solar Atmosphere††editors: A.C. Editor, B.D. Editor & C.E. Editor, eds.

1 Introduction

Over the past four centuries, the study of sunspots has advanced significantly from simple counts and drawings to detailed monitoring from space. A lot of effort has gone into the compilation of a coherent sunspot record (Clette et al.,, 2014). Despite the long-term record and the high spatial and temporal resolution nowadays available, the theoretical understanding of the processes leading to the formation of sunspots is still limited.

One of the first theoretical descriptions of sunspot formation and evolution was provided by Parker, (1955). He argues that the solar dynamo produces magnetic fields in the form of concentrated flux tubes, which are able to rise through the convection zone to the surface by magnetic buoyancy. An initially horizontal flux tube is expected to break through the surface, forming a typical bipolar sunspot pair. Originally, Parker, (1955) expected the flux tubes to be just $20\,{\rm Mm}$ below the surface, but later he argued that they should be near the overshoot layer (Parker,, 1975).

There has been a range of different approaches to this scenario: from thin flux tube models (e.g. Caligari et al.,, 1995) to detailed simulations of the emergence of horizontal field over the last few megameters (e.g. Rempel and Cheung,, 2014). Caligari et al., (1995) studied the process of the emergence from the tachocline, concluding that a critical magnetic field strength of $\approx 10^{5}\,{\rm G}$ in the subadiabatic layer is needed for the magnetic buoyancy to kick in and raise the tube. They also studied some of the observed properties of active regions, such as the tilt angle and the asymmetry of the inclination angle and the velocity, as well as the latitudinal emergence. On the other hand, Rempel and Cheung, (2014) studied in detail the rise of a pre-injected magnetic flux tube near the surface and the characteristics of bipolar spot formation in the photosphere. Turbulent flows tend to fragment and disperse the injected coherent magnetic flux rope, leading to the decay of the spot. Further examples of theoretical models of flux tube emergence can be found in Cheung and Isobe, (2014).

2 Change of paradigm

This traditional and most studied paradigm to explain the formation of sunspots is part of a more general picture of the Sun. In this framework, which is also referred to as the Babcock-Leighton dynamo (Babcock,, 1961; Leighton,, 1964, 1969), the poloidal magnetic field gets sheared into a toroidal one by the differential rotation of the Sun. This toroidal field resides at the tachocline. There, the layer can be stable enough to allow for a magnetic field of $\sim 10^{5}\,{\rm G}$, which it needs to survive its rise as a coherent tube through the convection zone and to break through the surface, forming spots. Therefore, these flux ropes of magnetic field are global features that feed the dynamo of the Sun, and the sunspots are deeply rooted within the convection zone. One of the problems of the stability of these tubes is convection itself. The turbulence tends to diffuse and destroy the tubes, so they have to be strong enough to survive and rise throughout the entire convective envelope.

We could also imagine a different framework, where the dynamo may still be of Babcock-Leighton type, but it operates in the bulk of the convection zone and causes most of the global features of the Sun. Sunspots form under certain favorable conditions (dynamo strength, stratification, scale separation of the convective cells, etc) as the result of an instability. On top of this, we could also picture a scenario where the turbulence itself is driving a large-scale instability and forms the magnetic field concentrations seen as sunspots (Kleeorin et al.,, 1989, 1990, 1996). We argue that the negative effective magnetic pressure instability (NEMPI) is the instability operating at this level; see Brandenburg et al., (2016) for a recent review. These two scenarios are summarized in Table 1.

3 Negative effective magnetic pressure instability (NEMPI)

3.1 MHD simulations

In order to introduce the mechanism responsible for the formation of magnetic field concentrations, we have to go back to the equations used to study the dynamics of the plasma: the magnetohydrodynamics (MHD) equations. We then use simulations to solve the relevant set of equations. One of the main problems in numerical simulations of the Sun is the inability to use realistic parameters. The solar Reynolds number, Re, related to the amount of turbulence in the system, can span from $10^{10}$ to $10^{15}$, requiring a numerical resolution not yet achievable in super-computers. We, then, have to use a much smaller value of Reynolds number and study how the system scales with increasing values of Re.

Here, we solve the MHD equations using both direct numerical simulations (DNS) and mean-field simulations (MFS). The DNS scheme solves the full set of equations included in the problem, for a given numerical resolution. A drawback of this approach is that it can be difficult to capture the effects of small-scale fluctuations. This is different using MFS, where the problem is solved for the large-scale quantities and the contribution of small-scale quantities are parameterized in terms of large-scale quantities. This allows us to make assumptions for the derivation of the parameterizations and consider only the relevant terms. Therefore, comparing DNS results with known MFS physics will lead to a better understanding of the physical system.

3.2 DNS

The full set of equations we solve within a DNS scheme is:

• •

Continuity equation: \partial\rho/\partial t=-\mbox{\boldmath\nabla}{}\cdot(\rho\mbox{\boldmathU}{}),

• •

Momentum equation: {\rm D}{}\mbox{\boldmathU}{}/{\rm D}{}t=-{1\over\rho}\mbox{\boldmath\nabla}{}P_{\rm gas}+{1\over\rho}\mbox{\boldmathJ}{}\times\mbox{\boldmathB}{}+\mbox{\boldmathf}{}+\mbox{\boldmathg}{}+\mbox{\boldmathF}{}_{\nu},

• •

Induction equation: \partial\mbox{\boldmathB}{}/\partial t=\mbox{\boldmath\nabla}{}\times(\mbox{\boldmathU}{}\times\mbox{\boldmathB}{}-\eta\mu_{0}\mbox{\boldmathJ}{}),

where: {\rm D}{}/{\rm D}{}t\equiv\partial/\partial t+\mbox{\boldmathU}{}\cdot\mbox{\boldmath\nabla}{} is the advective derivative, $P_{\rm gas}$ is the gas pressure, $\nu$ is the kinematic viscosity, $\eta$ is the magnetic diffusivity, \mbox{\boldmathB}{}=\mbox{\boldmathB}{}_{0}+\mbox{\boldmath\nabla}{}\times\mbox{\boldmathA}{} is the magnetic field, \mbox{\boldmathB}{}_{0}=(0,B_{0},0) is a weak imposed uniform field, \mbox{\boldmathJ}{}=\mbox{\boldmath\nabla}{}\times\mbox{\boldmathB}{}/\mu_{0} is the current density, $\mu_{0}$ is the vacuum permeability, \mbox{\boldmathF}{}_{\nu}=\rho^{-1}\mbox{\boldmath\nabla}{}\cdot(2\nu\rho\mbox{\boldmath{\sf S}}{}) is the viscous force, {\sf S}_{ij}={\textstyle{1\over 2}}(\partial_{j}U_{i}+\partial_{i}U_{j})-{\textstyle{1\over 3}}\delta_{ij}\mbox{\boldmath\nabla}{}\cdot\mbox{\boldmathU}{} is the traceless rate-of-strain tensor, $f$ is the forcing function (random, white-in-time, plane, nonpolarized waves with a certain average wavenumber $k_{\rm f}$), \mbox{\boldmathg}{}=(0,0,-g) is the gravitational acceleration.

We measure distance in terms of the density scale height, $H_{\rho}=c_{\rm s}^{2}/g$, magnetic field in terms of the equipartition field strength, $B_{\rm eq0}=\sqrt{\mu_{0}\rho_{0}}\,u_{\rm rms}$, and time in terms of the turbulent-diffusive time, $t_{\eta}=(\eta_{\rm t}k_{1}^{2})^{-1}$ (which is proportional to the eddy turnover time, $\tau_{t0}=(u_{\rm rms}k_{\rm f})^{-1}$, since the turbulent magnetic diffusivity is assumed to the equal to the turbulent viscosity: $\eta_{\rm t}=\nu_{\rm t}=u_{\rm rms}/3k_{\rm f}$).

Brandenburg et al., (2011) solved these equations in a small 3D box, with a scale separation of $k_{\rm f}/k_{1}=15$, strong stratification with a density contrast of $\sim 535$, an initial magnetic field along the horizontal $y$-axis of $B_{0}/B_{\rm eq0}=0.05$, a fluid Reynolds number $\mbox{\rm Re}\equiv u_{\rm rms}/\nu k_{\rm f}$ of 36, a magnetic Prandtl number $\mbox{\rm Pr}_{M}=\nu/\eta$ of 0.5, a magnetic Reynolds number, $\mbox{\rm Re}_{M}=\mbox{\rm Pr}_{M}\mbox{\rm Re}=36$, and detected, for the first time, the spontaneous concentration of magnetic field. This concentration is a consequence of the large-scale instability, i.e., NEMPI, which is described in terms of mean-field equations.

3.3 MFS

Mean-field equations emerge as a result of the decomposition of all quantities as averages (means) and fluctuations, \mbox{\boldmathW}{}=\overline{\mbox{\boldmathW}}{}+\mbox{\boldmathw}{}. Applied to the MHD equations, this decomposition yields:

• •

Continuity equation: \partial\overline{\rho}/\partial t=-\mbox{\boldmath\nabla}{}\cdot(\overline{\rho}\overline{\mbox{\boldmathU}}{}),

• •

Momentum equation: \partial\overline{\mbox{\boldmathU}}{}/\partial t+\overline{\mbox{\boldmathU}}{}\cdot\mbox{\boldmath\nabla}{}\overline{\mbox{\boldmathU}}{}=-{1\over\overline{\rho}}\mbox{\boldmath\nabla}{}\overline{p}+\mbox{\boldmathg}{}+\overline{\mbox{\boldmath{\cal F}}}{}_{\rm M}+\overline{\mbox{\boldmath{\cal F}}}{}_{\rm K},

• •

Induction equation: \partial\overline{\mbox{\boldmathB}}{}/\partial t=\mbox{\boldmath\nabla}{}\times(\overline{\mbox{\boldmathU}}{}\times\overline{\mbox{\boldmathB}}{}+\overline{\mathbf{u}\times\mathbf{b}})+\eta\mbox{\boldmath\nabla}{}^{2}\overline{\mbox{\boldmathB}}{},

where overbars denotes averages, \overline{\mbox{\boldmath{\cal F}}}{}_{\rm M} is the mean Lorentz force and \overline{\mbox{\boldmath{\cal F}}}{}_{\rm K} is the total (turbulent plus microscopic) viscous force, and the turbulent viscous force is determined by Reynolds stresses.

3.3.1 From induction equation: dynamo instability

As a well-known example of mean-field theory, the mean electromotive force \overline{\mbox{\boldmathu}{}\times\mbox{\boldmathb}{}} in the mean field induction equation is caused by turbulent effects which can result in the generation of a large-scale magnetic field (Steenbeck et al.,, 1966; Moffatt,, 1978; Krause and Rädler,, 1980). In the case of isotropic and homogeneous turbulence, the effect of fluctuations can be parameterized as

[TABLE]

Therefore, if $\alpha\neq 0$, there can be a generation of a large-scale magnetic field. This effect is thought to be responsible for generating the Sun’s large-scale magnetic field (Stix,, 1976; Dikpati and Gilman,, 2001; Brandenburg and Subramanian,, 2005).

3.3.2 From momentum equation: NEMPI

For the current discussion, we are not interested in explaining the generation of a magnetic field, but in the study of a concentration of an existing small large-scale magnetic field into the form of surface spots. Therefore, we use the momentum equation:

[TABLE]

The total pressure $P_{\rm Tot}$ in the system is the gas pressure, $P_{\rm gas}$, plus the mean magnetic pressure, \overline{\mbox{\boldmathB}}{}^{2}/2\mu_{0}, and the total (kinetic plus magnetic) turbulent pressure, $p_{\rm turb}=p_{\rm turb}^{(0)}+p_{\rm turb}^{(B)}$, where the total turbulent pressure has a contribution, $p_{\rm turb}^{(B)}$, which depends on the mean magnetic field and one, $p_{\rm turb}^{(0)}$, which is independent on \overline{\mbox{\boldmathB}}{}. We define the effective magnetic pressure {\cal P}_{\rm eff}=\overline{\mbox{\boldmathB}}{}^{2}/2\mu_{0}+p_{\rm turb}^{(B)}, so that the total pressure is:

[TABLE]

where $\tilde{P}_{\rm gas}=P_{\rm gas}+p_{\rm turb}^{(0)}$. The turbulent pressure for isotropic turbulence can be written as:

[TABLE]

Therefore, if the total energy $E_{\rm K}+E_{\rm M}$ is conserved, an increase in the turbulent magnetic energy, $\delta E_{\rm M}>0$, suppresses the turbulent pressure, $\delta p_{\rm turb}=-(1/3)\delta E_{\rm M}<0$, resulting in a negative contribution to the total pressure. For strongly anisotropic turbulence, $\delta p_{\rm turb}=-\delta E_{\rm M}<0$, and the suppression of the turbulent pressure is stronger (Rogachevskii and Kleeorin,, 2007; Brandenburg et al.,, 2016). The contribution to the total turbulent pressure, $p_{\rm turb}^{(B)}$, is parameterized as p_{\rm turb}^{(B)}=-q_{\rm p}(\overline{\mbox{\boldmathB}}{})\,\overline{\mbox{\boldmathB}}{}^{2}/2\mu_{0}, so that the effective magnetic pressure is (Kleeorin et al.,, 1989, 1990, 1993; Kleeorin and Rogachevskii,, 1994; Rogachevskii and Kleeorin,, 2007):

[TABLE]

Substituting the effective magnetic pressure of Equation (5) into Equation (3), we obtain the following expression for the total pressure:

[TABLE]

In cases where the function $q_{\rm p}>1$, the contribution of the effective magnetic pressure to the total pressure is negative, and the gas pressure should increase if the total pressure is balanced sufficiently rapidly with its surroundings. Indeed, a small increase in the magnetic field in the upper layer results in the increase of the gas pressure in this region Therefore, there is a positive gas pressure difference between the upper and lower layer, which drives the downflows. This eventually concentrates the vertical magnetic field within it at the expense of turbulent energy, and the NEMPI is excited.

Effective magnetic pressure computation and parameterization:

The effective magnetic pressure relates the mean magnetic energy to the turbulent one through the function $q_{\rm p}$. We can compute the ad-hoc form of this function by comparing a simulation with an imposed mean magnetic field to a simulation without this field, and extract the effect of the mean magnetic field in the effective magnetic pressure.

Kemel et al., 2012b fitted the DNS results, and approximated the function $q_{\rm p}(\beta)$ by:

[TABLE]

where \beta=\overline{\mbox{\boldmathB}}{}^{2}/B_{\rm eq}^{2}, and $q_{\rm p0}$, $\beta_{\rm p}$, and $\beta_{\star}=\beta_{\rm p}q_{\rm p0}^{1/2}$ are constants. Also, Kemel et al., (2013) established an approximation for the growth rate of the instability, in the case of an isothermal atmosphere and $\beta_{\star}\gg\beta\gg\beta_{\rm p}$:

[TABLE]

Studying the large-scale instability, we can infer the importance of stratification and the scale separation at the onset of the instability. Therefore, we need three main conditions to excite the large-scale instability:

• •

The gradient of effective magnetic pressure must be negative.

• •

Strong enough stratification (small density scale height) $kH_{\rho}$.

• •

Large enough scale separation in the system $k_{\rm f}/k$.

These conditions have already been demonstrated in the DNS of Brandenburg et al., (2012, 2013) and their results set up a starting point for further numerical experiments towards solar-like simulations.

3.4 Some results

To go a step further towards a solar-like model, we focus on two aspects that influence the development of the instability present in the Sun: rotation and a basic representation of a corona envelope.

3.4.1 Effects of rotation

We study the effects of rotation by adding the angular velocity term, -2\mbox{\boldmath\Omega}{}\times\mbox{\boldmathU}{} in the momentum equation, such as \mbox{\boldmath\Omega}{}=\Omega\left(-\sin\theta,0,\cos\theta\right), where $\Omega$ is the rotation rate and $\theta$ is the colatitude of the Sun ($\theta=0$ at the poles). Therefore, we use the following equations of motion in DNS and MFS:

[TABLE]

As explained in Losada et al., (2012, 2013), we set up MFS and DNS of boxes of size $L^{3}$, and impose a weak horizontal uniform magnetic field in the $y$ direction, \mbox{\boldmathB}{}_{0}. The computational domain has periodic boundary conditions in the $x$ and $y$ directions, while we impose stress-free perfect conductor boundary conditions in the $z$ directions.

We solve these equations in a small 3D box, with a scale separation of $k_{\rm f}/k_{1}$, a strong stratification of a density contrast of $\sim 535$, a initial magnetic field along the y-axes of $B_{0}/B_{\rm eq}=0.05$, a fluid Reynolds number $\mbox{\rm Re}\equiv u_{\rm rms}/\nu k_{\rm f}$ of 36, a magnetic Prandtl number $\mbox{\rm Pr}_{M}=\nu/\eta$ of 0.5, and a magnetic Reynolds number, $\mbox{\rm Re}_{M}=\mbox{\rm Pr}_{M}\mbox{\rm Re}=18$.

In the presence of rotation, we observe the formation of magnetic field concentrations, as seen in Figure 1. We quantify the growth rate of the large-scale instability, normalized by the quantity $\lambda_{\ast 0}\equiv\beta_{\star}u_{\rm rms}/H_{\rho}$. In the DNS we can directly compute the value of $\beta_{\star}$, whereas in the MFS it is an input parameter of the simulation. We study two different sets of values in MFS: (i) with $q_{\rm p0}=20$ and $\beta_{\rm p}=0.167$, motivated by the results of Kemel et al., (2013), and (ii) with $q_{\rm p0}=32$ and $\beta_{\rm p}=0.058$, motivated by the DNS fits. As one of the results, Figure 2 shows the decrease of the growth rate of the instability when we increase the Coriolis number in all these three cases. In a solar context, however, where the Coriolis number strongly decreases towards the surface (see Figure 3), the limitation in Coriolis number of the instability constrains the area where the instability can operate close to the surface, where the Coriolis number is small enough. As another result, we can see a disagreement between DNS and MFS results for values of the Coriolis number larger than $\approx 0.1$. Such values of the rotation rate enable the production of kinetic helicity in the system, and therefore an $\alpha$-effect. We see the development of a Beltrami-like magnetic field generated by the large-scale dynamo together with the NEMPI instability. Jabbari et al., (2014) demonstrated that an $\alpha^{2}$ mean-field dynamo is responsible for the onset of dynamo and the presence of the instability for even larger values of the Coriolis number. In the presence of rotation (angular velocity above a certain threshold) the system is unstable with respect to two different instabilities: the large-scale dynamo instability and NEMPI. In general, we expect these instabilities to couple and therefore make the calculation of growth rate much more challenging.

3.4.2 Two-layer stratified model

Our previous models used either vertical field or perfect conductor condition at the top boundary for the magnetic field, when we either impose a vertical or horizontal magnetic field, respectively. These constrain the dynamics of NEMPI. One way to avoid this problem and to reproduce solar-like conditions is to include a non-turbulent region on top of turbulent layer. We call this layer a coronal envelope as it mimic a idealized corona on top of the convection zone of the Sun. Because of the coronal envelope the magnetic field evolves and changes freely at the surface (boundary between turbulent and non-turbulent layer). This approach have been used successfully for dynamo simulation of forced turbulence (Warnecke and Brandenburg,, 2010; Warnecke et al.,, 2011) and turbulent convection Warnecke et al., (2012); Warnecke et al., 2013a ; Warnecke et al., 2016a . The first attempt to combine this two-layer with the study of NEMPI can be found in Warnecke et al., 2013b ; Warnecke et al., 2016b .

Now the initially uniform $y$ component of magnetic field is able to change orientation at the surface, resulting in the spontaneous formation of bipolar regions, which form, evolve and disappear at the surface. We study the effects of slow rotation in a two-layer stratified model, using different resolutions, rotation rates and colatitudes. Figure 4 shows the formation of such bipolar regions in the case of slow rotation at the surface (left panel), and a vertical average (right panel) at the time of the maximum. Now, with the corona envelope, the instability is able to change the orientation of the initial homogeneous $y$ directed imposed field, generating a field in the vertical direction that concentrates in the form of bipolar regions.

A good proxy to understand how the instability operates is the evolution of the magnetic energy spectrum shown in Figure 5. Initially, the magnetic energy is concentrated at the forcing scale $k_{\rm f}/k_{1}=30$. When the bipolar region starts forming, the energy is transferred to a larger scale, with a maximum at $k=2$, when we see the bipolar structure appear at the surface. The instability is then suppressed, and the structures decay, transferring back the energy to smaller scales in the system.

As in previous results, we also see the suppression of the instability as we increase the Coriolis number, but now the presence of the corona layer helps the instability survive even at higher rotation rates.

4 Local helioseismology

Helioseismology is the only technique that allows us to study the solar interior and infer its inner structure and flows. Furthermore, local helioseismology allows the study of small patches of the Sun, opening the possibility to compute flows and magnetic fields. Local helioseismology make use of different methods, including ring-diagram analysis, time-distance helioseismology, helioseismic holography, and the Fourier-Hankel spectral method. For a full review on local helioseismology, see Gizon and Birch, (2005). However, the physical interpretation of these techniques is problematic in the presence of strong magnetic fields. The variation and suppression of acoustic power in sunspots causes anisotropies in the properties of the wave propagation, and therefore uncertainties in the subsurface flows. Although time-distance helioseismology and ring-diagram analysis produce qualitatively similar results (Kosovichev et al.,, 2011), and although some uncertainties can be addressed using different algorithms and approximations, the effects of magnetized subsurface turbulence in the methods is still not well accounted for, so the separation of magnetic and thermal effects remains open (Kosovichev,, 2012). There, most of the approximations are invalid and we still lack a reliable technique proven to reproduce the impact of magnetic fields on the internal traveling waves. Different studies of local helioseismology measure the flows around and beneath sunspots, but their results are contradictory (Jensen et al.,, 2001; Komm et al.,, 2008; Zharkov and Thompson,, 2008; Zhao et al.,, 2010). Studies combining photospheric measurements of the observed magnetic field vectors and subphotospheric flows derived from time-distance helioseismology yield similar velocity patters: inward flows in the sunspot umbra, and outward flows in the surrounding areas. There are also some differences: the flux emergence-related surface flows, like the separate motion of leading and following polarity, the fast rotation, or the apparent shear, do not have their counterpart in the subphotospheric layer. Other studies try to find subsurface signatures of the emergence of sunspots before their appearance, for example Kosovichev et al., (2016) detected signatures of flux emergence at a depth of 62–75 Mm, 12 hours before the first bipolar magnetic structures were seen in the photosphere. Using the same active region, however, different analysis techniques still give different results (Ilonidis et al.,, 2011; Braun,, 2012). Moreover, the attempt to reproduce simulated MHD data from time-distance helioseismology demonstrates the difficulty to recover horizontal flows near active regions, and even the inability to recover vertical flows at any part of the simulation (DeGrave et al.,, 2014).

Another possible approach in trying to study active regions and their emergence characteristics is the analysis of the region just before the emergence of the strong magnetic field, and the statistical analysis of several such emergences. Leka et al., (2013) compiled a list of pre-emerging active regions and analyzed them using helioseismic holography (Birch et al.,, 2013; Barnes et al.,, 2014). They selected around 100 active regions using GONG and SOHO/MDI, and computed averages and statistical emergence properties, such as flows. Statistically, they could not find large flows of more than $15\,{\rm m}\,{\rm s}^{-1}$ in the top 20 Mm below the photosphere on the day prior to the visible emergence of the active region, but could find some signatures of a downflow right before the emergence; see also the recent studies by Birch et al., (2016) and Schunker et al., (2016).

An independent approach is to measure strengthening of the $f$-mode prior to active region formation (Singh et al.,, 2014). It turned out that isolated active regions show a strengthening some 1–2 days prior to magnetic flux emergence (Singh et al.,, 2016). Such a slow build-up of subsurface magnetic flux is incompatible with a buoyant rise and suggestive of an in situ process such as NEMPI.

Yet another idea is to infer the flows in sunspots and active regions by studying their surroundings. This is possible with the Fourier-Hankel spectral method, which allows the separation of waves into inward and outward traveling ones, while using only the region around sunspots. This method was first proposed by Braun et al., (1987), and applied to different active regions, pores and the quiet Sun by a number of authors (Braun et al.,, 1988, 1992; Bogdan et al.,, 1993; Braun,, 1995; Crouch et al.,, 2005; Couvidat,, 2013), and also to the computation of the meridional circulation (Braun and Fan,, 1998; Doerr et al.,, 2010). These authors agree on the absorption of the incoming waves in the sunspots regions, which can reach up to 50% of the $p-$modes power and depends on the depth, frequency, and ridge mode.

We aim to combine these two approaches, by studying pre-emergence active regions using Hankel analysis. Although a comparison between regions with sunspots, pores and quiet-Sun has been made in these different studies, no comprehensive study comprises the pre-emergence phase of an active region, and there is a lack of results for high spherical degree $l$. Therefore, we hope to complete the picture of Hankel analysis on sunspots by selecting a group of active region and studying their emergence properties. This will be subject of a future publication.

5 Conclusions

Even though sunspots and active regions are of vital importance in solar physics and play an important role in many solar phenomena such as for example surface magnetic fields, coronal mass ejections, and solar wind, we lack a complete theoretical understanding to explain their formation and evolution. Does the deeply-rooted flux tube model provide a complete description of sunspot formation? Can NEMPI give a possible explanation for a shallow formation process? Do we have to develop a mixed model where flux tubes emerge and concentrate via NEMPI? We still lack the means to figure out if NEMPI is working in MHD simulations like those of Rempel and Cheung, (2014). Or do other mechanisms play an important role or even explain totally the formation of active regions? There are a number of MHD simulations where magnetic structures of different sizes form spontaneously, like the thermo-hydromagnetic instability proposed by Kitchatinov and Mazur, (2000), or the small-scale vortices at intersections of the intergranular lanes can develop strong vortical downdrafts leading to a spontaneous magnetic flux concentration (Kitiashvili et al.,, 2010), or the pore formation in the Stein and Nordlund, 2012a simulations.

NEMPI may well be responsible for concentrating magnetic fields in the shallow layers of a solar-type simulation. In view of application to the Sun as a next step, it will be important to consider more realistic modeling and include the effects of turbulent convection. It was demonstrated that the relevant NEMPI parameter $q_{\rm p}$ is indeed much larger than unity Kemel et al., 2012a , favoring the possibility of NEMPI. Turbulent convection simulations with an imposed magnetic field (Käpylä et al.,, 2016) yielded structures that are strongly reminiscent of those found in realistic solar surface simulations in the presence of full radiative transport (Stein and Nordlund, 2012b, ). However, in many existing convection simulations, (see e.g. Kitiashvili et al.,, 2010; Käpylä et al.,, 2016; Masada and Sano,, 2016), unlike the case of forced turbulence, the scale separation between the integral scale of the turbulence and the size of the domain is not large enough for the excitation of NEMPI and the formation of sharp magnetic structures. Therefore, the direct detection of negative effective magnetic pressure in turbulent convection with dynamo-generated magnetic fields is a difficult problem. Thus, additional work is needed in order to obtain a complete picture.

Apart from theoretical developments, new and longer high resolution observations are needed to answer some of these questions. In particular, it is crucial to have subsurface maps of downflows/upflows before and during the emergence of the sunspots, or a better description of their structure deeper down, and an improved treatment of the waves around strong magnetic fields.

Acknowledgements.

We thank Emanuel Gafton and Ariane Schad for useful discussions and comments on the manuscript. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Program (FP/2007-2013)/ERC Grant Agreement no. 307117. This work was partly supported 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 (IRL, AB), and the Research Council of Norway under the FRINATEK grant No. 231444 (AB, NK, IR). J.W. acknowledges funding by the Max-Planck/Princeton Center for Plasma Physics and funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement No. 623609. We acknowledge the allocation of 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 National Supercomputer Centers in Linköping and the High Performance Computing Center North in Umeå. This interdisciplinary collaboration would not have been possible without the support of the Solarnet “Mobility of Young Researches” program, awarded to I.R. Losada.