Magnetic concentrations in stratified turbulence: the negative effective magnetic pressure instability

In the presence of strong density stratification, hydromagnetic turbulence attains qualitatively new properties: the formation of magnetic flux concentrations. We review here the theoretical foundations of this mechanism in terms of what is now called the negative effective magnetic pressure instability. We also present direct numerical simulations of forced turbulence in strongly stratified layers and discuss the qualitative and quantitative similarities with corresponding mean-field simulations. Finally, the relevance to sunspot formation is discussed.


Introduction
Magnetohydrodynamic (MHD) turbulence has been studied for a long time, starting with early work on the energy spectrum [1,2] in the 1960s. In many subsequent studies, the effect of gravity was either not considered or it was thought to being just part of the convective driving of the turbulence. The idea that gravity itself could be responsible for causing qualitatively new phenomena in turbulence hardly occurred. This has changed dramatically in just the last few years. An important example is the combined action of gravity with an imposed vertical magnetic field. This drives cross helicity [3], which is an important invariant in ideal MHD. Another example, that will form the main focus of this review, is the spontaneous production of large-scale magnetic flux concentrations in small-scale MHD turbulence by the negative effective magnetic pressure instability (NEMPI), which has its roots in early analytic work [4][5][6][7][8][9], and emerged recently as a pronounced effect in direct numerical simulations (DNS); see Refs. [10][11][12][13][14].
In MHD turbulence there is another important effect that leads to the formation of magnetic structures, namely the dynamo instability. In that case, gravity is unimportant for structure formation, although it does often play a role in driving turbulence, for example through convection. Thus, the dynamo effect must not be confused with the type of structure formation where gravity is a crucial ingredient. Furthermore, convection leads to converging downdrafts that enhance the magnetic field by compression and tend to expel it from diverging flow regions [15,16].
NEMPI, or some similar process, in conjunction with dynamo theory, is one of the contenders in explaining the surface activity of the Sun and other stars. The other main contender is the rising flux tube scenario by which strong coherent flux tubes are being built in the tachocline [17][18][19], which is the shear layer between the convection zone and the radiative interior. However, the helioseismic signatures of such a scenario [20] have not been detected [21]. Observations are more consistent with a gradual build-up of an active region on the timescale of one to two days [22]. An entirely different kinematic process that can form magnetic concentrations is flux expulsion, by which magnetic fields are expelled from regions of rapid motion. A classical example is a convection cell where magnetic field is swept away from the diverging upflows of granules into intergranular lanes and vertices [23,24]. Results from relatively weakly stratified numerical simulations of convection can be explained by this process [16,25,26], but its role in the presence of strong stratification has not yet been studied. Numerical simulations with realistic surface physics have successfully produced active region formation from an unstructured initial magnetic field [27], but it is still a large leap to modeling actual sunspots [28].
Meanwhile, several simulations have displayed spontaneous magnetic structure formation. Some of them involve turbulent convection [16,27,[29][30][31], or a stably stratified polytropic atmosphere [32], so it remains to be clarified, whether gravity plays a direct role, or whether the magnetic field concentrations are mainly the result of converging downdrafts. Other simulations involve forced turbulence in isothermally stratified domains, where no thermally driven convection is possible [10][11][12][13][14][33][34][35], and yet one finds the formation of large-scale magnetic structures. What is remarkable is that these structures extend over the scale of many turbulent eddies. This property suggests that they should be amenable to a mean-field treatment involving averaged, effective equations.
A mean-field approach relevant for describing magnetic effects on the mean flow was developed nearly three decades ago [4][5][6][7][8][9], but only in recent years, with the assistance of DNS, has it gained sufficient attention. In the following, we review the essential properties of NEMPI, discuss analytical approaches to the understanding of the behavior in the presence of either horizontal or vertical magnetic fields, and then turn to DNS whose results can be understood in terms of NEMPI.

Effective magnetic pressure
In the following, we discuss the formation of magnetic structures through a reduction of turbulent pressure by the large-scale magnetic field. For large magnetic Reynolds numbers this suppression of the turbulent pressure can be large enough so that the effective large-scale magnetic pressure (the sum of non-turbulent and turbulent contributions to the large-scale magnetic pressure) becomes negative. The essence of this effect is as follows. The momentum equation describing the plasma motions reads where g is the acceleration due to gravity, is the momentum stress tensor, U and B are the velocity and magnetic fields, p and ρ are the fluid pressure and density, δ ij is the Kronecker tensor, ν is the kinematic viscosity, and is the trace-free rate of strain tensor. We have adopted units where the vacuum permeability µ 0 is set to unity.
Neglecting correlations between velocity and density fluctuations for low-Mach number turbulence, the averaged momentum equation is where ρ is the mean fluid density, U is the mean fluid velocity, Π ij = Π m ij + Π f ij is the mean momentum stress tensor split into contributions resulting entirely from the mean field (indicated by superscript m) and those of the fluctuating field (indicated by superscript f). In the following, the fluctuations of velocity and magnetic field are defined as u = U − U and b = B − B, respectively, where B is the mean magnetic field. The tensor Π m ij has the same form as equation (2), but all quantities have now attained an overbar, i.e.
where p is the mean fluid pressure. The contributions Π f ij , which describe the effect of turbulence on the large-scale Lorentz force, are determined by The turbulent stress tensor, Π f ij , together with the stress tensor describing the mean field contributions, Π m ij , comprise the total mean momentum tensor. Let us first consider isotropic turbulence. The total (hydrodynamic plus magnetic) turbulent pressure p turb (i.e., the isotropic part of Π f ij ) is, in this case, given by [36,37] where E K = 1 2 ρ u 2 is the kinetic energy density of the turbulent (small-scale) motions, E M = b 2 /2 is the energy density of the magnetic fluctuations. The different coefficients in the "turbulent" equation of state, equation (6), for the total turbulent pressure are caused by the fact that the contribution of velocity fluctuations to p turb is determined by the Reynolds stresses, u i u j , which, for isotropic turbulence, are On the other hand, the contribution of magnetic fluctuations to p turb is determined by the Maxwell stresses, In homogeneous turbulence with a uniform large-scale magnetic field, the total turbulent energy density E T = E K + E M is constant, because dissipation is compensated by a continuous supply of energy [8,9], so This implies that a uniform large-scale magnetic field performs no work on the turbulence. It can only redistribute energy between hydrodynamic and magnetic fluctuations. Equation (7) is a steady-state solution for the budget equation for the total turbulent energy density for a time-independent energy source of homogeneous turbulence, I T , with a zero mean velocity and weakly non-uniform mean magnetic field where J = ∇×B is the mean current density, τ 0 is the correlation time of the turbulent velocity field at the integral scale ℓ 0 of turbulent motions, and η t is the turbulent magnetic diffusion. The last term, −E T /τ 0 , in the right-hand side of equation (8) determines the dissipation of the total turbulent energy for the large fluid and magnetic Reynolds numbers. Equations (6) and (7) allow us to determine the change of the total turbulent pressure δp turb in terms of the change of magnetic energy density δE M , Equation (9) implies that the total turbulent pressure is reduced when magnetic fluctuations are generated (δE M > 0). Let us now consider anisotropic turbulence with a preferred direction parallel to some unit vectorê. Specifically, we assume the velocity to be given in the form u = u ⊥ + u zê , where u z = u zê and u ⊥ = u − u z are the velocities parallel and perpendicular toê. We characterize the degree of anisotropy by the parameter σ = u 2 ⊥ /2u 2 z − 1. Thus, for isotropic three-dimensional turbulence we have σ = 0, while for strongly anisotropic turbulence σ is large. For two-dimensional turbulence (an extremely anisotropic case) the parameter σ → ∞. We can then write Substituting equations (10) and (11) into equation (5), the isotropic part of the tensor Π f ij , which yields the total turbulent pressure for anisotropic turbulence, becomes The change of the total turbulent pressure δp turb can then be written in terms of the change of magnetic energy density δE M as For strongly anisotropic turbulence with σ ≫ 1 we have Therefore, for strongly anisotropic turbulence the reduction of the total turbulent pressure by the large-scale magnetic field is six times larger in comparison with that for isotropic turbulence [see equations (9) and (14)]. The turbulent stress tensor, Π f ij , is split into parts that are independent of the mean magnetic field (they determine turbulent viscosity and background turbulent pressure), and parts that do depend on the mean magnetic field. In the presence of a non-zero mean magnetic field only the difference in the stress tensor, depends on the mean magnetic field B, where Π To parameterize the tensor ∆Π f ij we use symmetry arguments, which allow us to construct a symmetric tensor with two preferential directions along the mean magnetic fieldβ = B/|B| and the gravity fieldĝ = g/g. Such a symmetric tensor should be a linear combination of symmetric tensors δ ij ,β iβj andĝ iĝj . These arguments yield [9] where q s , q p , and q g are functions of magnetic Reynolds and Prandtl numbers as well as the modulus of the normalized mean field, β = B/B eq . Here, B = |B| and B eq = (ρ u 2 ) 1/2 is the equipartition field strength. Additional contributions to the tensor ∆Π f ij involving, for example, the mean current density are possible, but will not be considered here.
In summary, the effective mean Lorentz force, which takes into account the turbulence effects, reads Equation (17) has been derived using the spectral τ approach [5,8,9], the renormalization approach [7] and the quasi-linear approach [38,39]. The sum of non-turbulent and turbulent contributions to the large-scale magnetic pressure determines the effective magnetic pressure, When the turbulent contribution, q p (β), becomes large enough, i.e. q p (β) > 1, the effective magnetic pressure is negative, and a large-scale instability, namely NEMPI, can be excited in strongly stratified turbulence. This effect will be studied in the next subsections.

Growth rate of NEMPI for horizontal field
To elucidate the mechanism of NEMPI for a horizontal field, we follow Ref. [13] and consider an equilibrium with zero mean velocity U = 0 and a weak imposed constant mean horizontal magnetic field, B 0 = (0, B 0 , 0). We use the linearized MHD equations neglecting for simplicity the terms which are proportional to the turbulent viscosity and turbulent magnetic diffusivity caused by the electrical conductivity of the plasma. To grasp the essence of NEMPI, we assume that we can apply the anelastic approximation, where ∇ ln ρ = (0, 0, −1/H ρ ) and H ρ is the density scale height. We consider the case H ρ = const and ∇ y = 0 with a weak imposed magnetic field in the y direction. Since the mean magnetic field is independent of y, the mean magnetic tension vanishes in the equation of motion, so where p tot = p + p eff is the sum of mean gas pressure p and effective magnetic pressure p eff . We use the y-component of the linearized induction equation whereŨ andB y are small perturbations of the mean velocity and magnetic fields. To eliminate the gradient of the total pressure in equation (20), we take twice the curl of this equation and linearize, and obtain where P eff (β) = p eff /B 2 eq , and perturbations of the effective magnetic pressure arẽ We introduce a new variable: V z = √ ρŨ z , and use equations (19) and (21), which yield where v A (z) = B 0 / ρ(z) is the mean Alfvén speed. When the characteristic scale of the spatial variation of the perturbations of the magnetic and velocity fields is much smaller than the density scale height, H ρ , the growth rate of the instability is Equation (26) implies that a necessary condition for the large-scale instability is In summary, the mechanism of NEMPI with a horizontal weak imposed magnetic field is as follows. Equation (19) shows that a downward motion U z < 0 leads to a compression: ∇·U < 0. This enhances an applied magnetic field locally [see equation (24)] and results in a large-scale instability, i.e., NEMPI. This instability causes the formation of strongly inhomogeneous magnetic structures.

Growth rate of NEMPI for vertical field
To consider the mechanism of NEMPI for a vertical imposed magnetic field, B 0 = (0, 0, B 0 ), we follow Ref. [40]. The MHD equations for small perturbations are We study an axisymmetric problem in cylindrical coordinates, (r, ϕ, z), apply the anelastic approximation, use the magnetic vector potential A, and introduce the stream function Ψ, i.e., Using the radial components of equations (28) and (29), we obtain for Φ(t, r, z) = ρ −1 ∇ z Ψ the expression where ∆ s is the radial part of the Stokes operator, The density profile in an isothermal layer is ρ = ρ 0 exp(−z/H ρ ). The solution of equation (32) is given by where k r is a suitably defined radial wavenumber and J 1 (x) is the Bessel function of the first kind, which is an eigenfunction of the radial part of ∆ s , i.e., ∆ s J 1 (k r r) = −k 2 r J 1 (k r r). Substituting equation (33) into equation (32), we obtain forΦ(z) the equation For |∇ 2 zΦ (z)|/(k 2 rΦ ) ≪ 1, the growth rate of NEMPI for a vertical imposed field is given by A necessary condition for NEMPI with vertical field is, again, given by equation (27).
In summary, the mechanism of NEMPI for a vertical imposed field is as follows. The downflow removes gas from the upper parts of the turbulent region so that the pressure decreases, which results in a return flow that draws with it more vertical field. This can lead to magnetic field amplification to a strength that exceeds the equipartition field strength in the upper part of the turbulent region [11].

Early simulations toward NEMPI
Numerical simulations of NEMPI began in 2010 [41]. An immediate goal of those early simulations was the demonstration of an instability in a strongly stratified layer using a mean-field parameterization. Another goal was the verification of the foundations of NEMPI using DNS. The main effect that the effective magnetic pressure can be negative, is demonstrated even in the absence of stratification. Thus, DNS in triply periodic domains with an imposed magnetic field have been used. As the imposed magnetic field is increased, the turbulence becomes suppressed (i.e., the turbulent kinetic energy E T decreases), while the energy of the magnetic fluctuations E M increases; see the left panel of figure 1. Note that the total turbulent energy is actually not quite constant as theoretically expected (see Sect. 2.1), but it shows a small dip. This is because for larger fields the turbulent correlation time becomes weakly dependent on the mean magnetic field. However, even if the energy were constant, the applied magnetic field would always have a negative effect on the total (magnetic plus kinetic) turbulent pressure (see Sect. 2.1). The resulting dependence of q p on B 0 /B eq together with the contributions from velocity fluctuations, q K p , and magnetic fluctuations, q M p , is shown in the right panel of figure 1. The measured DNS parameter q s has always turned out to be compatible with zero [29,38,41]. We shall therefore ignore this term in future considerations. The value of q g is usually also found to be negligible. Simulations of convection suggest that q g is positive and of the order of 100, but such a seemingly large values has still only a minor effect in mean-field simulations (MFS) [29]. A horizontal magnetic field B 0 = (0, B y , 0) is imposed by writing the field as B = B 0 + ∇ × A and evolving only the magnetic vector potential A. The appropriate boundary condition in that case is the perfect conductor boundary condition, A x = A y = A z,z = 0. This condition is applied both in the present MFS as well as in the corresponding DNS discussed below.
Having estimated q p , we can now use this term in a mean-field model in which the effective (mean-field) magnetic pressure is parameterized appropriately by replacing If the value of the turbulent magnetic diffusivity is small enough, the mean-field system shows an instability (i.e., NEMPI) [38] which manifests itself mainly in a growth of the mean flow, because the equilibrium mean flow velocity is zero; see left panel of figure 2.
Perturbations to the mean magnetic field also grow exponentially, but they are harder to detect directly because of the much stronger imposed magnetic field. The instability develops concentrations with a certain horizontal wavelength λ ⊥ ≈ 10 H ρ0 or horizontal wavenumber k ⊥ = 2π/λ ⊥ ; see right panel of figure 2. Here, H ρ0 is the density scale height at some reference depth slightly below the surface. Thus, we have k ⊥ H ρ0 ≈ 0.6, which is roughly consistent with subsequent studies [12,13,40].
During the nonlinear evolution of NEMPI, the magnetic structures shown in figure 3 take the form of a droplet and look like a balloon hanging upside down filled with water. In earlier papers we referred to the associated downflows as the "potato sack" effect [10,41]. We discuss this effect in more detail later in Sect. 5.

Detection of NEMPI in DNS
Let us now turn to the explicit verification of NEMPI using DNS in strongly stratified forced turbulence. Early attempts to detect NEMPI failed and it was clear only afterwards why no NEMPI developed: the scale separation ratio, k f /k 1 , was only about five [12]. Here, k f is the forcing wavenumber and k 1 is the smallest wavenumber that fits into the domain. Simulations with k f /k 1 = 15 did finally show NEMPI [10], but the effect was still relatively weak and became clearly noticeable only after having averaged along the direction of the imposed mean magnetic field. Nevertheless, as seen in figure 3, there is a clear sign of the typical droplet shape associated with subsequent downward motion.
Another reason for not having noticed NEMPI in earlier DNS studies could be related to the fact that the field strength was not in the right range. The understanding of this aspect came as a benefit of having used idealizing condition such as an isothermal equation of state. In that case the scale height is constant and independent of z, so the system is similar at all depths except that the density changes. For a given imposed field  strength B 0 , there will always be one particular height where B 0 /B eq takes the preferred value for NEMPI to develop (B 0 /B eq is in the range 0.03-0.2; see table 1 of Ref. [32]). If the field strength is increased, NEMPI develops simply at a larger depth [42]. This is shown in figure 4.
In figure 5 we show the magnetic field for a run with k f /k 1 = 30, so NEMPI is now stronger than before and the flux concentrations can clearly be seen in snapshots even without averaging. Here, time is given in turbulent-diffusive times, τ td = (η t0 k 2 1 ) −1 , where η t0 = u rms /3k f is the estimated turbulent magnetic diffusivity. Again, the field concentration develops first near the top of the domain and then sinks downward. Clearly, this is just opposite to the usual magnetic buoyancy instability [43,44], where magnetic fields rise toward the surface. This underlines the physical reality of a negative effective magnetic pressure.
Making the domain wider results in the replication of statistically similar flux concentrations. In figure 6 we see that the wavelength is approximately 8 H ρ , so k ⊥ H ρ = 2π/8 ≈ 0.8. Similar behavior is also found in MFS, both in two and three dimensions [40].
In Ref. [13], an expression for the growth rate of NEMPI was derived using the anelastic approximation under the assumption H ρ = const (see for details Sect. 2). Remarkable agreement with the numerical calculations has been reported in Ref. [13], where several examples were shown that demonstrated quantitative agreement between DNS and MFS; see also Ref. [40].

NEMPI with vertical field
In the presence of a vertically imposed magnetic field, B 0 = (0, 0, B 0 ), the appropriate boundary condition is the so-called vertical field boundary condition, A x,z = A y,z = A z = 0. It turns out that with a vertical field, the effect of NEMPI is much stronger; see figure 7 for a visualization of B z on the periphery of the computational domain at different times. In this case, B z /B eq can reach and even exceed unity. One reason for this is that the downflows associated with NEMPI lead to a converging return flow in the upper layers, which pinches the field further. Cross sections of the resulting magnetic field look quite different from standard visualizations of buoyant flux tubes that rise and pierce the surface. Here, the magnetic field seems to diffuse out as one Evolution from a uniform initial state toward a circular spot for B z0 /B eq0 = 0.02. Here, B z /B eq0 is shown on the periphery of the domain. Dark shades correspond to strong vertical fields. Time is in units of τ td . An animation is available on http://youtu.be/Um 7Hs 1RzA. Adapted from Ref. [11]. Figure 8. Cuts of B z /B eq (z) in the xy plane at the top boundary (z/H ρ = π) and the xz plane through the middle of the spot at y = 0. In the xz cut, we also show magnetic field lines and flow vectors obtained by numerically averaging in azimuth around the spot axis. Adapted from Ref. [11].
goes deeper down, see figure 8, where the field lines have been computed after taking an axisymmetric average of the magnetic field obtained from the DNS.
The resulting flux concentrations from NEMPI with a vertical imposed magnetic field are much stronger compared to the case with a horizontal field because of the apparent absence of the aforementioned potato sack effect for vertical field. This effect has been observed in turbulence with a horizontal imposed magnetic field; see figures 3 and B eq . Adapted from Ref. [40]. and 4. The potato sack effect is a direct consequence of the negative effective magnetic pressure, making such magnetic structures heavier than their surroundings [10,29,41]. The potato sack effect removes horizontal magnetic field structures from regions in which NEMPI is excited and pushes them downward. For vertical magnetic field, the heavier fluid moves downward along the field without affecting the flux tube, so that NEMPI is not stabilized prematurely by the potato sack effect.
In a mean-field framework, it is quite straightforward to produce an axisymmetric model of a magnetic spot [40]. However, it is important to make sure that the outer radius of the domain is chosen in a suitable manner. If it is too big, downflows will develop on the rim of the cylinder; see figure 8 of Ref. [40]. Furthermore, for an isothermal gas it is straightforward to extend the domain arbitrarily in the vertical direction upward and downward. Figure 9 shows the resulting flux concentration in a domain tall enough so that the field becomes uniform both far above and far below the flux concentration.

Bipolar regions and dynamo-generated magnetic fields
In reality, there will never be a purely vertical nor a purely horizontal imposed field, but a dynamo-generated one which always has some natural horizontal variability. This would allow the field to penetrate the surface. Thus, when allowing for a layer above which there is no turbulence, one can see the development of bipolar regions [45,46]. This layer is meant to represent the effects of a free surface or a photosphere. An example of this is shown in figure 10. Here, a Cartesian domain of isothermally stratified gas was divided into two layers. In the lower layer, turbulence was forced with transverse nonhelical random waves, whereas in the upper layer no flow was driven. A weak uniform magnetic field was imposed in the entire domain at all times. Formation of bipolar magnetic structures was found over a large range of parameters. The magnetic structures became more intense for higher stratification until a density contrast of around 100 across the turbulent layer was reached. The magnetic field in bipolar regions was found to increase with higher imposed field strength until the field became comparable to the equipartition field strength of the turbulence. A weak imposed horizontal field component turned out to be necessary for generating bipolar structures. In the case of bipolar region formation, an exponential growth of the large-scale magnetic field was found, which is a clear indication of a hydromagnetic instability [46]. Additionally, the flux concentrations were correlated with strong large-scale downward-oriented converging flows. These findings strongly suggest that NEMPI is indeed responsible for magnetic flux concentrations found in this system. Many subsequent studies of NEMPI have been undertaken in the meantime. Of particular interest is the case where the main field is not an imposed one, but the result of a dynamo. This situation was first studied in MFS in spherical shells [47], and later also in Cartesian domains [34]. An example is shown in figure 11. Here, the large-scale dynamo is the result of the combined action of stratification and rotation giving rise to kinetic helicity of the turbulence. On the other hand, if rotation is too strong, NEMPI will be suppressed [48,49].
The suppression of NEMPI by rotation came as a surprise, especially because the critical values of the maximum permissible rotation speed were found to be rather low. In dynamo theory the importance of rotation on the flow is usually measured in terms of the Coriolis number, Co = 2Ωτ td , but those values are only around 0.1 when NEMPI begins to be suppressed; see figure 2 of [49]. They argued that the reason for this is the   fact that the growth time for NEMPI is longer than the turnover time of the turbulence τ td . If one normalizes instead with the typical growth rate of NEMPI, λ 0 * the critical values of 2Ω/λ 0 * are found to be slightly above unity.
Although NEMPI does not appear to be excited when rotation is too strong, simulations show that magnetic flux concentrations are still being produced when there is strong stratification and a dynamo is operating preferentially in the deeper parts [34]. In figure 12 we show a case where a sub-equipartition strength magnetic field is produced in the deeper parts and leads to a super-equipartition strength magnetic fields in the surface layers. The different polarities can then be driven together to form sharp structures in the form of an inverted Y-shaped pattern in a vertical cross-section [35,50]. This is an example where the phenomenon of turbulent magnetic reconnection has been seen to occur in a natural setting. Similar behavior has also been seen in global spherical shell dynamos [51]. In simulations [35,50], the generated magnetic field reached superequipartition levels so rapidly that it was not possible to detect NEMPI during the growth of the magnetic field. It should be noted that NEMPI cannot be excited for super-equipartition magnetic fields.
To investigate the role of magnetic reconnection, the flow around the sharp interface was zoomed in, and the dynamics of the current sheet in this region was studied [35]. The reconnection rate was determined independently through the inflow velocity in the vicinity of the current sheet and via the electric field in the reconnection region. For large Lundquist numbers (Lu > 1000), the reconnection rate was found to be nearly independent of the value of Lu [35], where Lu = v A L c /η with L c being the length of the current sheet. This agrees with earlier studies of turbulent reconnection [52], which also showed independence of Ohmic resistivity [53,54], as well as results of recent numerical simulations performed by other groups in simpler settings [52,[55][56][57][58].

Turbulent convection
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 convectively-driven turbulence. It is not obvious that NEMPI applies straightforwardly to this case. It was shown already early on that-even in convection-the relevant NEMPI parameter q p is indeed much larger than unity [29], favoring the possibility of NEMPI. Subsequent convection simulations with an imposed magnetic field [30] yielded structures that are strongly reminiscent of those found in realistic solar surface simulations in the presence of full radiative transport [27].
The spontaneous formation of surface magnetic structures from a large-scale α 2 dynamo by strongly stratified thermal convection in Cartesian geometry has recently also been studied by [31]. They found that large-scale magnetic structures are formed at the surface only in cases with strong stratification. The presence of rapid uniform rotation was an argument in [31] that NEMPI seems not to be responsible for these structures. On the other hand, the combined effect of rapid uniform rotation and stratification can produce helicity and an α effect, which causes a large-scale α 2 dynamo. A similar situation was encountered in connection with the dynamo simulations of forced turbulence, where rotation also did not suppress the formation of structures [35]. Thus, the question of the origin of these structures remains unsettled.
The direct detection of negative effective magnetic pressure in turbulent convection with dynamo-generated magnetic fields is a difficult problem. In many existing convection simulations, 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.

Cross helicity effect
The role of NEMPI is not always evident, especially when the magnetic field strongly exceeds the equipartition value. However, all these simulations have in common that there is a vertical magnetic field such that g · B = 0. Interestingly, g · B is not only a pseudoscalar, but it is odd in the magnetic field. In MHD, there is an important invariant of the ideal equations, namely the cross helicity u · B [59]. It is often not important unless it was present in the initial conditions. However, it is not too surprising that u · B is generated whenever g · B is nonzero. This was explored in Ref, [3], where it was found that If the turbulence intensity is nonuniform, there is yet another contribution to u · b that is proportional to η t (∇ ln u rms ) · B, which was first obtained in Ref. [60]. All the simulations that produce large-scale magnetic field structures display what can be characterized as inverse cascade or inverse transfer behavior of magnetic field from the driving scale to large scales with a horizontal wavenumber k ⊥ such that k ⊥ H ρ ≈ 0.8 or even less. The role of the conservation property of u · B still needs to be explored, but it is clear that there is a remarkable analogy between inverse transfer seen in figure 13 and that of a large-scale dynamo, where the conservation of magnetic helicity is known to lead to an inverse cascade [61][62][63].

Concluding remarks
The purpose of this review has been to highlight the fact that the presence of gravitational stratification introduces a qualitatively new phenomenon in MHD turbulence, namely the formation of large-scale magnetic structures via excitation of NEMPI. Turbulence causes a modification of the large-scale magnetic pressure, so that the effective magnetic pressure becomes negative for large fluid and magnetic Reynolds numbers, and this results in the excitation of NEMPI. DNS demonstrate that the effective magnetic pressure can be negative in any kind of turbulence, e.g., in nonstratified and stratified isothermal turbulence, polytropic stably stratified turbulence, turbulent convection, and in turbulence with an upper non-turbulent layer. However, the actual instability is excited only in stratified turbulence when the initial mean magnetic field is less than the equipartition field. For very large fluid and magnetic Reynolds numbers, NEMPI weakly dependent on the level of turbulence. In some cases, where there is locally a vertical magnetic field, NEMPI causes the formation of magnetic spots.
In astrophysics, there are many examples where MHD turbulence is accompanied by strong density stratification. A prediction from our studies reviewed in this paper would be that such systems should exhibit magnetic spots. We do not know whether there is a direct relation to sunspots, which are generally hypothesized as being the result of deeply rooted thin magnetic flux tubes. Whether such isolated flux tubes really exist and how they are formed is an open question.
We do know of magnetic flux tubes in MHD turbulence [64,65], that are analogous to vortex tubes in hydrodynamic turbulence [66], but these tend to scale with the resistive scale [67], so such tubes would become smaller as the magnetic Reynolds number is increased. Global simulations of convective spherical shell dynamos have been used to visualize magnetic flux tubes [68,69]. Those may well be the type of tubes seen in earlier Cartesian simulations, but they could also be local field enhancements resulting from the large-scale dynamo. It is hard to imagine that these flux structures alone can explain the formation of sunspots, unless there was some kind of reamplification. Clearly, as future simulations of global dynamos gain in resolution, they would eventually display spots, just like the Sun and other stars do. It will then be important to have possible frameworks in place for understanding such spots. We hope that the present review has provided some relevant inspiration beyond the standard paradigm.
Further steps toward more realistic modeling of the formation of magnetic spots and bipolar regions include replacing forced turbulence by self-consistently driven convective motions that are influenced by the radiative cooling at the surface together with partial ionization. Including more realistic physical processes at the solar surface might also help to reproduce the surrounding spot structures.