NUMERICAL SUNSPOT MODELS: ROBUSTNESS OF PHOTOSPHERIC VELOCITY AND MAGNETIC FIELD STRUCTURE

Numerical simulations of sunspot fine structure have seen substantial improvements over the past five years. Schüssler & Vögler (2006) presented the first realistic MHD simulations of umbral dots, followed by a series of simulations in "slab-geometry" (Heinemann et al. 2007; Scharmer et al. 2008; Rempel et al. 2009b) addressing the transition from umbra to penumbra and gave the first insight into the origin of penumbra fine structure and the Evershed flow. Kitiashvili et al. (2009) focused on a magneto-convection simulation in strongly inclined field representative of penumbral conditions and studied the influence of field strength and inclination on the presence horizontal outflows. Rempel et al. (2009a) and Rempel (2011a, 2011b, 2011c) focused on full circular sunspots with increasing realism, with the best resolved simulations such as those shown in Rempel (2011a) reaching and exceeding the richness of detail seen in the best available observations today. As described in detail in Rempel (2011a, 2011b) these simulations reproduce the key aspects of penumbral fine structure as inferred from high-resolution observations (see, e.g., Scharmer et al. 2002, 2007; Langhans et al. 2005, 2007; Rimmele & Marino 2006; Ichimoto et al. 2007a, 2007b; Rimmele 2008; Franz & Schlichenmaier 2009; Bellot Rubio et al. 2010; Franz 2011) and recent reviews by Solanki (2003), Thomas & Weiss (2004, 2008), Rempel & Schlichenmaier (2011), and Borrero & Ichimoto (2011).

Simulations are still far from perfect. The numerical grid resolution might be marginal for resolving relevant details of the penumbra, the initial condition in terms of a monolithic axisymmetric field is very artificial, and there is a strong influence from the (unavoidable) imposed boundary conditions. Therefore, in this paper we investigate numerical resolution and boundary conditions. We also deepen our discussion on the flow field and magnetic structure in the observable layers of the penumbra. At the bottom boundary most of the sunspot simulations impose a fixed magnetic field to prevent a rapid decay of the sunspot. Rempel (2011c) showed that this constraint can be relaxed in sufficiently deep domains where convective timescales reach several days, but not in shallow ∼6 Mm deep domains typically used in simulations that address sunspot fine structure. Here we investigate how the magnetic top boundary condition as well as grid resolution influences details of penumbral fine structure, while we keep the initial state as well as the bottom boundary unchanged. We deepen the analysis presented in Rempel (2011b) to test the numerical robustness of their findings. In addition, we expand the discussion to address three aspects that have been subject to major controversy in recent years.

• 1.  
  The convective nature of the penumbra. Vigorous convection (about half the strength of granulation) is the key process in MHD simulations leading to penumbral fine structure. Indirect observational evidence for overturning motions is derived from "twisting" motions (Ichimoto et al. 2007a; Bharti et al. 2010) or correlation tracking that shows the convergence of tracers toward the edge of filaments similar to granulation (Márquez et al. 2006). Zakharov et al. (2008) found evidence for convective roles in filaments oriented parallel to the solar limb and estimated that the observed velocities are sufficient to explain the penumbral brightness. Direct evidence for convective motions has been questioned by some authors (Franz & Schlichenmaier 2009; Bellot Rubio et al. 2010), while it was found by others in high-resolution observations (Sánchez Almeida et al. 2007; Joshi et al. 2011; Scharmer et al. 2011; Scharmer & Henriques 2011). Joshi et al. (2011) found convective downflows mostly in the inner penumbra using the rather deep forming CI 5380 line. Scharmer et al. (2011) and Scharmer & Henriques (2011) found consistent results using the CI 5380 and Fe i 6301 lines in terms of an anisotropic convection pattern with a correlation between intensity and vertical flow velocity not very different from quiet Sun. Such correlations were also previously indicated in the analysis of Sánchez Almeida et al. (2007).
• 2.  
  The presence of opposite polarity magnetic flux throughout the penumbra and its relation to downward directed mass flux. While inverse polarity flux is certainly present in the outer penumbra (del Toro Iniesta et al. 2001), it is debated if this is also the case in the inner penumbra. Sánchez Almeida (2005) found substantial amount of opposite polarity flux throughout the penumbra, accounting to about 58% of the unsigned magnetic flux of the penumbra, Langhans et al. (2005) did not find strong evidence in magnetogram data, possibly due to resolution effects. Recently, Franz (2011) found with Hinode data that about 40% of downflows in the outer penumbra are in regions with opposite polarity flux and consider this a lower limit.
• 3.  
  The magnetization of the Evershed flow. Spectropolarimetric observations point to a flow in regions with a substantial magnetic field (see, e.g., Ichimoto et al. 2008; Borrero 2009; Borrero & Solanki 2008, 2010); see also discussions in Brummell et al. (2008) and Thomas (2010). Some models suggest that penumbral filaments are mostly field free "gaps" (Spruit & Scharmer 2006; Scharmer & Spruit 2006). A similar conclusion was also derived from an interpretation of "twisting filament motions" as fluting instability (Spruit et al. 2010), leading to an upper bound for magnetic field in flow channels of about 300 G. MHD simulations however point toward strong 1–2 kG field in flow channels, in the outer penumbra the field is even enhanced in flow channels compared to the background (Rempel 2011b).

The numerical models were computed with an expanded version of the MURaM radiative MHD code (Vögler et al. 2005). We introduced an artificial limitation of Alfvén velocities to 60 km s⁻¹ in order to prevent severe Courant–Friedrichs–Lewy time step constraints and implemented a new numerical diffusivity approach; both are described in Rempel et al. (2009b). For the simulations presented here with grid sizes of up to 4.8 billion grid points we did additional performance enhancements including a major rewrite of Input/Output routines and improved the scalability of the code to up to 24, 576 cores. Most of the simulations presented here were computed in the 1024–12, 288 core range.

We present here a series of numerical sunspot models which were all computed in a domain with a size of 49.152 × 49.152 × 6.144 Mm³. All models were initialized with an axisymmetric self-similar magnetic field configuration described in Appendix A. For the models presented here we used the parameters B₀ = 6.4 kG, R₀ = 8.2 Mm, and z₀ = 6.4185 Mm, leading to a sunspot with an initial flux of Φ₀ = 1.2 × 10²² Mx, a field strength of B₀ = 6.4 kG at the bottom of our domain (z = 0 Mm), dropping to 2.56 kG at the top (z = 6.144 Mm).

While keeping this initial condition fixed we explore the dependence of the resulting sunspots on the top boundary condition and grid resolution.

We use in the horizontal direction periodic boundaries; the bottom boundary condition is open for convective flows in regions with |B| < 2.5 kG, but closed otherwise. The open boundary condition imposes a symmetric mass flux (all three components) in the ghost cells and extrapolates the gas pressure such that its value at the boundary is fixed. In outflow regions the entropy is determined by the upstream values from within the computational domain, in inflow regions the entropy is set to a fixed value that leads to the correct solar energy flux under quiet-Sun conditions. The closed boundary condition in regions with |B| > 2.5 kG is implemented through an antisymmetric mass flux. The top boundary condition is closed (vertical mass flux antisymmetric) and stress-free for horizontal motions.

For the magnetic field we use a boundary condition that allows us to manipulate the inclination angle of the magnetic field. We focus on a single sunspot, which implies that the horizontal boundary imposes the same polarity spots nearby and one thus produces only a very subdued penumbra (Rempel et al. 2009a). A relaxation of the horizontal periodicity is non-trivial and we thus decided to use a boundary that allows us to control the inclination angle of the magnetic field to mimic different global field configurations, which is described in detail in Appendix B. A free parameter, α, gives for α = 1 a potential field extrapolation subject to horizontal periodicity, i.e., the field becomes vertical asymptotically. Values of α > 1 lead to more horizontally inclined field, i.e., the horizontal field component is approximately enhanced by a factor of α compared to the potential field reference. The resulting field extrapolations are for α ≠ 1 not force-free outside the computational domain. They are not intended to give a realistic description for the global field topology above a sunspot, but rather to explore the influence of different field geometries on the structure of a sunspot penumbra.

All of the simulations presented here were computed with gray radiative transfer. To allow for a better comparison with observations through forward modeling of spectral lines we computed two non-gray models that are derived from the gray simulations. We have a non-gray version of the simulation with 32 [16] km grid resolution that was advanced for 26 minutes and we have a non-gray model with 12 [8] km resolution that was restarted from the upper half of our 16 [12] km case and advanced for 15 minutes. The non-gray model with 32 [16] km was analyzed by Bharti et al. (2011); more results will be presented in forthcoming publications.

In the following analysis, we will focus entirely on the gray simulations and discuss penumbral structure mostly through quantities extracted on constant τ surfaces. In Figure 1 we display a snapshot from our highest resolution gray simulation (16 [12] km) that we consider in this paper. The simulation was computed with the α = 2 boundary condition and evolved for 1 hr at the highest resolution.

Zoom In Zoom Out Reset image size

3.1. Influence from Magnetic Top Boundary Condition

In Figure 2 we show four sunspot models computed with the boundary conditions α = 1, 1.5, 2, 2.5. The case α = 1 (lower right quadrant) leads to a very subdued penumbra with a few isolated filaments, which is consistent with the penumbra structure that was found in Rempel et al. (2009a) in the direction where the periodicity imposes the same polarity sunspots. The other three cases show more developed penumbrae; the overall radial extent increases with the value of α. In Figure 3 we show azimuthal averages of (a) intensity, (b) Evershed flow, (c) vertical and radial magnetic field, and (d) inclination angle. These quantities are in addition averaged for about 1 hr in time. Apart from the monotonic increase of the penumbra extent with α, it is also evident that the most dramatic differences occur between the potential field case (α = 1) and the other three cases with α > 1. For α = 1 the azimuthally averaged radial flow velocity stays around 500 m s⁻¹, while all other cases have outflows with more than 3 km s⁻¹ on average. Similarly, the curves for the radial magnetic field as well as inclination angle do not differ as much between the α > 1 cases as compared to α = 1 and the rest. While the potential field case reaches only an average inclination of 50°, all the other solutions with more extended penumbra and Evershed flow reach about 65°.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This leads to the surprising conclusion that the potential field case (which is the only physical, i.e., force-free, solution outside the computational domain) is a clear outlier compared to the rest. This result has to be seen in the context of horizontal periodicity, which is not the proper horizontal boundary condition to describe the magnetic field above sunspots. A more reasonable (but computationally more difficult to implement) boundary condition would be the free expansion of magnetic field into a half-room (perhaps considering in addition spherical geometry), which would lead naturally to a magnetic field with stronger horizontal field components. The boundary conditions with α > 1 emulate that behavior while maintaining horizontal periodicity. Note that simply increasing the horizontal extent of the domain for α = 1 does not alleviate this problem, since the field remains asymptotically vertical. For the following investigation we will use a top boundary condition with α = 2.

3.2. Influence from Numerical Grid Resolution

We modify the numerical grid resolution in the range from 96 [32] to 16 [12] km (horizontal [vertical] resolution). The sunspot models computed with different grid resolutions are not fully independent from each other. The model with a resolution of 48 [24] km resolution was started from a snapshot of the 96 [32] km resolution simulation evolved for 1 hr and ran for an additional 5.5 hr. The model with a resolution of 32 [16] km resolution was started from a snapshot of the 48 [24] km resolution run after 3.3 hr and ran for an additional 2.2 hr. The highest resolution case with 16 [12] km resolution was started from the last snapshot of the 32 [16] km case and evolved for an additional hour. Overall, we compare sunspot models that have been evolved for about 6 hr from the initial state. We did not evolve the higher resolution cases for the full length of time because of their computational expense (1 hr at 16 [12] km resolution, 3072 × 3072 × 512 grid points, costs about 800,000 CPU hours on a CRAY XT-5). For the small length and associated short timescales in the photosphere about 1 hr at the highest resolution is sufficient to allow for the solution to adapt to the grid spacing. Structures in the deeper parts of the domain are already well resolved at lower resolution.

Figure 4 shows the (gray) intensity and magnetogram at τ = 1; Figure 5 shows the corresponding azimuthal averages of quantities at τ = 1. Overall there is no significant influence from grid resolution on the radial extent of the penumbra as well as the magnetic structure at τ = 1, displayed in panels (c) and (d). Significant differences occur with regard to the average intensity profile and Evershed flow shown in panels (a) and (b). Only the higher resolution cases show the formation of a plateau-like feature at around 0.7 I_☉ in the inner penumbra. The average Evershed flow speed is increasing monotonically with resolution from a peak at around 2.3 km s⁻¹ for 96 [32] km resolution to about 4.6 km s⁻¹ for 16 [12] km resolution. However, we see also here a clear sign of convergence: while the relative resolution changes from 96 [32] to 48 [24] km and from 32 [16] to 16 [12] km are the same, the relative changes in the peak flow velocities are 1.55 and 1.1, respectively. The strong dependence of the Evershed flow velocity on grid resolution is a consequence of the driving mechanism that is concentrated in a thin boundary layer just beneath τ = 1 (Rempel 2011b), which is difficult to resolve numerically (see also Section 7). The umbra is, with 0.3 I_☉, rather bright for a sunspot with an umbral field strength exceeding 3 kG. We return to this point in Section 8.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 presents the penumbral fine structure at the τ = 1 level for our highest resolution case (16 [12] km resolution). This figure is very similar to Figure 3 in Rempel (2011b), in which we displayed the fine structure at 32 [16] km resolution for the double sunspot simulation of Rempel et al. (2009a). The large degree of similarity underlines that the details of fine structure are not very sensitive to the numerical setup. Penumbral filaments show an enhancement of the radial magnetic field component (panel (b)), while the vertical field component is strongly reduced (panel (c)). In the outer penumbra we find a substantial amount of opposite polarity flux, which will be characterized further in Section 4.3. The combination of an enhanced radial magnetic field component with a reduced vertical field component results in the uncombed structure of the penumbra with a strong variation of the field inclination angle (panel (d)). Fast horizontal (Evershed) outflows are found along the almost horizontal stretches of magnetic field (panel (e)). Overturning convection (panel (f)) is the underlying driving mechanism, which will be characterized further in Section 5. We provide an animation of Figure 6 in the online version of the journal.

Zoom In Zoom Out Reset image size

Figure 7 presents the height dependence of the azimuthally averaged Evershed flow and the vertical rms velocity. As already pointed out in Rempel (2011b), the Evershed flow peaks in the deep photosphere and falls off rapidly with height; we see the transition to an inverse Evershed flow at about τ = 0.01. This figure is qualitatively very similar to the observations reported by Bellot Rubio et al. (2006; see Figure 8 therein) and by Borrero et al. (2008). In Figure 7(b) we show the height dependence of the vertical rms velocity. We find essentially a drop by a factor of two between the τ = 1 and τ = 0.1 levels, but no further drop toward higher layers.

Zoom In Zoom Out Reset image size

4.1. Robustness of Penumbral Fine Structure

To characterize fine structure further and quantify the resolution dependence of quantities, we present in Figure 8 correlations between the Evershed flow velocity and (a) intensity, (b) magnetic field strength, as well as (c) radial and (d) vertical magnetic field components. All correlations are computed based on fluctuations around a smooth background that was obtained through a convolution with a Gaussian having an FWHM of 1.5 Mm. At τ = 1 the Evershed flow is found in the inner penumbra in bright features, whereas the correlation vanishes toward the outer penumbra (panel (a)). A similar behavior was also found by Schlichenmaier et al. (2005) and Ichimoto et al. (2007b). The correlation between Evershed flow and magnetic field strength is negative in the inner half, but positive in the outer half of the penumbra. An Evershed flow in stronger magnetic field regions in the outer penumbra was also inferred by Tritschler et al. (2007) and Ichimoto et al. (2008) from net circular polarization variations with the viewing angle. The sign change in the v_R − |B| correlation results from two effects: a positive correlation with the radial field component and a negative correlation with the vertical field component throughout the penumbra (Figures 8(c) and (d)). While the negative correlation with the vertical field dominates the inner penumbra, the positive correlation with the radial field component dominates the outer penumbra (see also Figure 6). Except for the correlation between v_R and B_z that is quite significant with values about −0.7, all others reach only moderate levels of up to 0.5 in the deep photosphere. All correlations decay very rapidly with height.

Zoom In Zoom Out Reset image size

To quantify the robustness of photospheric fine structure, we focus now on the correlations at τ = 1 and study the resolution dependence in Figure 9. The lowest resolution case (96 [32] km) somewhat misses the v_R − I correlation in the inner and the v_R–|B| correlation in the outer penumbra; all other grid resolutions produce comparable levels of correlations. We can conclude that the magneto-convection process underlying penumbral fine structure is well captured starting from a resolution of about 48 [24] km. This resolution also produces an average Evershed flow amplitude at this resolution that differs less than 20% from our highest resolution case. As we will discuss in Section 7, the trend of increasing magnetic field in flow channels (see v_R − B_R correlation) and the resolution dependence of the Evershed flow have a common origin.

Zoom In Zoom Out Reset image size

4.2. Magnetized Flow Channels

The positive v_R − B_R correlation throughout the penumbra increases in amplitude with increasing resolution, which strongly points toward an active process supporting strong horizontal field in the flow channels as opposed to remnant magnetic field due to unavoidable numerical diffusion effects as speculated by Nordlund & Scharmer (2010). It was pointed out by Rempel (2011b) that the horizontal field originates mostly from the horizontal shear in the subphotospheric Evershed flow, i.e., the term B_z∂_zv_R in the induction equation (see Section 7). The robustness of this effect points toward a strongly magnetized Evershed flow in photospheric layers as it has been inferred from spectropolarimetric observations (see, e.g., Ichimoto et al. 2008; Borrero & Solanki 2008, 2010; Borrero 2009). This result is not necessarily in contradiction to the "gappy" penumbra model of Spruit & Scharmer (2006) and Scharmer & Spruit (2006) since the strong horizontal magnetic field is confined to the thin boundary layer at τ = 1 (Rempel 2011b). The value of |B_R| is typically about a few 100 G stronger in flow channels than in the background. In Figure 10 we present a bivariate probability density function (PDF) for v_R and |B| in the inner penumbra (R = 10–13 Mm) in panel (a) and the outer penumbra (R = 13–16 Mm) in panel (b). In the inner penumbra we find field strength of up to 4 kG in regions with low radial velocity, which correspond to the spines. Fast outflows have field strength of more than 2 kG, mostly due to B_R. We see a similar behavior in the outer penumbra with reduced overall field strength and less pronounced spines, but also here fast outflows are associated with field strengths around 2 kG. The extension toward low field strength is caused by a few granules present in the region 13 Mm < R < 16 Mm. In both, inner and outer penumbrae, there is also a population of inflows with substantial field strength. We see also a trend for increasing |B| with increasing outflow velocities (v_R > 0).

Zoom In Zoom Out Reset image size

4.3. Inverse Polarity Flux

Inverse polarity magnetic flux embedded in the sunspot penumbra is an integral part of penumbral fine structure. Due to the substantial amount of overturning convection (see Section 5) a certain amount of magnetic field lines has to turn back into the photosphere within the penumbra. Return flux along the outer penumbra rim was first found by Westendorp Plaza et al. (1997) and interpreted in the context of predictions from flux tube models such as Thomas & Montesinos (1993) and Montesinos & Thomas (1997). Fast (supersonic) downflows and return flux in the interior of the penumbra were found by del Toro Iniesta et al. (2001). Substantial amounts of inverse polarity flux were deduced from observations by Sánchez Almeida (2005), which we will use as a reference for the analysis presented here. This subject is not without controversy; for example, Langhans et al. (2005) did not find much evidence in magnetograms at high spatial resolution.

In Figures 11 and 12 we quantify the amount of inverse polarity flux found in the penumbra. Following Sánchez Almeida (2005), we show in the left panels of both figures the quantities R〈|B_z|〉 (solid) and R〈B_z〉 (dashed) as a function of radius, where 〈...〉 indicates the azimuthal average. The right panels show the ratio of the radially integrated signed and unsigned magnetic flux. Figure 11 shows how these quantities change with optical depth in the range from τ = 1 to τ = 0.001. Significant amounts of inverse polarity flux are only present in the deep photosphere, where we find integrated over the whole sunspot up to 11%, already at τ = 0.1 only half of that remains. Figure 12 shows how these quantities vary with grid resolution as well as boundary conditions. Comparing the integrated values at R = 18 Mm in panel (b), we see an increase of the inverse polarity flux from abut 8% to 11% over the resolution range considered here (about half of the spread is due to our lowest resolution run). We also find that more extended penumbrae tend to have more opposite polarity flux than less extended ones (yellow dashed/dotted lines). Given the fact that the effective extent of the penumbra computed with α = 1.5 (2.5) is 16 (19) Mm, the corresponding fractions of opposite polarity fluxes are 6 (12)%. For a better comparison with the results from Sánchez Almeida (2005), we present in Table 1 quantities that are comparable to those in their Table 1. Here the quantity Φ_z is defined analogously to theirs as

$$\begin{equation} \Phi _z(R_1<R<R_2)=2\pi \int _{R_1}^{R_2}R \langle B_z\rangle {d}R\;. \end{equation} \tag{ 1 }$$

To compute the values in Table 1, we used R_s = 18 Mm for the radius of the sunspot. We present here values computed for the highest resolution case. Our simulated sunspot is about a factor of 1.6 larger in terms of the unsigned flux and a factor of two larger in terms of the signed flux. While we find a large degree of qualitative agreement with Sánchez Almeida (2005), the overall amount of inverse polarity flux falls short by a factor of about three (first line in the table). We find a substantial amount of opposite polarity flux in the region 12.7 Mm < R < 18 Mm, where the signed flux is 50% of the unsigned flux. In contrast to Sánchez Almeida (2005) we do not find opposite polarity flux at photospheric levels in the umbra (although, some amount of opposite polarity field is present in umbral dots beneath the τ = 1 surface in our simulations).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Quantification of Inverse Polarity Flux

Location	Description	Unsigned	Signed
Φz(0 < R < Rs)	Full spot	Φ0 = 1.18 × 1022 Mx	0.89 Φ0
Φz(Rs/2 < R < Rs)	Penumbra	Φ1 = 0.44 Φ0	0.76 Φ1
$\Phi _z(R_s/\sqrt{2}<R<R_s)$	Outer penumbra	Φ2 = 0.19 Φ0	0.49 Φ2
Φz(0 < R < Rs/2)	Umbra	Φ3 = 0.54 Φ0	1.00 Φ3

Download table as:  ASCIITypeset image

The fact that the amount of return flux in simulations shows only little resolution dependence might come as a surprise, since most of this flux is present at small scales which vary somewhat with grid resolution. However, return flux is a consequence of vigorous overturning convection, which itself shows only little resolution dependence (see Section 5). Overturning convection is a direct consequence of energy flux constraints that are the same in all cases considered regardless of resolution.

4.4. Relation between Opposite Polarity Regions and Downflows

Figure 13(a) presents the filling factors of regions with opposite polarity field B_z < −25 G (blue) and B_z < −250 G (red). Both peak at about R = 16 Mm, which has also been indicated in Figure 6 through the dotted circle. The peak filling factors are 0.4 and 0.2 for these two thresholds. The filling factor of supersonic downflows (green) reaches about 1%. In Figure 13(b) we evaluate the downward directed mass flux in these regions relative to the total downward directed mass flux as a function of radius. Up to 60% (44%) of the returning mass flux is found in regions with B_z < −25 G (B_z < −250 G), up to 5% in supersonic downflows. Integrated over the entire penumbra (from R = 8 to R = 18 Mm), about 40% (27%) of the downward directed mass flux is found in regions with B_z < −25 (B_z < −250). Supersonic downflows contribute only 2%–3%. Solid, dotted, and dashed lines indicate the respective quantities for the resolution levels of 16 [12], 32 [16], and 48 [24], respectively.

Zoom In Zoom Out Reset image size

At τ = 1 the average velocity in the supersonic downflow regions is −9.6 km s⁻¹, with the fastest flows reaching −15 km s⁻¹. The average magnetic field strength is 2.7 kG. The average gas pressure is, with 6.8 × 10⁴ dyne cm⁻², about 30% lower than the typical photospheric pressure (at a constant height level the difference would be even larger), leading to an average intensity that is clearly enhanced with 1.1 I_☉, in particular for a strong downflow. This has consequences for the I–v_z correlation in the outer penumbra, which we will discuss further in Sections 5 and 6.

The fact that most of the returning mass flux is not associated with opposite polarity magnetic flux indicates that it returns mostly by submerging magnetic field rather than by flowing along downward directed field lines. To quantify this aspect further, we present in Figure 14(a) the inclination of flows and magnetic field relative to the horizontal radial direction, i.e., the quantities $\gamma _v=\arctan (v_z/v_R)$ and $\gamma _B=\arctan (B_z/B_R)$. We average these expressions in Evershed flow channels separated into upflow and downflow regions, i.e., we consider regions with v_R > 0 and v_z > 250 or v_z < −250 m s⁻¹. In upflow regions both field and flow show a very similar inclination dropping from about 30° at R = 8 Mm to 10°–15° at R = 15 Mm. For R > 15 Mm the inclination for the magnetic field stays at around 10°, while the flow inclination increases to 40°. The situation is different in downflow regions. Values for the flow inclination are −40° in the inner and outer penumbra, while in the center of the penumbra values around −20° are typical. The magnetic field starts with an angle of 35° at R = 8 and continues to point upward to about R = 14.5 Mm. The minimum inclination the field reaches is −10° at R = 16 Mm. Except for the cores of penumbral filaments (equivalent to upflow regions) in the innermost 2/3 of the penumbra, the inclination angles of flow and field are quite different. Flows in the penumbra are not just flows along the magnetic field, in particular not in downflow regions. To illustrate this point more clearly, we present in Figure 14(b) the average angle between the flow and field given by $\arccos (\vert \vec{v}\cdot \vec{B}\vert /(\vert \vec{v}\vert \, \vert \vec{B}\vert))$. In the innermost 2/3 of the penumbra, the angle between flow and field is 10°–20° in upflow and 25°–60° in downflow regions. These values are larger than those evident from panel (a) since we measure here the average of the local misalignment, which is different from the misalignment of the average field and flow. The inclination of flow and field agrees qualitatively with the findings of Scharmer et al. (2011, supporting online material). They found that in locally bright features (upflows) inclination angles of flow and field agree very well, while the magnetic field was close to horizontal or weakly upward pointing in downflow regions (locally dark features).

Zoom In Zoom Out Reset image size

Franz (2011) concluded from Hinode data that at least 40% of penumbral downflows contain magnetic field with opposite polarity, similar to our finding. This agreement might be coincidental as a substantial fraction of upflowing and downflowing mass flux still might be hidden in observations. Franz (2011) speculated that this is a lower limit and possibly all downflows have opposite polarity flux. We see only a moderate increase of this fraction with resolution (for the resolutions of 48 [24], 32 [16], 16 [12] km, we find 33%, 36%, 40%, respectively). In the inner penumbra, mass is returning beneath the photosphere by submerging entire field lines (that still connects to the upper boundary) along the lateral downflow lanes of filaments without requiring opposite polarity flux to form (see also Figure 19 in Rempel 2011b). We see a moderate increase in this number with the overall extent of the penumbra for a fixed resolution of 48 [24] km from 29% to 36% (comparing the α = 1.5 and α = 2.5 cases). If we extrapolate our results assuming that 100% of downflows have opposite polarity flux (for which we do not see any indication here), we would obtain values similar to those found by Sánchez Almeida (2005).

From our models we expect that more than half of the returning mass flux in the penumbra is found in regions that have the same polarity as the umbra of the sunspot. Examples for such downflows were recently detected by Katsukawa & Jurčák (2010), although their connection to overturning motions is not evident from their observations.

4.5. Mass Flux Associated with Large-scale Flow Component

The mass flux in the penumbra can be formally separated into a large-scale mean flow with upflows in the inner penumbra and downflows in the outer penumbra and laterally overturning motions. The former corresponds in the photosphere to the Evershed flow. Note that this separation is only meaningful in terms of the overall mass flux balance. Individual flow elements always have a combination of radial outflow and lateral overturning motions. To quantify both components of the mass flux we define (see also Rempel 2011b) m⁺_z = 〈v_zϱ > 0〉 and m⁻_z = 〈v_zϱ < 0〉, where 〈...〉 denotes the azimuthal average. The unsigned vertical mass flux is then given through m_tot = m⁺_z − m⁻_z, while the locally unbalanced mass flux associated with radially separated regions of up and downflows (and corresponding horizontal flows in between) is given by m_mean = m⁺_z + m⁻_z. As a relative measure for the mass flux contained in the azimuthal mean flow we consider

$$\begin{equation} \varepsilon =\frac{\int _{R_1}^{R_2}R\big| m_z^++m_z^-\big| {d}R}{ \int _{R_1}^{R_2}R\big| m_z^+-m_z^-\big| {d}R}. \end{equation} \tag{ 2 }$$

Evaluating this expression on a constant height surface about 350 km beneath the τ = 1 level in the plage region (to make sure we stay below τ = 1 in the inner penumbra), we obtain with R₁ = 6 and R₂ = 18 Mm values of 15%–16% for grid resolutions from 48 [24] km to 16 [12] km. Similarly, we obtain values within this range for the two simulations at 48 [24] km resolution with different extents of the penumbra due to changes in the top boundary condition. Here we used a value of R₂ = 16 Mm and R₂ = 19 Mm for the α = 1.5 and α = 2.5 cases to consider the different extent of the region occupied by the Evershed flow according to Figure 3. This is also in agreement with the value found previously by Rempel (2011b) for the double sunspot simulation of Rempel et al. (2009a).

In order to evaluate this quantity on a constant τ surface, we have to replace m_z with the mass flux normal to τ levels and include a geometric factor considering projection effects due to inclined τ surfaces. Both effects are considered by using the expression $m_{\tau }=\vec{m}\cdot \nabla \bar{\tau }/\vert \partial _z\bar{\tau }\vert$ instead of m_z at the respective τ level. Here $\bar{\tau }$ denotes a horizontally smoothed τ, since we did not find a sufficiently balanced mass flux if we use τ at maximum resolution. We compare here results that were obtained after convolving τ in the horizontal direction with a Gaussian having an FWHM of 192 km. We again find comparable values for the resolution levels from 48 [24] km to 16 [12] km. Values of ε for the τ levels of 1.0 (0.1) are 14%–16% (8%–11%) (note that this is relative to the unsigned total mass flux at the respective level, which is at the τ = 0.1 level about 25% of that at τ = 1). The value at τ = 1 is very close to what we found on a constant height surface about 300 km deeper, pointing toward robustness of ε with regard to the height (or τ level) as well as grid resolution.

With the advent of magneto-convective models of sunspot penumbrae, the role of overturning convection in the penumbra has been the subject of a controversial debate.

On the theoretical side, overturning convection is the key process responsible for energy transport and filamentation (see, e.g., Spruit & Scharmer 2006; Scharmer & Spruit 2006; Heinemann et al. 2007; Scharmer et al. 2008; Rempel et al. 2009a, 2009b; Kitiashvili et al. 2009; Rempel 2011b).

On the observational side, the evidence for overturning convection is controversial. While the existence of downflows and upflows is evident (see, e.g., the direct evidence by del Toro Iniesta et al. 2001; Franz & Schlichenmaier 2009), their nature is debated: do upflows and downflows correspond to the radial endpoints of penumbral flow channels and are just the vertical component of the Evershed flow, or do individual penumbral filaments show a flow pattern similar to that of granules with a substantial amount of mass flux turning over in lateral directions? With regard to the latter, direct and indirect evidence has been presented by Márquez et al. (2006), Sánchez Almeida et al. (2007), Ichimoto et al. (2007a), Rimmele (2008), Zakharov et al. (2008), Bharti et al. (2010), Joshi et al. (2011), Scharmer et al. (2011), and Scharmer & Henriques (2011). Other investigations claim that these motions do not exist (Franz & Schlichenmaier 2009; Bellot Rubio et al. 2010).

Figure 15 quantifies the vertical rms velocity and the relative rms intensity contrast for different grid resolutions. In panel (a) the dashed line indicates an approximate relation of the form

$$\begin{equation} I\propto \sqrt{v_z^{\rm rms}(\tau =1)}\,, \end{equation} \tag{ 3 }$$

which was already found by Rempel (2011b). This relation means that about 50% of the level of overturning convection found in the plage surrounding the sunspot is required to maintain an intensity of about 0.7 I_☉ typical for the inner penumbra. This is possible since the rms contrast is about a factor of 1.5 larger with 25% than the contrast in the plage (17%). Note that the contrast values are a few percent higher in our gray simulations compared to a non-gray run.

Zoom In Zoom Out Reset image size

We find again only a very moderate dependence on the grid resolution as well as on penumbra extent (dashed and dotted yellow lines). Only the contrast of the lowest resolution case shows significant deviations in regions with I < 0.7 I_☉. If there is any systematic variation at all, it is a very moderate increase of the rms velocity with increasing grid resolution. Overall, the predictions about the amplitude of overturning convection in the penumbra are very robust. Lower values for the amplitude of overturning motions would be only possible in combination with a substantially higher contrast to satisfy energy flux constraints.

While the overall amount of overturning convection as characterized by the vertical rms velocity is robust, there are differences in how much upflow and downflow contribute. To characterize their contribution further, we show in Figure 16(a) the filling factors of regions with more than 250 m s⁻¹ upflow (blue) and less than −250 m s⁻¹ downflow velocity (red). Filling factors of upflow and downflow are comparable in the plage region outside the sunspot with values of about 45%. The outer penumbra shows a clear preference for downflows with a filling factor of 55%, while upflows occupy only 30% of the area. In the inner penumbra the upflow filling factor exceeds the downflow filling factor in the highest resolution case and is comparable with the other cases. In Figure 16(b) we show the respective average velocities in upflow and downflow regions. In the plage region downflows dominate over upflows by a factor of about 1.5. This is mostly an optical depth effect since we show quantities on the τ = 1 surface. Since upflow and downflow velocity amplitude decline with height, the elevated τ levels in upflow regions sample lower velocities than the depressed τ levels in downflow regions, while the difference on a constant height surface is typically less pronounced. This "apparent redshift" is not to be confused with the well-known convective blueshift that originates in unresolved observations from the different intensity weighting of upflow and downflow regions dependent on the temperature sensitivity of the considered spectral line. In the inner penumbra upflow, velocities dominate over downflows by a factor of up to 1.7 (despite the fact that also here the elevated τ levels in upflow regions produce a bias toward downflows). For the currently explored resolution range we do not yet see a convergence of filling factors as well as mean upflow and downflow velocities. It is likely that the trend toward upflow dominance in the inner penumbra will continue somewhat further with increasing resolution.

Zoom In Zoom Out Reset image size

In Figure 17 we present the widths of intensity, upflow and downflow features as a function of radius. To this end we determine the average azimuthal extent of regions with intensity larger than 0.05 I_☉ compared to the background that is defined through a smoothed intensity profile (Gaussian with an FWHM of 1.5 Mm). Upflow and downflow are again defined through regions with more than 250 m s⁻¹ velocity amplitude. We consider here the following measures for the widths:

$$\begin{equation} \langle s(r)\rangle _n=\frac{\sum _i s_i(r)^{n+1}}{\sum _i s_i(r)^n}\;, \end{equation} \tag{ 4 }$$

where s_i(r) are the detected widths at a given radial position. For n = 0 we obtain the regular average, while n = 2 can be considered as an "area-weighted average" if we make the assumption that the area of a feature scales with s²_i. Figure 17 shows 〈s(r)〉₀ in panel (a) and 〈s(r)〉₂ in panel (b). Regardless of the adopted measure we find that the widths of intensity and flow features are comparable in the plage region, while flow features are systematically smaller in the penumbra by a factor of up to two. We do not see a convergence of 〈s(r)〉₀ with grid resolution (we show here the 16 [12] and 32 [16] km cases by solid and dotted lines) in either penumbra or plage due to small-scale features that dominate in both regions by numbers. The area-weighted average 〈s(r)〉₂ shows better convergence, in particular for the width of intensity features. The most dramatic changes occur in flow features in the mid to outer penumbra. The latter indicates that resolution even higher than 16 km might be required to capture the width of flow features in the penumbra properly and that the current width of about 100 (200) km has to be considered as an upper limit for the average (area-weighted average) width.

Zoom In Zoom Out Reset image size

In Figure 18 we present the v_z–I correlation, which is very often used to characterize convective energy transport. The v_z–I correlation reaches in the inner penumbra values similar to those found outside the sunspot (peak at around 0.6). Interestingly, the correlation vanishes in the outer penumbra. We further find that a significant correlation is only present with the velocity at τ = 1 inside the penumbra, whereas τ = 0.1 leads also to a very moderate correlation (0.3–0.4) outside the sunspot. We find here again only little resolution dependence in the range from 48 [24] to 16 [12] km. A v_z–I correlation in the inner penumbra comparable to the plage surrounding the sunspot and a very weak correlation in the outer penumbra are captured well starting from a resolution of 48 [24] km.

Zoom In Zoom Out Reset image size

Note that the drop of the correlation in the outer penumbra is due to photospheric effects. If we compute instead the correlation between m_z and (e_int + p)/ϱ + 1/2v² (here e_int and p denote the internal energy and gas pressure) on a constant height surface located about 350 km beneath average τ = 1 in the plage region, we find a constant value of 0.65 for R > 8 Mm. This indicates that rather moderate values of the I–v_z correlation in the photosphere cannot be taken as an argument against convective energy transport, in particular not if strong magnetic field is present.

In Section 5 we discussed the correlation between intensity and vertical velocity. We found comparable values (∼0.6) in the inner penumbra and the strong plage surrounding the sunspot, while the correlation disappears in the outer penumbra. In Figure 19 we show for the highest resolution case (16 [12] km) bivariate PDFs for vertical velocity and intensity, computed (left to right) for the inner penumbra (R = 10 to R = 14 Mm), the outer penumbra (R = 14 to R = 18 Mm), plage (R > 20 Mm), as well as plage regions with |B_z| < 500 G. We use here wider masks for the inner and outer penumbrae than in Section 4.2 since we want to include all of the fast downflows at the outer edge of the penumbra. The PDFs are computed for absolute intensity values (rather than fluctuations about the background) to allow for a direct comparison of intensity values between penumbra and plage.

Zoom In Zoom Out Reset image size

Comparing first the plage region and the plage with |B_z| < 500 G, we see that the moderate correlation of 0.59 is mostly due to a population of very bright downflows that are associated with strong field. Masking out these regions increases the correlation to 0.73. The correlation for the inner penumbra has a value of 0.53 very close to that of plage; however, the shape of the PDF is substantially different. While the correlation in the plage regions results in about equal parts from bright upflows and dark downflows, it is in the inner penumbra mostly due to bright upflows, some with brightness values reaching 1.4 I_☉. Similarly to the plage region, the decorrelation in the outer penumbra is due to a strong population of bright downflows. In particular, most supersonic downflows have a brightness comparable to I_☉. Most of these features also have strong opposite polarity magnetic field, indicating that the mechanism for the brightness enhancement is likely similar in the penumbra and plage.

The less pronounced tendency for downflows to be darker in the inner penumbra is in part due to the fact that the typical width of flow features is about half the width of intensity features according to Figure 17. While the upflows are mostly located centrally in bright patches, downflows are located near the edge and are partially found in bright and dark features.

Computing the PDFs based on the vertical velocity at τ = 0.1 does not lead to any significant correlation in the penumbra, while the plage region shows values of 0.38 and 0.55 excluding strong field regions. The core of the PDFs leads in most cases to a stronger correlation. This implies that smoothing the data increases the overall level of correlation found for both penumbra and plage. For example, convolving the data by a Gaussian with FWHM of 150 km leads to correlations of 0.63, 0.19, 0.72, and 0.82 for τ = 1. These values are comparable to those recently found by Scharmer & Henriques (2011) for the interior penumbra and quiet Sun. We see a good agreement with their results overall, although differences in the shape of the I–v_z PDFs for the penumbra and plage between their observations and our simulations exist. They did find quite comparable results for the CI 5380 and Fe i 6301 lines, while we see a rather strong height dependence for the PDFs and resulting correlations between τ = 1 and τ = 0.1. The bivariate PDFs do not show substantial differences between the 48 [24] km and 16 [12] km cases in terms of the features discussed above.

PDFs for v_z, v_R, and v_Φ are presented in Figure 20. At τ = 1 the PDF for v_z peaks in the plage region near v_z = 1.5 km s⁻¹ (the "typical" upflow velocity), while the PDFs for the inner and outer penumbrae peak near v_z = 0 km s⁻¹. The PDF in the inner penumbra is skewed toward upflows; the PDF in the outer penumbra has a far extending tail of fast supersonic downflows. Toward lower τ values the shape of the PDF for all three regions becomes similar, while the PDF in the inner and outer penumbrae remains narrower compared to the plage region. Note that the fast downflow wing in the outer penumbra extends toward higher layers than the fast upflow wing in the inner penumbra, which is only present at τ = 1. A similar behavior was seen by Scharmer & Henriques (2011), who found that downflows have a similar strength in CI 5380 and Fe i 6301, while upflows are stronger in the deeper forming CI 5380 line.

Zoom In Zoom Out Reset image size

The PDF for v_R is broader in the penumbra compared to the plage at all τ levels shown. For τ = 1 and 0.1, it is skewed toward outflows, while the inflow part matches up with that of the plage region; toward τ = 0.001 it is skewed toward inflows, while the outflow part matches up that of the plage region. In the inner penumbra the PDF shows a secondary peak at τ = 1 at an outflow velocity of about 7–8 km s⁻¹, corresponding to the "typical" Evershed flow velocity within flow channels (the azimuthal average is about a factor of two smaller).

At τ = 1 the PDF for laterally overturning motions in the outer penumbra and plage are similar, while the inner penumbra shows a narrower PDF. In higher layers the PDFs for the penumbra are more narrow than the PDF for the plage.

Observed distribution functions for the vertical velocity in the range from −2 to 2 km s⁻¹ were shown by Franz & Schlichenmaier (2009). They found in the quiet Sun that upflows dominate over downflows for all velocity intervals. In the penumbra (they did not distinguish between the inner and outer penumbrae) upflows dominate for |v_z| < 400 m s⁻¹, downflows for faster velocities. The overall width of the PDF in the penumbra was about half of the value found in the quiet Sun. While our results certainly show about a factor of two difference in the widths (see also the vertical rms velocity in Figure 15), other details might also depend on the spectral line utilized as well as on resolution and likely require a comparison through forward modeling.

In Figure 21 we analyze subphotospheric processes responsible for the Evershed flow as well as the horizontal magnetic field found in filament flow channels. Similarly to Rempel (2011b) we focus here on regions with v_z > 0 and v_R > 0 in the inner penumbra, where the strongest driving of the Evershed flow takes place.

Zoom In Zoom Out Reset image size

The left panels ((a), (d), (g)) present averages of magnetic field, velocity, and temperature. Apart from the dramatic change of the Evershed flow speed we also find that the enhancement of B_R at τ = 1 is less pronounced at lower resolution, consistent with Figure 9.

The middle panels ((b), (e), (h)) show energy conversion terms that play a central role in the maintenance of the Evershed flow. As discussed in detail in Rempel (2011b) pressure buoyancy driving takes place in upflow regions, while downflows are close to hydrostatic, which is opposite to the situation in non-magnetic convection. The pressure driving in upflows is diverted by the magnetic field into the horizontal direction where the acceleration takes place, while the pressure force itself does not have a strong component in the horizontal direction. To quantify the influence from grid resolution on this process we show the work by the vertical component of the pressure/buoyancy force (black), work by the vertical component of the Lorentz force (red), work by the radial component of the Lorentz force (blue), and radial component of the acceleration force (green). Dashed lines correspond to simplified expressions that were discussed further in Rempel (2011b). The dotted line indicates the work by the radial component of the pressure force multiplied a factor of 10 for better visibility. Regardless of resolution we find the same balance, only the overall amplitudes of the various terms increase with increasing resolution. A quadratic increase with the peak Evershed flow velocity is expected, since the acceleration term is given in leading order by −〈ϱv_Rv_z∂_zv_R〉 and v_z itself does not change with resolution (tied to the fixed amount of energy transported).

The panels on the right ((c), (f), (i)) show the contributions from the stretching term (red), advection term (black), and divergence term (blue) in the induction equation for the radial field component. Here it is most striking that the stretching term and in particular the contributions from B_z∂_zv_R (red, dashed) triples in amplitude over the resolution range explored. As discussed in detail in Rempel (2011b), there is a strong connection between the Evershed flow and the horizontal magnetic field in a thin boundary layer beneath the τ = 1 surface. The subsurface shear profile of the Evershed flow induces strong horizontal field, which is distributed along the τ = 1 level by overturning convection and leads to a thin "outer shell" of the filament with strongly enhanced horizontal field. This boundary layer is also the place where strong magnetic tension forces are in balance with vertical pressure and horizontal acceleration forces (see also Figure 22). This strong linkage between B_R and v_R is also evident from Figures 5 and 9, in which both the Evershed flow and horizontal magnetic field show the strongest resolution dependence with a trend toward faster flows and stronger field at higher resolution.

Zoom In Zoom Out Reset image size

8.1. Role of Top Boundary Condition

8.2. Convergence with regard to Numerical Resolution

8.3. Comment on Diffusion

8.4. Visibility of Convective Signatures in Real Observations

In deepening our previous studies in Rempel et al. (2009b, 2009a) and Rempel (2011a, 2011b), we investigated how boundary conditions and grid resolution affect the structure of penumbrae in simulated sunspots and worked out observational consequences that can help to further constrain the magneto-convective structure of the penumbra. The main findings are the following.

• 1.  
  The magnetic top boundary condition appears to determine the radial extent of the penumbra. We do not get an extended penumbra if we impose a potential field subject to horizontal periodicity. Thus, simulations in periodic domains cannot predict the extent of penumbrae from first principles at this point.
• 2.  
  Magneto-convection forming penumbral fine structure and driving the Evershed flow is robust over the explored resolution range. While all explored grid resolutions reproduce basic aspects of a penumbra, a resolution of at least 48[24] km (horizontal [vertical]) is required for a quantitative agreement.
• 3.  
  The overall amount of convective motions in the photosphere as characterized by the vertical rms velocity shows little dependence on resolution and extent of the penumbra. It is linked to the intensity through an approximate relation of the form $I\propto \sqrt{v_{z}^{\rm rms}(\tau =1)}$.
• 4.  
  We find asymmetries between upflow/downflow filling factors and velocities that favor the detection of upflows in the inner and downflows in the outer penumbra. These asymmetries become more pronounced at higher resolution. Forward modeling likely overestimates the visibility of a convective flow pattern in penumbrae at this point.
• 5.  
  Evershed flow channels are strongly magnetized at photospheric levels, values around 1.5–2.5 kG are typical. The tendency for enhanced horizontal magnetic field within flow channels by a few 100 G compared to the background is a robust feature.
• 6.  
  Opposite polarity flux is present in the penumbra at a level of about 10% of the total flux of the sunspot, with no resolution dependence for an otherwise fixed setup. The overall amount of opposite polarity flux increases with the extent of a penumbra from 6% to 12%.
• 7.  
  Integrated over the entire penumbra, about 40% of the downward directed mass flux is associated with magnetic flux having a polarity opposite to that of the umbra. We see only a moderate increase of this number with resolution and extent of the penumbra. Thus, there should be a significant fraction of downflows in the penumbra with the same field polarity as the umbra, i.e., a substantial amount of mass flux returns by subduction of upward pointing magnetic field rather than by flowing along downward pointing field lines.
• 8.  
  While strong supersonic downflows are present in the outer penumbra (we found up to 15 km s⁻¹ flow amplitude), their overall contribution to the downward directed mass flux is, with 2%–3%, negligible.
• 9.  
  In the outer penumbra we find a substantial population of bright downflows that lead to a vanishing I–v_z correlation, while the correlation in the inner penumbra and surrounding plage region are comparable (reaching moderate values around 50%–60%). Nevertheless, energy is transported in all these regions by convection. Thus, low or moderate I–v_z correlations cannot be taken as an indication for the absence of convective energy transport, in particular not when magnetic field is present in the photosphere.
• 10.  
  Most of the vertical mass flux in the penumbra is balanced within each radial shell. The unbalanced component of the mass flux (upflows in the inner and downflows in the outer penumbra) is about 15% of the total overturning mass flux integrated over the whole penumbra.

We conclude this paper by summarizing the magnetic field and flow structure of penumbral filaments found in numerical simulations to date in Figure 22. This sketch is mostly based on the conclusions from this paper as well as Rempel (2011b; see Figure 17 therein).

The side view shows the central upflow regions of filaments where the average field and average flow are well aligned. The pressure driving in the upflows is almost orthogonal to the mostly horizontal acceleration of the Evershed flow, requiring the Lorentz force to facilitate the energy exchange (the net work by the Lorentz force is essentially zero). In the lateral downflow regions the magnetic field continues to point upward, leading to a substantial misalignment between flow and field and submergence of field lines. The side view describes mostly the situation in the inner penumbra, where we do not find a substantial amount of inverse polarity magnetic flux.

The top view shows the central upflow region with the adjacent lateral downflows. While mass is moving outward (Evershed flow) it moves to the edge of filaments where it returns beneath the solar surface mostly through submergence of field lines. Only near the outer end of the filament does a significant fraction of the mass flux return through flows along downward pointing field lines.

The cross section shows a structure with reduced magnetic field strength in the subphotospheric layers that is filled with an overturning convection pattern maintaining the penumbral brightness. While this structure has a lot of similarity with the idealized field free gap model of Spruit & Scharmer (2006) and Scharmer & Spruit (2006), it is not field-free. The filament has a core with a substantial vertical field component (see also Figure 17 of Rempel 2011b) that is spreading out laterally near the top of the filament. A fraction of it passes through the curved τ = 1 surface (indicated by a dashed line) and continues upward. Overturning convection can bend some of these field lines back downward, leading to opposite polarity flux along the side of filaments. If the τ = 1 level is located deep enough, this can lead to the appearance of opposite polarity flux along the edge of filaments in magnetograms. The Evershed flow (shown in green) increases rapidly toward the τ = 1 level in a region with a substantial vertical field component, leading to induction of a strong horizontal field component along the axis of the filament. Overturning motions distribute the magnetic field along the τ = 1 level in a thin boundary layer shown in red, resulting in strongly magnetized Evershed flow channels at photospheric levels that might appear like horizontal flux tubes to the observer. Note that we simplified the picture in the sketch with regard to the distribution of the horizontal field component. Horizontal field is present throughout most of the filament; it is however about a factor of two stronger in the boundary layer indicated in the sketch and mostly there dynamically relevant (the acceleration of the Evershed flow takes place in this boundary layer; see indicated force balance).

The author thanks M. Knölker for insightful discussions and very valuable comments on the manuscript. This research has greatly benefited from discussions that were held at the International Space Science Institute (ISSI) in Bern (Switzerland) in 2010 February as part of the international working group Extracting Information from Spectropolarimetric Observations: Comparison of Inversion Codes. The National Center for Atmospheric Research is sponsored by the National Science Foundation. Computing resources were provided by the National Institute for Computer Sciences (NICS) under grant TG-AST100005. We thank the staff at the supercomputing center for their technical support.

For initialization of our simulations we use an axisymmetric self-similar magnetic field configuration. This approach has been used previously for magnetohydrostatic sunspot models by Schlüter & Temesváry (1958), Deinzer (1965), and more recently by Schüssler & Rempel (2005):

$$\begin{eqnarray} B_z(R,z) &=& B_0\, f(\zeta)\, g(z)\,, \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} \vspace{5.69054pt}B_R(R,z) &=& {-}B_0\, {R\over 2}\,f(\zeta)\,g^{\prime }(z)\,. \end{eqnarray} \tag{ A2 }$$

Here, R, z correspond to the radial and vertical direction in cylindrical coordinates, g(z) is the vertical magnetic field profile along the symmetry axis, and $\zeta = R\sqrt{g(z)}$. The function f(ζ) describes the height-independent radial profile of the vertical field component and can be freely chosen. For initializing our simulations we used the following functions to describe the radial and vertical field profile:

$$\begin{eqnarray} f(\zeta) &=& \exp [-(\zeta /R_0)^4]\, \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} g(z)&=& \exp [-(z/z_0)^2]\,, \end{eqnarray} \tag{ A4 }$$

leading to a total of three parameters that determine the field structure: B₀, R₀, and z₀. The total flux is given by Φ₀ = 2.7833B₀R²₀. Note that the most important parameters for our problem are Φ₀ and B₀; the detailed choices of f(ζ) and g(z) are secondary. We have chosen g(z) in this form to ensure g'(z = 0) = 0 for consistency with our vertical field bottom boundary condition and f(ζ) to initially concentrate most field into the sunspot. The typically used Gaussian profile would lead to a smaller sunspot and a stronger plage region.

We insert the above magnetic field structure into a thermally relaxed hydrodynamic (or non-magnetic convection) simulation. We do not correct the thermal structure to account for magnetic forces, instead we allow this adjustment to happen as part of the dynamical evolution. We set however the entropy in regions with B > 2 kG to the value found at the bottom of the domain, while keeping the pressure unchanged and initially quench velocities in magnetized regions. The latter two have a rather minor effect on the overall evolution.

For the simulations presented in this publication we used a modified top boundary condition, which allows us to arbitrarily to change the inclination angle of the magnetic field, while maintaining horizontal periodicity. To explain this boundary condition we start from a potential field extrapolation in a horizontally periodic domain, which can be derived from a potential of the form (Cheung 2006)

$$\begin{equation} \Phi (x,y,z)=\int \tilde{\Phi }(k_x,k_y)e^{-\vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y\;, \end{equation} \tag{ B1 }$$

with $\vert k\vert =\sqrt{\vphantom{A^A}\smash{\hbox{${k_x^2+k_y^2}$}}}$. The magnetic field components are given by $\vec{B}=-\nabla \Phi$, together with the relation $\tilde{B}_0(k_x,k_y)=\vert k\vert \tilde{\Phi }(k_x,k_y)$. Here $\tilde{B}_0$ is the Fourier transform of B_z(x, y, 0) at the boundary z = 0:

$$\begin{eqnarray} B_x &=&{-}\int \tilde{B}_0(k_x,k_y) \frac{i k_x}{\vert k\vert } e^{-\vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y\nonumber \\ B_y &=&{-}\int \tilde{B}_0(k_x,k_y) \frac{i k_y}{\vert k\vert } e^{-\vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y\nonumber \\ B_z &=&{+}\int \tilde{B}_0(k_x,k_y) e^{-\vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y. \end{eqnarray} \tag{ B2 }$$

From this potential field extrapolation we derive a generalized non-potential field boundary condition by introducing an additional parameter α controlling the field inclination angle:

$$\begin{eqnarray} B_x^{\alpha } &=&{-}\alpha \int \tilde{B}_0(k_x,k_y) \frac{i k_x}{\vert k\vert } e^{-\alpha \vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y\nonumber \\ B_y^{\alpha } &=&{-}\alpha \int \tilde{B}_0(k_x,k_y) \frac{i k_y}{\vert k\vert } e^{-\alpha \vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y\nonumber \\ B_z^{\alpha } &=&{+}\int \tilde{B}_0(k_x,k_y) e^{-\alpha \vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y. \end{eqnarray} \tag{ B3 }$$

While using a value of α = 1 reproduces the potential field, a value of 0 ⩽ α < 1 leads to a field more vertical than potential (vertical field for α = 0), a value of α > 1 to a field extrapolation more horizontal than potential.

The electrical current distribution associated with this field extrapolation is given by

$$\begin{eqnarray} j_x^{\alpha } &=&{+}(1-\alpha ^2) \int \tilde{B}_0(k_x,k_y) i k_y e^{-\alpha \vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y\nonumber \\ j_y^{\alpha } &=&{-}(1-\alpha ^2) \int \tilde{B}_0(k_x,k_y) i k_x e^{-\alpha \vert k\vert z} e^{i k_x x}e^{i k_y y} dk_x\, dk_y\nonumber \\ j_z^{\alpha } &=& 0. \end{eqnarray} \tag{ B4 }$$

While the vertical current remains zero in the region z > 0, closed horizontal current loops induce magnetic field which changes the effective field inclination in z ⩽ 0. The magnetic field in the region z > 0 is not force-free, but the magnetic field in the computational domain z < 0 remains close to force-free due to the low β conditions we consider.