1 Introduction

The solar convection zone is the ultimate driver of activity in the solar chromosphere and corona. It is the only available source of mechanical energy. Upper atmosphere activity and heating is empirically intimately connected with the presence of magnetic fields. Hence the need to understand the behavior of magnetic fields at the solar surface. The solar magnetic field is produced by dynamo action within the convection zone. Thus, to understand the energy source of the chromosphere and corona we need to understand the solar dynamo, magneto-convection, and the transport of magnetic flux through the convection zone. For recent reviews see Fan (2009) and Charbonneau (2010). For probing the subsurface layers of the Sun, our best tools are the various techniques of local helioseismology (Gizon and Birch, 2005). Accurate modeling of the rise through the convection zone and emergence of magnetic flux, of sunspots and active regions is needed for improving helioseismic probing of solar subsurface structure.

Convection is the transport of energy by bulk mass motions. In a convection zone, energy is transported as thermal energy, except in layers where hydrogen is only partially ionized, in which case most of the energy is transported as ionization energy. Typically, the motions are slow compared to the sound speed so that approximate horizontal pressure balance is maintained. As a result, warmer fluid is less dense and buoyant while cooler fluid is denser and gets pulled down by gravity. For a detailed review of solar surface convection see Nordlund et al. (2009).

The topology of convection is controlled by mass conservation (Stein and Nordlund, 1989). Convection has a horizontal cellular pattern, with the warm fluid ascending in separate fountainlike cells surrounded by lanes of cool descending fluid. In a stratified atmosphere, with density decreasing outward, most of the ascending fluid must turn over and be entrained in the downflows within a density scale height (ignoring gradients in velocity and filling factor). Fluid moving a distance Δr in an atmosphere with a density gradient dln ρ/dr would, if its density remained constant, be overdense compared to its surroundings by a factor Δρ/ρ = -(dlnρ/drr. This is unstable and produces a pressure excess in the upflow cell interiors that pushes the fluid to turn over into the surrounding downflow lanes. Since the fluid velocity decreases inward from the top of the convection zone, its derivative has opposite sign to that of the density, so the length scale for entrainment is increased. Warm upflows diverge and tend to be laminar, while cool downflows converge and tend to be turbulent. Temperature in stellar convection zones increases inward, so the scale height and, as a result, the size of the horizontal convective cellular pattern also increase inward. Think of the rising fluid as a cylinder. As described above, most of the fluid entering at the bottom of the cylinder must leave through its sides within a scale height. If the ratio of vertical to horizontal velocities does not change much with depth, then the radius of the cylinder can increase in proportion to the scale height and still maintain mass conservation (Stein and Nordlund, 1998).

The solar surface is covered with magnetic features with spatial sizes ranging from unobservably small to hundreds of megameters. Their distribution is featureless (Parnell et al., 2009; Thornton and Parnell, 2011). Large-scale magnetic structures, sunspots and active regions, possess some well defined global properties (Hathaway, 2010). The main observed properties of small scale magnetic structures are (de Wijn et al., 2009): Strong fields tend to be vertical and weaker fields horizontal. The strongest vertical fields are in pressure equilibrium with their surroundings and tend to occur in the magnetic network and the intergranular lanes. Horizontal fields are found predominantly inside granules and near the edges of granules. Horizontal field properties are similar in the quiet Sun, plage, and polar regions (Ishikawa and Tsuneta, 2009). Three orders of magnitude more magnetic flux emerges in the quiet Sun than emerges in active regions (Thornton and Parnell, 2011). This new flux is first seen as horizontal field inside granules followed by the appearance of vertical field at the ends of the horizontal field (Centeno et al., 2007; Martínez González and Bellot Rubio, 2009; Ishikawa et al., 2010; Guglielmino et al., 2012).

In the presence of magnetic fields, convection is altered by the Lorentz force, while convection influences the magnetic field via the curl (v × B) term in the induction equation.

Where the conductivity high, the magnetic field is frozen into the ionized plasma. Where the magnetic field is weak (magnetic energy small compared to kinetic energy), convective motions drag it around. To maintain force balance, locations of higher field strength (higher magnetic pressure) tend to have smaller plasma density and lower gas pressure. Diverging, overturning motions quickly sweep the field (on granular times of minutes) from the granules into the intergranular lanes (Tao et al., 1998a; Emonet and Cattaneo, 2001; Weiss et al., 2002; Stein and Nordlund, 2004; Vögler et al., 2005; Stein and Nordlund, 2006). In hours (mesogranular times), the field tends to collect on a mesogranule scale. In days (supergranule times), the slower, large scale supergranule motions collect the field in the magnetic network at the supergranule boundaries. Convective flows produce a hierarchy of loop structures in rising magnetic flux. Slow upflows and buoyancy raise the flux, while fast downflows pin it down, which produces Ω- and U-loops (Cheung et al., 2007). The different scales of convective motion produce loops on these different scales, with smaller loops riding piggy-back in a serpentine fashion on the larger loops (Cheung et al., 2007; Stein et al., 2010b). Dynamo action occurs in the turbulent downflows where the magnetic field lines are stretched, twisted, and reconnected, increasing the field strength (Nordlund et al., 1992; Cattaneo, 1999; Vögler and Schüssler, 2007; Schüssler and Vögler, 2008; Pietarila Graham et al., 2010).

Magnetic fields influence convection via the Lorentz force, which inhibits motion perpendicular to the field. As a result, the overturning motions that are essential for convection are suppressed and convective energy transport from the interior to the surface is reduced. Radiative energy loss to space continues, so regions of strong field cool relative to their surroundings. Being cooler, these locations have a smaller scale height. Plasma drains out of the magnetic field concentrations in a process called “convective intensification” or “convective collapse” (Parker, 1978; Spruit, 1979; Unno and Ando, 1979; Nordlund, 1986; Bercik et al., 1998; Grossmann-Doerth et al., 1998; Bushby et al., 2008). This process can continue until the magnetic pressure (plus a small gas pressure) inside the flux concentration equals the gas pressure outside, giving rise to a field strength much greater than the equipartition value with the dynamic pressure of the convective motions. These magnetic flux concentrations are cooler than their surroundings at the same geometric layer. However, because they are evacuated, their opacity is reduced so photons escape from deeper in the atmosphere (Wilson depression, Maltby, 2000). Where the magnetic concentrations are narrow, there is heating from their hotter side walls and they appear as bright points (Spruit, 1976). Where the concentrations are wide, the side wall heating is not significant and the flux concentrations appear darker than the surroundings as pores or sunspots.

Magnetic fields alter granules’ properties - producing smaller, elongated, lower intensity contrast, “abnormal” granules (Muller et al., 1989; Title et al., 1992; Bercik et al., 1998; Vögler, 2005; Cheung et al., 2007). Strong magnetic flux concentrations typically form in convective downflow lanes, especially at the vertices of several such lanes, due to the sweeping of flux by the diverging convective upflows (Vögler et al., 2005; Stein and Nordlund, 2006). They are surrounded by downflows which sometimes become supersonic.

Magneto-convection simulations have been very useful in understanding and interpreting observations. Sánchez Almeida et al. (2003), Khomenko et al. (2005), Shelyag et al. (2007), and Bello González et al. (2009) have used simulations to calibrate the procedures for analyzing and interpreting Stokes spectra in order to determine the solar vector magnetic field. Fabbian et al. (2010) has shown that magnetic fields alter line equivalent widths by altering the temperature stratification and by Zeeman broadening. These two effects act in opposite directions, but still leave a net result and hence alter abundance determinations. Zhao et al. (2007), Braun et al. (2007), Kitiashvili et al. (2009), Birch et al. (2010), and Braun et al. (2012) have used convection and magneto-convection simulation results to analyze local helioseismic inversion methods.

Magneto-convection is highly non-linear and non-local, so it needs to be modeled using numerical simulations. Two complementary approaches are being used to study magneto-convection, which we will call “idealized” and “realistic”. “Idealized” studies ignore complex physics by assuming a fully ionized, ideal plasma and energy transport by thermal conduction. Magneto-convection in the deep, slow moving, adiabatic portion of the convection zone satisfies these idealized assumptions and, in addition, can use the anelastic approximation whereby acoustic waves are eliminated from the calculation, which permits larger time steps. “Idealized” calculations are important for isolating and studying fundamental physical phenomena as well as for exploring parameter space (because they run fast). “Realistic” studies include complex physics - an equation of state for partially ionized gas, non-grey radiation transport and, in some cases, even some non-equilibrium effects. “Realistic” calculations are necessary to make quantitative comparisons with observations in order to understand the observations and to provide artificial data for evaluating data analysis procedures. In this review we focus on the “realistic” numerical modeling of solar surface magneto-convection. It updates and extends the section on magneto-convection from the review by Nordlund et al. (2009) of solar surface convection.

It is organized as follows: Section 2 states the equations that need to be solved. Section 3 describes the solar observations. Section 4 describes the simulation results for: dynamo action (4.1), flux emergence (4.2), flux concentrations (4.3), and pores and sunspots (4.4).

2 Equations

To simulate magneto-convection, the conservation equations for mass, momentum, energy and the induction equation for the magnetic field must be solved, together with Ohm’s law for the electric field and an equation of state relating pressure to the density and energy. For a detailed discussion of the equations governing convection see Nordlund et al. (2009).

Mass conservation controls the topology of stratified convection,

$$\frac{{\partial p}} {{\partial t}} = - \nabla \cdot (\rho u),$$
((1))

where ρ is the density and u the velocity.

Momentum conservation controls the plasma motions. In the presence of magnetic fields, convection is altered by the Lorentz force in the momentum equation, which becomes:

$$\frac{{\partial (\rho u)}} {{\partial t}} = - \nabla \cdot (\rho uu), - \nabla P - \rho g + J \times B - 2\rho \Omega \times u - \nabla \cdot \tau _{visc} .$$
((2))

Here P is the pressure, g is the gravitational acceleration, B is the magnetic field, J = ∇ × B/μ is the current, μ is the permeability (magnetic constant), and τvisc is the viscous stress tensor,

$$\tau _{ij} = \mu \left( {\frac{{\partial u_i }} {{\partial x_j }} + \frac{{\partial u_j }} {{\partial x_i }} - \frac{2} {3}\nabla \cdot u\delta _{ij} } \right).$$
((3))

The Lorentz force inhibits motion perpendicular to the field. As a result, the overturning motions that are essential for convection are suppressed and convective energy transport from the interior to the surface is reduced. When large depths are included where the fluid motions become slow, the coriolis force, -2ρΩ× u, needs to be included. Angular momentum conservation then produces a surface shear layer with the surface rotating slower than the bottom of the domain.

Kinetic energy is changed by energy transport and work against the forces acting on the plasma. The equation for kinetic energy is

$$\frac{\partial } {{\partial t}}\left( {\frac{1} {2}\rho u^2 } \right) = - \nabla \cdot (\frac{1} {2}\rho u^2 u) - u \cdot \nabla P + \rho u \cdot g + u \cdot J \times B + u \cdot \nabla \cdot \tau _{visc} .$$
((4))

Note that if there is no net mass flux, 〈ρuxyt = 0, then there is no net buoyancy work by gravity. The vertical convective velocity and density are correlated. Upflows have lower density and cover a larger area while downflows have higher density and cover a smaller area. Downflows are pulled down by gravity and upflows are buoyant. Gravity drives the convection, doing positive work on both the upflows and downflows. But the total work by gravity vanishes. Hence, the positive work on the convective motions is balanced by an equal but negative work on the mean flow. There is necessarily a horizontally averaged mean flow in the opposite direction to gravity. Such mean flows do exist in the simulations.

Internal energy is changed by transport, by PdV work, by Joule heating, by viscous dissipation, and by radiative heating and cooling. It is the fluid version of the 2nd law of thermodynamics and (together with the density) determines the plasma temperature, pressure, and entropy.

$$\frac{{\partial e}} {{\partial t}} = - \nabla \cdot (eu) - P(\nabla \cdot u) + Q_{rad} + Q_{visc} + \eta J^2 ,$$
((5))

where e is the internal energy per unit volume. The radiative heating/cooling is:

$$Q_{rad} = \int_\nu {\int_\Omega {\rho \kappa } } _\nu (I_\nu - S_\nu )d\Omega d\nu .$$
((6))

Here κv is the opacity (1/ρκv = lv is the mean free path of photons of frequency v), \(I_\nu (r,\hat n,t)\) is the radiation intensity (energy at frequency v, at location r, travelling in direction \(\hat n\) , at time t, per unit area, per unit solid angle, per unit frequency, per unit time), and Sv = v/κv is the source function (v is the rate of energy emission at frequency v per unit frequency, per unit mass, per unit time, per unit solid angle). The viscous dissipation is:

$$\begin{array}{*{20}c} {Q_{visc} = \sum\limits_{ij} {t_{ij} \frac{{\partial u_i }} {{\partial x_j }}} } \\ { = \frac{\mu } {2}\sum\limits_{ij} {\left[ {\frac{{\partial u_i }} {{\partial x_j }} + \frac{{\partial u_j }} {{\partial x_i }} - \frac{2} {3}\delta _{ij} \sum\limits_k {\frac{{\partial u_k }} {{\partial x_k }}} } \right]^2 .} } \\ \end{array}$$
((7))

Magnetic energy changes due to transport by the Poynting flux, E × B/μ, work by the Lorentz force, u · J × B and joule dissipation, J · E,

$$\frac{\partial } {{\partial t}}\left( {\frac{{B^2 }} {{2\mu }}} \right) = - \nabla \cdot [E \times B/\mu ] - u \cdot J \times B - J \cdot E.$$
((8))

Adding kinetic, internal and magnetic energy equations gives the equation for the total energy, ET = 1/2ρu2 + e + B2/2μ,

$$\frac{\partial } {{\partial t}}\left( {\frac{1} {2}\rho u^2 + e + B^2 /2\mu } \right) = - \nabla \cdot \left[ {\left( {\frac{1} {2}\rho u^2 + e + P + B^2 /2\mu } \right)u + E \times B/\mu + u \cdot \tau _{visc} } \right] + \rho u \cdot g + Q_{rad} .$$
((9))

Convection influences the magnetic field via the curl (u × B) term in the induction equation,

$$\frac{{\partial B}} {{\partial t}} = - \nabla \times E,$$
((10))

where the electric field is given by Ohm’s Law. In a one-fluid MHD system, it is

$$E = - u \times B + \eta J + \frac{1} {{en_e }}(J \times B - \nabla P_e ),$$
((11))

where η is the resistivity, ne is the electron number density, Pe is the electron pressure, and e is the electron charge. The last two (Hall and pressure) terms are usually neglected, but the Hall term may be important in the weakly ionized photosphere. Where the magnetic field is weak and the resistivity low, the field is frozen into the ionized plasma. Convective motions drag the field around, stretching and twisting it.

To make “realistic” models of solar surface magneto-convection it is necessary to include all the significant physical processes occurring near the solar surface. In the photosphere and upper convection zone Local Thermodynamic Equilibrium (LTE) is a good approximation. For models that extend into the chromosphere non-local thermodynamic equilibrium (NLTE) effects must also be included, as is being done in the bifrost code (Gudiksen et al., 2011).

Ionization energy accounts for 2/3 of the energy transported near the solar surface and must be included to obtain the observed solar velocities and temperature fluctuations (Stein and Nordlund, 1998). An equation of state (EOS) is used to determine the pressure and temperature for the partially ionized plasma. Typically this is in tabular form and includes LTE ionization of the abundant elements as well as hydrogen molecule formation as a function of log(density) and internal energy per unit mass.

Radiation from the solar surface cools (and heats) the plasma and produces the low entropy, high density fluid whose buoyancy work drives the convective motions. Since the optical depth is of order unity in these regions, neither the diffusion approximation (Ustyugov, 2010) nor the escape probability approach (Abbett and Fisher, 2010) are sufficiently accurate. Radiative heat/cooling is calculated by explicitly solving the radiation transfer equation in both continua and lines,

$$\frac{{\partial I_\nu }} {{d\ell }} = \varepsilon _\nu - \chi _\nu I_\nu ,$$
((12))

where \(I_\nu (r,t;\hat n)\) is the intensity at frequency v, position r, time t, in direction \(\hat n\), v is the radiation emission at that frequency, location, time, and in that direction, and χv is the absorption probability of radiation at that frequency, location, time, and direction. Since time steps in the solar case are short, the plasma state is known from the previous time step, so the emission and absorption can be calculated and the radiation transfer can be solved explicitly. The main problem is the large number of frequencies necessary to accurately represent the radiation. A multi-group approximation is used to drastically reduce the number of frequencies for which the transfer equation is solved. The opacity and emissivity are grouped into a few bins (typically 4 – 12) according the magnitude of the opacity and its frequency (Nordlund, 1982; Skartlien, 2000; Stein and Nordlund, 2003; Vögler et al., 2004; Vögler, 2004). The number of directions is also limited.

3 Observations

The solar surface is covered with magnetic features with spatial scales from smaller than can currently be resolved (∼ 70 km with the Swedish 1 m Solar Telescope) to active regions covering up to 100 Mm (Figure 1). These evolve on a correspondingly wide range of time scales, from seconds for the smallest observed features, to months for some active regions. If one counts as a single feature any contiguous collection of the same polarity with magnitude above some cutoff, then the magnetic flux distribution is a power law of slope -1.85 (Figure 2). Alternatively, if one identifies features as individual flux peaks, then the distribution is log-normal (Parnell et al., 2009; Thornton and Parnell, 2011). These power laws are featureless, they have no peaks or valleys.

Figure 1
figure 1

Stokes V in blue wing of 630.2 nm line from Hinode. Around a sunspot (a) and in the quiet Sun (b), showing the wide range in size of magnetic structures on the Sun. The dimensions of both figures are 110 Mm and the pixel size is 108 km. Image reproduced by permission from Parnell et al. (2009), copyright by AAS.

Figure 2
figure 2

Histogram of feature frequency vs. magnetic flux. Image reproduced by permission from Parnell et al. (2009), copyright by AAS.

The large-scale magnetic structures, sunspots, and active regions possess some well defined global properties: all active regions in a given hemisphere have the same polarities of leading/following spots, but reversed between the northern and southern hemispheres (Hale’s polarity law); the polarities reverse in a semi-periodic 22 year cycle; in each cycle spots first appear at mid latitudes and then their appearances migrate toward the equator; active regions are tilted with the leading spot closer to the equator (Joy’s law); trailing spot fields migrates toward the poles and sunspots tend to reappear at certain active longitudes (Figure 3). See review by Hathaway (2010). These properties imply the existence of a global dynamo process.

Figure 3
figure 3

Magnetic butterfly diagram from longitudinally averaged radial magnetic field. This illustrates Hale’s polarity law, Joy’s law, transport of flux toward the poles and migration of active region emergence sites toward the equator. Image reproduced by permission from Hathaway (2010).

There is an excellent review of small-scale solar magnetic field observations (network and internetwork quiet Sun) by de Wijn et al. (2009). The main observational results are: strong fields tend to be vertical and weaker fields horizontal. Vertical kilogauss fields (in pressure equilibrium with their surroundings) are found in the magnetic network and as isolated, intermittent concentrations in intergranular lanes. Horizontal magnetic fields are found all over the Sun (Trujillo Bueno et al., 2004; Harvey et al., 2007), predominantly inside and near the edges of granules. They are transient, intermittent and have granule-scale sizes and lifetimes and strengths in the hectogauss range (generally less than equipartition with the convective dynamic pressure) (Ishikawa and Tsuneta, 2010). Weaker horizontal fields have no preferred orientation. Stronger ones tend to align with the active regions. The horizontal field properties are similar in the quiet Sun, plage, and polar regions (Ishikawa and Tsuneta, 2009). The spatially averaged horizontal magnetic field strength is 50 – 60 G, while the spatially averaged vertical field strength is only 11 G (Lites et al., 2008). This may be due to the larger area covered by horizontal fields compared to the isolated vertical field concentrations. There is no characteristic size or lifetime for the horizontal fields (they have an exponential distribution both in size and lifetime) (Danilovic et al., 2010). There is some question of the accuracy of Hinode determinations of quiet Sun transverse magnetic fields due to s/n problems (Borrero and Kobel, 2011).

Bright points in the G-band (dominated by CH molecular transitions) have been used as proxies for the magnetic field. At disk center small magnetic concentrations appear as bright points in the intergranular lanes, while larger concentrations are dark. The increased brightness in magnetic concentrations is due to their lower density compared with their surroundings. At a given geometric height, granules are hotter than the intergranular lanes, which are, in turn, hotter than G-band bright points. Although at a given geometric height the magnetic elements are cooler than the surrounding medium, one sees into deeper layers, to where the temperature is higher, due to the reduced opacity and to heating from the hot surrounding granules. Locations of large magnetic concentrations are cooler than even the G-band bright points because both convective heat transport and side wall heating are reduced.

The presence of strong magnetic fields enhances the pillow appearance of granules because their low density and resulting low opacity allow one to see deeper into the hot granules behind them (the “hot wall” effect Spruit, 1976, 1977; Keller et al., 2004; Shelyag et al., 2004; Carlsson et al., 2004). Where the fields are strong, the intergranular lanes are depressed up to 350 km below the mean height. Thus the τ = 1 surface is extremely corrugated. Toward the limb, where the surface is viewed at an angle, the low density and opacity in the strong magnetic elements allows one to see the hot granule walls behind. These are the faculae (Figure 4) (Keller et al., 2004; Carlsson et al., 2004). The excess brightness comes from a thin layer (∼ 30 km) of steep density gradient at the interface between the magnetic and nonmagnetic atmospheres. Typically there is a dark lane just centerward of the bright faculae. As the line of sight moves limbward from granule to faculae, it first intersects a granule top and is bright, then intersects cool material above the granule and inside the magnetic concentration, and finally intersects the hot granule wall on the far side of the magnetic concentration. Variations in the field strength produces variations in the density and opacity which leads to a striated appearance in the bright granule walls. Where the field is weaker, the density is higher, so the opacity larger. This effect is enhanced by a higher CH concentration also due to the higher density. Thus, where the magnetic field is weaker, the radiation emerges from higher, cooler layers, so these locations appear darker.

Figure 4
figure 4

Comparison of G-band intensity at viewing angle μ = 0.63 of observations (left) and at μ = 0.6 simulated (right).

High resolution observations of solar faculae show that they have an asymmetric contrast profile, with some brightness extending up to one arcsecond in the limbward direction from their peak in brightness (Hirzberger and Wiehr, 2005). The wide contrast profile cannot be explained solely by the “hot wall” effect, as was noted by Lites et al. (2004). The works by Keller et al. (2004) and Steiner (2005) address this issues, with somewhat conflicting but broadly consistent explanations. One conclusion is that the limbward extension of brightness comes from seeing the granule behind the facular magnetic field through the rarefied facular magnetic flux concentration; a circumstance that observers suspected decades ago. The explanation is corroborated by the direct comparisons of observations and simulations by De Pontieu et al. (2006).

Three orders of magnitude more magnetic flux is observed to emerge as small scale loops in the quiet Sun than emerges in active regions (Thornton and Parnell, 2011). This new flux is first seen as horizontal field (linear polarization in Stokes spectra) inside granules followed by the appearance of vertical field at the ends of the horizontal field (circular polarized Stokes spectra) (Centeno et al., 2007; Martínez González and Bellot Rubio, 2009; Ishikawa et al., 2010; Guglielmino et al., 2012) (Figure 5).

Figure 5
figure 5

Emergence of a small magnetic loop in the quiet solar photosphere. Top: continuum intensity at 630 nm. Center: total linear polarization (Stokes Q,U) in the 630.25 line. Bottom: total circular polarization (Stokes V). Signal is clipped at ± 0.1 pm. Red contours are linear polarization > 0.22 pm, while black and white contours are circular polarization > ± 0.1 pm. Distances are in acrsec. Image reproduced by permission from Martínez González and Bellot Rubio (2009), copyright by AAS.

These Ω-loop footpoints get quickly swept into the intergranular lanes and the horizontal field to the edges of the granules. They do not show a helical structure. Transient horizontal fields also appear briefly where new downflow lanes form (Danilovic et al., 2010). The flux in these emerging bipoles is small, 1016 × few × 1017 Mx, but their rate of appearance is large, around a few × 10-10 km2 s, hence their dominant contribution to the emerging flux of the Sun (Martínez González and Bellot Rubio, 2009; Ishikawa and Tsuneta, 2009; Ishikawa et al., 2010). Most of these small loops are low lying, with only about a quarter reaching up to chromospheric heights.

4 Simulations

Two types of numerical studies of magneto-convection are being undertaken: “idealized” and “realistic”. Both approaches give valuable, but different, insights into the properties of magnetoconvection. Idealized simulations were pioneered by Weiss (1966) and extensively used by Tao et al. (1998b), Cattaneo (1999), Abbett et al. (2000), Hurlburt and Rucklidge (2000), Emonet and Cattaneo (2001), Weiss et al. (2002), Cattaneo et al. (2003), and Bushby et al. (2008). See reviews by Weiss (1991) and Schüssler (2001). They are especially useful for gaining physical insights into convective properties. In these calculations an ideal gas equation of state is assumed and energy transport is assumed to be only by convection and thermal conduction. For modeling magneto-convection in the solar interior anelastic or reduced sound speed calculations with an ideal gas equation of state and diffusive radiation transport are appropriate (Miesch, 2005; Miesch et al., 2008; Fan, 2009; Miesch et al., 2011; Hotta et al., 2012). An alternative approach applicable to the deep convective layers is to reduce the sound speed (Hotta et al., 2012). This, as in the anelastice approximation, allows larger time steps. “Realistic” simulations were pioneered by Nordlund (1982) and have been extensively developed by Stein and Nordlund (1998, 2006), Steiner et al. (1998); Vögler et al. (2005), Schaffenberger et al. (2005), Hansteen et al. (2007), Abbett (2007), Jacoutot et al. (2008), Carlsson et al. (2010), and Gudiksen et al. (2011). A tabular equation of state is used, which includes the partial ionization of hydrogen, helium and other abundant elements, because below 40,000 K in the Sun ionization energy dominates over thermal energy in convective energy transport. The radiation transfer equation is solved to determine the radiative heating and cooling, because the optical depth is of order unity near the visible solar surface, so that neither the diffusion nor optically thin approximations are valid. Such detailed physics is necessary to make quantitative comparisons with observations. Here we restrict ourselves to the more “realistic” surface simulations. Magneto-convection and dynamo action in the deeper layers of the convection zone are reviewed by Miesch (2005) and Miesch and Toomre (2009).

4.1 Turbulent convection and dynamo action

Meneguzzi et al. (1981) and Cattaneo (1999) were the first to demonstrate, via magneto-convection simulations, that dynamo action will occur in turbulent convection even in the absence of rotation. These calculations were for closed, Boussinesq systems. Questions were raised whether local dynamo action is possible in the highly stratified solar convection zone (Stein et al., 2003) because in a stratified atmosphere with much stronger downflows than upflows, magnetic flux is pumped down (Tobias et al., 2001; Dorch and Nordlund, 2001). Abbett (2007), Vögler and Schüssler (2007), and Pietarila Graham et al. (2010) showed that a local surface dynamo is indeed possible. Vögler and Schüssler (2007) and Pietarila Graham et al. (2010) used a shallow, very high resolution, magneto-convection simulation with no Poynting flux in or out of the domain, but with a high magnetic diffusivity in the bottom boundary layer to mimic the loss of magnetic flux to the deeper convection zone. Pietarila Graham et al. (2010) demonstrated that during the kinematic (linear) growth phase, the primary dynamo process (95%) was stretching of magnetic field lines against the magnetic tension component of the Lorentz force by convective motions at subgranule scales (0.1 – 1 Mm) in the turbulent downdrafts, which generates still smaller scale (20 – 200 km) magnetic field. The other 5% was from work against magnetic pressure by fluid motions at granule scales (Figures 6 and 7). In addition, magnetic pressure also produces a cascade of magnetic energy from the dynamo generated scales to still smaller scales. In the saturated phase, generation by stretching is almost balanced by the compressive cascade with a little energy also going into MHD wave generation.

Figure 6
figure 6

Dynamo energy transfers in the kinematic phase. The dominant process is turbulent stretching of magnetic field lines against the magnetic tension component of the Lorentz force. Image reproduced by permission from Pietarila Graham et al. (2010), copyright by AAS.

Figure 7
figure 7

Dynamo net energy transfer rates as a function of horizontal spatial scale in the kinematic dynamo phase. Left: work against magnetic tension (pink dashed line) and work against magnetic pressure (green dotted line) as a function of the fluid motion spatial frequency. Right: dynamo stretching (blue dot-dash line), dynamo compression (black solid line), and magnetic energy removed by compression (red dotted line) as a function of the magnetic field spatial frequency. Image reproduced by permission from Pietarila Graham et al. (2010), copyright by AAS.

Abbett (2007) and Schüssler and Vögler (2008) showed that such small scale dynamo action produces many low-lying loops with large amounts of horizontal field overlying the granules (Figures 8 and 9). Steiner (2010) argues that the preponderance of horizontal over vertical field is an inherent consequence of the fact that granules are wider than a scale height. Consider an area of length L and height h. The horizontal (ΦH) and vertical (ΦV) fluxes for a loop are the same, so that ΦH = 〈BHLh = ΦV = 〈BVL2, where BH is the horizontal and BV is the vertical field and L and h are the horizontal and vertical extents of the field. Hence 〈BH〉/〈BV〉 ≈ L/h and low lying loops connecting opposite sides of granules must have larger average horizontal than vertical field.

Figure 8
figure 8

Photospheric magnetic field lines showing many low-lying, horizontally directed magnetic structure from a simulation from the upper convection zone to the corona. Image reproduced by permission from Abbett (2007), copyright by AAS.

Figure 9
figure 9

Image of Log horizontal B in vertical slice through saturated phase of dynamo simulation of Vögler and Schüssler (2007). Mean optical depth unity is at approximately 0.9 Mm. Note many loops of different sizes closing in the photosphere. Image reproduced by permission from Schüssler and Vögler (2008), copyright by ESO.

Global dynamos have been simulated by Brun et al. (2004), Miesch (2005), Dobler et al. (2006), Browning et al. (2006), and Brown et al. (2007, 2010). See reviews by Brandenburg and Dobler (2002), Ossendrijver (2003), Miesch and Toomre (2009), and Charbonneau (2010). Note that the fact that both the rate of magnetic flux emergence and the probability distribution of magnetic flux magnitudes are featureless power laws from 1016 – 1023 Mx suggests that the solar dynamo has no preferred scale, that it acts throughout the convection zone with each scale of convective motions generating new flux on that scale (Parnell et al., 2009; Thornton and Parnell, 2011). That is, all the surface magnetic features are produced by a common process (which can not be all dominated by surface effects since the sunspots and active regions clearly are not).

4.2 Subsurface rise and emergence of magnetic flux

Magnetic fields are produced by dynamo action throughout the solar convection zone. Their emergence through the visible surface is driven by two processes: advection by convective upflows and buoyancy (to maintain approximate pressure equilibrium with their surroundings the density inside the concentration is smaller than in its surroundings). Fan (2009) has reviewed the rise of magnetic flux through the deep convection zone. Simulations of magnetic flux emerging from the surface layers of the convection zone have been initiated in three ways: either from coherent twisted flux tubes forced into the computational domain through the bottom boundary or made buoyant by lowering their density (Yelles Chaouche et al., 2005; Cheung et al., 2007, 2008; Martínez-Sykora et al., 2008, 2009; Cheung et al., 2010), or by inflows at the bottom advecting minimally structured, uniform, untwisted, horizontal field advected by inflows into the domain through the bottom boundary (Stein et al., 2010a,b), or produced locally by dynamo action (Abbett, 2007; Abbett and Fisher, 2010). These very different approaches, using several different computer codes, show several, robust, common features.

Convective flows produce a hierarchy of loop structures in rising magnetic flux (Figure 10). Magnetic flux rises through the convection zone because it is advected by broad upflows and because it is buoyant. Along the way, it encounters convective downflows piercing the upflows on smaller and smaller scales, with downflow speeds significantly larger than the upflow speeds. The portions of the magnetic concentration in the downflows will be dragged down, or at least have their upward motion slowed, while the portions still in the upflows or that have large density deficits and so large buoyancy continue to ascend rapidly (Figures 10, 11). The different scales of motions produce a hierarchy of magnetic Ω- and U-loops with small loops riding piggy-back on larger loops in a serpentine structure (Cheung et al., 2007, 2008; Kitiashvili et al., 2010; Stein et al., 2010b) (Figure 12). As a result, emergence of large, undulated Ω-loops occurs as a collection of small-scale, mixed polarity, emergence events. In general, the asymmetry of upflow and downflows (amplitudes and topology) leads to a tendency for downward transport of magnetic flux; a process known as “turbulent pumping” (Drobyshevskii et al., 1980; Nordlund et al., 1992; Petrovay and Szakály, 1993; Tobias et al., 1998; Dorch and Nordlund, 2000, 2001; Tobias et al., 2001) (Figure 10).

Figure 10
figure 10

mpg-Movie (37927878 KB) Still from a movie showing Log B showing multiple loops and several vertical flux concentrations, one of which has become a pore (see Section 4.4). The movie shows that most of the magnetic features are being pushed down by convective downflows, but some of the loops are rising toward the surface and in places loops have opened out through the top boundary leaving vertical “flux tubes” behind. Movie produced by Sandstrom, CSC, NASA Ames Res. Ctr. (For video see appendix)

Figure 11
figure 11

mpg-Movie (194745724 KB) Still from a movie showing Rise and emergence of initially uniform, untwisted horizontal magnetic field continuously being advected into the computational domain by inflows in the centers of supergranule cells at the bottom. Left: |B| image and velocity vectors. Right: |B| image and magnetic field vectors both in a vertical plane. The boundary field strength was gradually increased from 0.2 kG with an e-folding time of 5 hrs until it reached 5 kG and thereafter held constant. (For video see appendix)

Figure 12
figure 12

Several magnetic field lines showing large scale loops with smaller serpentine loops riding piggy-back on them. Shading shows downflows.

_

As the magnetic flux rises it expands (Figure 13). For horizontal flux tubes, the horizontal expansion is much larger than the vertical expansion. Consider a purely horizontal field \(B = B\hat x\) in the x-direction. Suppose the rates of expansion in the horizontal and vertical directions are α = ∂vx/∂x = ∂yy/∂y and ∂vz/∂z respectively. The rate of change of the magnetic field following the fluid motion is given by the Lagrangian derivative

$$\frac{{DB}} {{Dt}} = - (\nabla \cdot u)B + (B \cdot \nabla )u,$$
((13))

which becomes,

$$\frac{{D\ln B}} {{Dt}} = - (1 + \varepsilon )\alpha .$$
((14))

The continuity equation (1) becomes

$$\frac{{D\ln \rho }} {{Dt}} = - \nabla \cdot u = - (2 + \varepsilon )\alpha .$$
((15))

So for emerging horizontal fields, Bρ(1+)/(2+) (Cheung et al., 2010). For isotropic expansion ( = 1), Bρ2/3. For vertical expansion small compared to horizontal expansion ( ≪ 1) Bρ1/2.

Figure 13
figure 13

Time sequence of vertical cross-sections perpendicular to an initial coherent twisted flux tube. The grey scale is temperature and the color coding is magnetic field strength |B|. The purple line is the τ500 = 1. Image reproduced by permission from Cheung et al. (2007), copyright by ESO.

The fields first appear at the surface in localized regions as small bipoles with a small-scale, mixed pepper and salt polarity. The emergence region spreads in time (Figure 14). As the bipoles begin to emerge, horizontal and vertical fields have similar strengths, but horizontal fields are more common (cover more area) than vertical fields, except for the strongest fields. The mixed polarity fields collect into separated unipolar regions due to the underlying large scale magnetic structures (Figure 14).

Figure 14
figure 14

mpg-Movie (169942682 KB) Still from a movie showing Time sequence of emergent continuum intensity (left), vertical magnetic field at τRoss = 0.1 (right), for an initial twisted flux half torus driven into the computational domain at 7.5 Mm depth. The initial central field strength was 21 kG and total flux 7.6 × 1021 Mx. The domain is 92 × 49 Mm wide. Image reproduced by permission from Cheung et al. (2010), copyright by AAS. (For video see appendix)

Diverging, overturning convective motions quickly sweep magnetic fields (on granular times of minutes) from the granules into the intergranular lanes (Figure 15) (Weiss, 1966; Hurlburt and Toomre, 1988; Tao et al., 1998a; Emonet and Cattaneo, 2001;Weiss et al., 2002; Steiner et al., 1998; Stein and Nordlund, 2004; Vögler et al., 2005; Stein and Nordlund, 2006). In hours (mesogranular times) the field tends to collect on a mesogranule scale. Observations averaged over several hours reveal this magnetic pattern (de Wijn et al., 2005; Ishikawa and Tsuneta, 2010), even though there is no distinct convective meso-granule scale. In days (supergranule times) the slower, large scale supergranule motions sweep the fields to the supergranule boundaries. Eventually a balance is reached where the rate of emergence of new flux balances the rate at which flux is swept to larger horizontal scales. This balance empirically occurs at supergranulation scales and produces the magnetic network at the supergranule boundaries. Here the new flux encounters existing magnetic flux, which it either cancels or augments (Simon et al., 2001; Krijger and Roudier, 2003; Isobe et al., 2008).

Figure 15
figure 15

Sweeping of magnetic field into intergranular lanes. Initially the entire surface is covered with 30 G horizontal field. Surface magnetic field strength is shown at 0.5 min (left), 10 min (center) and 30 min (right). Black contours are zero velocity contours which outline the granules. Fields stronger than 0.5 kG appear as white and black. Field magnitudes less than 30 G are shown in grey. Diverging upflows first sweep the granules clear of strong fields and on a longer time scale sweep the interiors of mesogranules free of strong fields. Image reproduced by permission from Stein and Nordlund (2006), copyright by AAS.

The rising magnetic loops are not coherent, but rather have a filamentary structure (Figure 10). Some individual filaments rise more rapidly than others. The small-scale crenulation of the loops produces the “pepper and salt” pattern as the flux emerges through the visible surface. As the bulk of the magnetic loop reaches the surface, the different polarities concentrate in unipolar regions accompanied by flux cancellation where opposite polarities come in contact (Figures 14 and 17). This happens because the mixed polarity emergence is due to the undulating magnetic field lines produced by convective upflows and downdrafts distorting the large loops rising from below. The underlying large-scale magnetic structures organize the mixed polarity fields when they approach the surface. In order for the like polarity branches to collect, the mass trapped in the U-loops between the peaks of the small Ω-loops must be removed. This occurs by the U-loop getting pinched off and forming plasmoids which may either sink below the surface or get ejected into space (Figure 16) (Lites, 2009; Cheung et al., 2010).

Figure 16
figure 16

Schematic scenario of how granular motions pinch off U-tubes to produce O-loops in 2D and plasmoides in 3D. Image reproduced by permission from Cheung et al. (2010), copyright by AAS.

Figure 17
figure 17

mpg-Movie (108657227 KB) Still from a movie showing Emergent continuum intensity (left) and vertical magnetic field at τcont = 0.01 (right) from simulation with initial/boundary condition of convective inflows advecting 1 kG uniform, untwisted, horizontal field into the computational domain at 20 Mm depth. The intensity range is I/〈I〉 = [0.13, 2.5] and the magnetic field range is ± 3.5 kG. The pores may form spontaneously in vertical flux tubes from magnetic loops that have reached the surface and opened out through top boundary. Compare this with Figure 14 for the rise of a coherent twisted flux tube. (Movie shows the initial “pepper and salt” emergence, the horizontal advection of the field, its concentration into unipolar regions with cancellation where opposite polarities meet and merging of like polarities to form pores. Resolution was increased from 48 km to 24 km horizontally at time 51.7 hrs.) (For video see appendix)

Magnetic flux emergence simulations starting with horizontal, uniform, untwisted field at 20 Mm depth is shown in Figures 10, 11 and 18. Figure 18 shows magnetic field lines in the simulation box viewed from the side and slightly above. The red line in the lower left is horizontal field being advected into the domain. In the lower center is a loop like flux concentration rising toward the surface. In the upper right is a vertical flux concentration or “flux tube” through the surface (Stein and Nordlund, 2006). While the field lines form a coherent bundle near the surface, below the surface they become tangled and spread out in many different directions (Vögler et al., 2005; Schaffenberger et al., 2005). A “flux tube” at the surface forms by a loop-like flux concentration rising up through the surface and opening up through the upper boundary where the condition is that the field tends toward a potential field. This leaves behind the legs of the loop. Typically one leg is more compact and coherent than the other and persists for a longer time as a coherent entity while the other is quickly dispersed by the convective motions. Cattaneo et al. (2006) have studied the existence of flux tubes using an idealized simulation of a stably stratified atmosphere with shear in both the vertical and one horizontal direction driven by a forcing term in the momentum equation. They find that in the absence of symmetries, even in this laminar flow case, there are no flux surfaces separating the inside of a flux concentration from the outside, so that the magnetic field lines in the concentration connect chaotically to the outside and the “flux tube” is leaky (Figure 19). The fact that magnetic fields that are concentrated close to the surface tend to tangle and spread out in many directions below the surface has been demonstrated earlier by Grossmann-Doerth et al. (1998) — see also Vögler et al. (2005), Schaffenberger et al. (2005), and Stein and Nordlund (2006).

Figure 18
figure 18

Magnetic field lines in a simulation snapshot viewed from an angle. The red line in the lower left is horizontal field being advected into the domain. In the lower center is a loop like flux concentration rising toward the surface. In the upper right is a vertical flux concentration or “flux tube” through the surface with its field lines connecting chaotically to the outside below the surface. Image reproduced by permission from Stein and Nordlund (2006), copyright by AAS.

Figure 19
figure 19

Magnetic flux concentration at the solar surface and magnetic field lines showing the complex field line connections below the surface. The “flux tube” is a local surface phenomena. Image reproduced by permission from Schaffenberger et al. (2005).

Magnetic fields alter the granule properties - producing smaller, lower intensity contrast, “abnormal” granules (Bercik et al., 1998; Vögler, 2005). Strong magnetic flux concentrations typically form in convective downflow lanes, especially at the vertices of several such lanes, due to the sweeping of flux by the diverging convective upflows (Vögler et al., 2005; Stein and Nordlund, 2006). They are surrounded by downflows which sometimes become supersonic. The normal convective downflows are enhanced surrounding the flux concentrations by baroclinic driving due to the influx of radiation into the concentration (Deinzer et al., 1984; Knölker et al., 1991; Bercik, 2002; Vögler et al., 2005).

Granules become larger and darker as the field first emerges (due to the suppression of vertical motions by the horizontal section of the bipoles and adiabatic cooling due to their expansion) and elongate in the direction of the horizontal component of the field (Figure 20).

Figure 20
figure 20

Emergent continuum intensity as a twisted flux tube emerges through the solar surface. Image reproduced by permission from Yelles Chaouche et al. (2005).

The main differences between these two approaches are that a coherent initial flux tube leads to a more coherent symmetrical structure when it emerges through the surface and field line connections below the surface are more localized. In the minimally structured approach organized magnetic field concentrations develop spontaneously when sufficient flux is present, instead of being imposed as initial and boundary conditions, so the emergent structures are less coherent.

4.3 Small scale flux concentrations

Magnetic concentrations arise from loops emerging into the upper solar atmosphere and leaving their legs behind and from the diverging convective upflows which sweep magnetic field into the intergranular lanes and concentrate the field into sheets and at the vertices of the lanes into “tubes” of magnetic flux (Vögler et al., 2005; Schaffenberger et al., 2005; Stein and Nordlund, 2006). To maintain force balance, locations of higher field strength (higher magnetic pressure) tend to have smaller plasma density and lower gas pressure. Strong magnetic fields, through the Lorentz force, inhibit overturning convective motions and hence the transport of energy toward the surface. Radiative energy loss to space continues, so regions of strong field cool relative to their surroundings at the same geometric layer. Being cooler, these locations have a smaller scale height. Plasma drains out of the magnetic field concentrations in a process called “convective intensification” or “convective collapse” (Parker, 1978; Spruit, 1979; Unno and Ando, 1979; Nordlund, 1986; Bercik et al., 1998; Grossmann-Doerth et al., 1998; Bushby et al., 2008). This process continues until the magnetic pressure (plus a small gas pressure) inside the flux concentration approximately equals the gas pressure outside, giving rise to a field strength much greater than the equipartition value with the dynamic pressure of the convective motions.

The opacity of magnetic flux concentrations is reduced, because they are evacuated, so photons escape from deeper in the atmosphere, that is, optical depth surfaces are depressed into the interior (Wilson depression, Maltby, 2000). Radiatively, they are holes in the surface. The temperature structure in these concentrations is nearly in radiative equilibrium with radiative heating from fluid flowing down along their sides and cooling from emission in vertical rays (Bercik, 2002, Figure 21). Where the flux concentration is narrow, heating from the side walls raises the internal temperature at optical depth unity and the concentration appears bright (Spruit, 1976). Small magnetic flux concentrations may appear especially bright in the continuum (Bercik, 2002; Keller et al., 2004; Carlsson et al., 2004; Steiner, 2010). This enhanced brightness extends for all the sight lines that pass through the low density, optically thinner, magnetic concentration where photons escape from deeper, hotter layers (Figures 22 and 23). Where the concentrations are wide, the side wall heating is not significant and the flux concentrations appear darker than the surroundings as pores or sunspots.

Figure 21
figure 21

Radiative heating and cooling (1010 erg/g/s) in a vertical slice through a magnetic flux concentration. The top two panels show the net heating (yellow & red)/cooling (green & blue) with superimposed contours of temperature (top) and magnetic field (next to top). The bottom three panels show the net heating/cooling for vertical (cos θray,vertical = μ = 1), slanted (μ = 0.5), and nearly horizontal rays (μ = 0.05). Image reproduced by permission from Bercik (2002).

Figure 22
figure 22

Temperature, density, and magnetic field strength along a vertical slice through magnetic and non-magnetic regions, with the average formation height for the G-band intensity for a vertical ray (black line) and at μ = 0.6 (white line). Axes are distances in Mm. The bottom panel shows temperature as function of lg τ500. Image reproduced by permission from Carlsson et al. (2004), copyright by AAS.

Figure 23
figure 23

Schematic sketch of a magnetic flux concentration (region between the thin lines) and adjacent granules (thick lines). The dashed lines enclose the region where 80% of the continuum radiation is formed. Bright facular radiation originates from a thin layer at the hot granule wall behind the limbward side of the optically thin magnetic flux concentration. The line of sight for the dark centerward bands is shown by the dark shaded region. Image reproduced by permission from Keller et al. (2004), copyright by AAS.

In the G-band there is an additional, smaller, effect — the CH molecule becomes dissociated in the low density magnetic concentrations (Steiner et al., 2001; Carlsson et al., 2004; Shelyag et al., 2004; Steiner, 2005). The bottom panel of Figure 22 shows the temperature as function of lg τ500. The contrast in temperature between magnetic concentrations and non-magnetic areas increases with decreasing optical depth giving larger intensity contrast with increasing opacity (e.g., Ca H,K). The G-band has its mean formation height (black line in bottom panel) at lg τ500 = -0.48 corresponding to a mean formation height 54 km above where τ500 = 1, therefore giving a larger contrast than in the continuum. The contrast enhancement by the destruction of CH is seen as a dip in the curve showing the mean formation optical depth in the bottom panel. Note also that the G-band intensity has its peak contribution at similar heights as the continuum (that is why the granulation pattern looks similar). Bright points in the G-band have been used as a proxy for magnetic field concentrations. While G-band bright points are a good proxy for strong magnetic fields, there are many more regions of strong field that appear dark in the G-band, typically because they cover a larger area (Figure 24). In simulations, all the bright points correspond to locations of large field magnitude, but not all large field locations correspond to bright points (Vögler et al., 2005; Stein and Nordlund, 2006). Further, the field has a longer lifetime than the bright points. The contrast in the G-band has also been studied by Rutten et al. (2001), Sánchez Almeida et al. (2001), and Steiner et al. (2001).

Figure 24
figure 24

G-band brightness vs. magnetic field strength at continuum optical depth unity for a snapshot of magneto-convection with a unipolar magnetic field. Note that while all bright points correspond to strong magnetic fields there are many locations of strong field that appear dark in the G-band.

4.4 Pores and sunspots

Recently, “realistic” radiative-convective MHD models of pores and sunspots have become possible. Bercik et al. (2003), Stein et al. (2003), and Vögler et al. (2005) found micropores forming spontaneously in magneto-convection simulations. Cameron et al. (2007a) modeled pores in magneto-convection simulations by imposing them as initial and boundary conditions. Stein et al. (2010b) found that large pores formed spontaneously in deep magneto-convection simulations. Schüssler and Vögler (2006) simulated a sunspot umbra where convection produced umbral dots. Heinemann et al. (2007), Scharmer et al. (2008), and Rempel et al. (2009b) started with flaring, rectangular, slab magnetic concentrations. Rempel et al. (2009a) and Rempel (2011) modeled sunspots in a magneto-convection simulation starting from a pair of axisymmetric, self-similar magnetic funnels. Cheung et al. (2010) modeled sunspots formed from an emerging, twisted half torus magnetic “flux tube”. An excellent review, especially of helioseismic applications, is Moradi et al. (2010). For more details on sunspots see the review by Rempel and Schlichenmaier (2011).

In magneto-convection simulations with initial vertical fields, micropores form spontaneously in vertices of the intergranular lanes where several lanes meet (Bercik et al., 2003; Stein et al., 2003; Vögler et al., 2005). In the typical formation scenario a small bright granule is surrounded by strong magnetic fields in the intergranular lanes. The upward velocity in the small granule reverses and it disappears with the area it occupied becoming dark. The surrounding strong fields move into the dark micropore area (Figure 28).

Figure 28
figure 28

Micropore formation sequence. Left panels: images of the magnetic field strength, center panels: emergent intensity, and right panels: mask showing low intensity, strong field locations. Image reproduced by permission from Bercik (2002) and Bercik et al. (2003).

As the upflow velocity in a flux concentration slows and reverses, the upward heat flux decreases and the plasma inside the concentration cools by radiation through the surface (Figure 25). As a result, the density scale height decreases and the plasma settles lower. Initially the material piles up below the surface until a new hydrostatic structure is established (Figure 25). The micropores are also heated by radiation from their hotter sidewalls (Figure 21, Spruit, 1976, 1977). On the order of a granular timescale the magnetic field is dispersed and the micropore disappears.

Figure 25
figure 25

Magnetic field (filled contours at 250 G intervals from 0 G to 3500 G) and temperature (top) (1000 K intervals from 4000 K to 16,000 K), and ln ρ (bottom) (in 0.5 intervals from -2 to 4). The τ = 1 depth is shown as the thick line around z = 0 Mm. The flux concentrations are significantly cooler than their surroundings and so have a smaller scale height. The established, strong “flux tube” in the center has been evacuated and is in equilibrium. The smaller flux concentrations on either side are in the process of being evacuated, starting above the surface and piling up plasma below the surface. Image reproduced by permission from Bercik (2002); Bercik et al. (2003).

Figure 26
figure 26

Image of vertical velocity (red/yellow down, blue/green up in km s-1) with magnetic field contours at 0.5 kG intervals at the surface (left) and 1.5 Mm below the surface (right). Image reproduced by permission from Bercik (2002).

Pores have developed spontaneously in magneto-convection emerging flux simulations when rising Ω-loops emerge through the surface and the upper boundary, leaving behind vertical magnetic field concentrations (Stein and Nordlund, 2006). The pores grow by accumulating flux from their surroundings. The pore pictured in Figure 27 has reached a flux of 2.4 × 1020 Mx and occupies an area of 6 Mm2. The flux concentration develops first near the surface. It cools and quickly becomes partially evacuated and flux concentration extends downward, reaching all the way to the bottom of the domain (at 20 Mm depth) — Figure 29 and Figure 30, see also Kitiashvili et al. (2010). Most magnetic field lines in the pore connect to the end of a large scale loop rising from the bottom of the domain, although some connect to various other structures. Additional flux is being transported into the pore by horizontal flows along the intergranular lanes. These flows feeding the pore extend to depths of several megameters. The simulated pores have sometimes lasted for a long time — greater than 8 hrs (Kitiashvili et al., 2010) and 12 hrs in our case.

Figure 27
figure 27

Spontaneously formed simulated pore. Clockwise from upper left: emergent intensity, vertical magnetic field at 〈τ〉 = 1, horizontal magnetic field at the same level, inclination angle to the vertical. The pore is edge brightened in part of its circumference. The field is vertical in the pore interior and becomes inclined more than 45° at the pore edge.

Figure 29
figure 29

mpg-Movie (24956644 KB) Still from a movie showing Time sequence of temperature and density fluctuations during pore formation viewed from below (surface is at the bottom). Note the cooler temperatures extending downward from surface, followed by lower densities. Movie produced by Sandstrom, CSC, NASA Ames Res. Ctr. (For video see appendix)

Figure 30
figure 30

mpg-Movie (34719413 KB) Still from a movie showing Time sequence of log |B| during pore formation. The flux concentration forms first at the surface and then extends downward. Near the surface the pore’s field is filamentary, but at large depths it becomes mostly coherent. Movie produced by Sandstrom, CSC, NASA Ames Res. Ctr. (For video see appendix)

Pores, like micropores, are surrounded by downflows concentrated into a few downdrafts. The ubiquitous occurrence of downflows in the close vicinity but outside magnetic flux concentrations (see, for example, also Steiner et al., 1998) has been explained in terms of baroclinic flows by Deinzer et al. (1984). The effect has been observationally verified by Langangen et al. (2007). Pores are edge brightened (Figure 27). Cameron et al. (2007b) explain this as due to the fact that the surface of unit optical depth occurs at slightly higher temperature at the edges of pores, possibly due to decreased overlying density because of the spreading magnetic field.

No one has yet produced a sunspot ab initio. Several “realistic” magneto-convection simulations of sunspots have been made starting with idealized, imposed initial magnetic field configurations. See review by Rempel and Schlichenmaier (2011). Cheung et al. (2010) has come closest, starting from an emerging, twisted, half torus magnetic “flux tube”. Others have started with monolithic, self-similar magnetic configurations. Schüssler and Vögler (2006) investigated magneto-convection in a uniform, vertical field representing a sunspot umbra. Heinemann et al. (2007), Scharmer et al. (2008), and Rempel et al. (2009b) studied flaring rectangular slab field configurations. Rempel et al. (2009a) and Rempel (2011) started with a pair of axisymmetric, self-similar funnels (Schüssler and Rempel, 2005), with the same flux but slightly different field strengths (Figure 31). All these simulations develop thin upflow plumes with surrounding downflows that are the observed umbral dots. The most challenging property of spots to model has been their penumbra, which are found to depend crucially on the existence of very inclined magnetic fields in the outer parts of the spots.

Figure 31
figure 31

Still from a movie showing Emergent intensity and |B| (kG) in opposite polarity spot pair initiated from a pair of axial symmetric, self-similar flaring magnetic field funnels. Each spot has the same flux, but the one on the left has a slightly weaker field. The simulation domain is 98 × 49 Mm by 6 Mm depth. The vertical dimension has been stretched by a factor of 2 in the bottom panel (from Rempel et al., 2009a). (To watch the movie, please go to the online version of this review article at http://www.livingreviews.org/lrsp-2012-4.)

In the Cheung et al. (2010) simulation, spots form from an emerging Ω-loop (Figure 14). The field first emerges with mixed polarities. The opposite polarities then counterstream to collect into the opposite polarity sunspots. This counterstreaming of opposite polarities is due to the underlying large-scale structure of the coherent subsurface roots of the emerged “flux tube”, which influence the surface dynamics via the Lorentz force, especially magnetic tension (Cheung et al., 2010).

Although the strong magnetic fields in sunspots inhibit convection, they do not shut it down entirely. Umbral convection is observed as umbral dots and has been simulated by Schüssler and Vögler (2006). In such strong fields, convection manifests itself as very narrow upflow plumes of hot plasma with neighboring, narrow cool return downflows. The tendency of the magnetic field to expand as the gas pressure declines toward the surface pinches the rising plumes and accelerates the upward flow. As in normal convection, the upflows are braked rapidly near the surface where the plasma loses buoyancy due to radiative cooling. The plasma piles up, the gas pressure increases and makes the plasma expand latterly, which reduces the magnetic field strength. As a result of the enhanced density, the optical depth increases and photons can only escape from higher, cooler layers producing a dark lane through the bright umbral dot (Figure 32). Above the plume, which has been decelerated, the magnetic field again closes in, arching over the gap in a cusp shape.

Figure 32
figure 32

Vertical slice through an umbral dot. Image is density fluctuation with respect to the surroundings. The solid line is Rossland optical depth unity. The dotted lines are isotherms. The arrows are velocity (longest is 2.7 km/s). Image reproduced by permission from Schüssler and Vögler (2006), copyright by AAS.

Heinemann et al. (2007), Scharmer et al. (2008), and Rempel (2011) have modeled sunspot penumbra (Figures 31, 33, 34). They find that penumbra are produced by overturning convective motions that occur in an inclined magnetic field and that the observed Evershed outflows (Evershed, 1909) are the horizontal flow of overturning convection channeled along the penumbral magnetic field.

Figure 33
figure 33

Penumbral field lines color coded by velocity at several locations in a penumbra. Set 1 is in the umbra just adjacent to the penumbral head. Sets 2 – 6 are at increasing radial distance in the penumbra. Velocities are 0 – 2 km/s (red), 2 – 4 km/s (yellow), 4 – 8 km/s (green), > 8 km/s (blue). The cutting plane shows the vertical field magnitude at 〈τ〉 = 1. Image reproduced by permission from Rempel (2011), copyright by AAS.

Figure 34
figure 34

Work as a function of radius and depth in the penumbra. The contours are where the average outflow is more than 2 and 4 km/s. The Lorentz force accelerates the fluid in the Evershed flow along the penumbral filaments, while in the inner penumbra below the surface there is an approximate balance between pressure and Lorentz forces. Image reproduced by permission from Rempel (2011), copyright by AAS.

However, unlike normal convection, it is the pressure force that that is pushing the upflow as well as the overturning horizontal flow. The downflows are in nearly hydrostatic equilibrium. Near the lower pressure photosphere the nearly vertical field lines of the penumbra spread outward (tending toward a potential field structure) and become more horizontal. The mass loading by the overturning, horizontal convective motions bends the magnetic field lines downward even more, when the initial inclination is more than 45°, which produces the nearly horizontal penumbral field (Figure 33). Cooling near τ = 1 increases the plasma density and field line bending. The Lorentz force turns the flow along the magnetic field to produce the Evershed outflow (Figure 34).

5 The Future

Rapid progress is currently occurring in solar magneto-convection simulations. What can we expect in the near future? More physics will be included: more accurate representation of the frequency dependence of the opacity, non-equilibrium ionization, partial ionization and non-LTE radiation. Such work is already begun in the bifrost (Gudiksen et al., 2011; Martínez-Sykora et al., 2012), stagger and MuRAM (Cheung and Cameron, 2012) codes. The next big step for numerical simulations of the upper photosphere and chromosphere is the inclusion of time-dependent and partial ionization effects (generalized Ohm’s Law and ambipolar diffusion) and the extension from single to multifluid equations (neutrals, ions and electrons) (Khomenko and Collados, 2012; Cheung and Cameron, 2012). The most time consuming part of “realistic” convection calculations is the radiative transfer, even with the drastic approximations currently made. Alternative numerical solutions of the transfer equation are possible (Hayek et al., 2010). Observations from larger ground (GREGOR, the Big Bear New Solar Telescope and the Advance Technology Solar Telescope) and new space based observatories (Solar Dynamics Observatory, Interface Region Imaging Spectrograph (IRIS) and Solar Orbiter) will allow more detailed comparisons between observations and simulations, which will assist in clarifying the significant physical processes that determine the solar magneto-dynamics.

The biggest unanswered questions are: exactly how does the solar dynamo work, the details of the process of mass and energy transport through and energy dissipation in the chromosphere and corona, and the origins of eruptive events. We have qualitative models of these processes. We now need a more quantitative understanding. We would like to know: how the large scale regularities of the solar cycle are produced, the relation between global and local surface dynamo action, the origin of supergranulation, the role of weak fields in energizing the chromosphere and corona, the triggers of eruptive events, and the relation between global and local coronal behavior.