Issue 
A&A
Volume 551, March 2013



Article Number  A92  
Number of page(s)  12  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201117361  
Published online  28 February 2013 
Structures in compressible magnetoconvection and the nature of umbral dots
^{1}
Department of Physics, Gustaf Hällströmin katu 2a PO Box 64,
00014
University of Helsinki,
Finland
email:
chunlin.tian2012@gmail.com
^{2}
Eötvös University, Department of Astronomy,
Budapest, Pf. 32, 1518
Hungary
Received:
27
May
2011
Accepted:
19
December
2012
Context. Structures seen in idealized numerical experiments on compressible magnetoconvection in an imposed strong vertical magnetic field show important differences from those detected in observations or realistic numerical simulations of sunspot umbrae.
Aims. To elucidate the origin of these discrepancies, we present a series of idealized 3D compressible magnetoconvection experiments that differ from previous such experiments in several details, bringing them closer to realistic solar conditions.
Methods. An initially vertical magnetic field B_{0} is imposed on a time snapshot of fully developed solarlike turbulent convection in a layer bounded by a stable layer from above. Upon relaxation to a statistically steady state, the structure of the flow field and magnetic field is examined.
Results. Instead of the vigorous granular convection (GRC) well known to take place in magnetized or weakly magnetized convection, for high values of B_{0} heat is transported by smallscale convection (SSC) in the form of narrow, persistent convective columns consisting of slender upflows accompanied by adjacent downflow patches, which are reminiscent of the “convectons” identified in earlier semianalytic models. For moderate field strengths, flux separation (FXS) is observed: isolated fieldfree inclusions of GRC are embedded in a strongly magnetized plasma with SSC. Between the SSC and FXS regimes, a transitional regime (F/S) is identified where convectons dynamically evolve into multiply segmented granular inclusions and back.
Conclusions. Our results agree in some aspects more closely with observed umbral structures than earlier idealized models, because they do reproduce the strong localized, patchy downflows immediately adjacent to the narrow convective columns. Based on recent observations of umbral dots, we suggest that in some cases the conditions in sunspot umbræ correspond to the newly identified F/S transitional regime.
Key words: magnetohydrodynamics (MHD) / sunspots / convection
© ESO, 2013
1. Introduction
The spectacular increase in the resolution of solar observations experienced in the past decade has resulted in unprecedentedly detailed images and movies of the smallscale structures seen in sunspots. This in turn has rekindled interest in theoretical studies of magnetoconvection in order to understand the origin of these structures. These studies have resulted in some spectacular breakthroughs, such as the first ever numerical simulation of whole sunspots (Rempel et al. 2009a,b; Rempel 2012). Despite these successes, however, there are still many questions to be answered. They include the exact nature of umbral substructures, such as umbral dots, light bridges and dark nuclei, and their relation to structures seen in different numerical magnetoconvection experiments.
High resolution imaging and spectropolarimetric observations with the Swedish Solar Tower (Sobotka & Hanslmeier 2005; Riethmüller et al. 2008b; Sobotka & Puschmann 2009; Ortiz et al. 2010), Hinode (Kitai et al. 2007; Riethmüller et al. 2008a; Bharti et al. 2007b, 2009; Sobotka & Jurčák 2009) and the Dunn Solar Telescope (Bharti et al. 2007a; Rimmele 2008) have shown that umbral dots (UDs) have characteristic sizes of ~200 km and that they are mostly shortlived, the majority having lifetimes of a few times ten minutes or less. Larger UDs often split or coalesce, and in high resolution images, dark lanes are often seen inside them. In the bright parts of UDs, an upflow of ~1 km s^{1} has been unambiguously detected in many instances. Strong localized downflow patches have been detected on the periphery of UDs. This association with up and downflows clearly demonstrates that UDs are a manifestation of magnetoconvection in the strongly magnetized, compressible umbral plasma. Kilcik et al. (2012) have studied the morphological and statistical properties of bright UDs by comparing the high resolution data from New Solar Telescope at the Big Bear Solar Observatory and 3D magnetohydrodynamic (MHD) simulations. Based on good seeing, Watanabe et al. (2012) studied the temporal evolution of different kinds of UDs in detail. These very recent studies show not only agreement but also discrepancies among observations made by different authors and numerical simulations, such as lifetime, size, shape, downflows, and correlations between different quantities. We note that these studies looked at different sunspots, and the simulated sunspots are also different from the observed ones. So far, no clear picture has emerged either from observations or simulation study regarding how the detailed properties of UDs depend on the properties of the sunspots (i.e. field strength, overall size, etc.). In this sense it is not clear to what degree discrepancies are simply due to the natural variation of sunspot properties and to what degree there are true discrepancies between realistic simulations and observed sunspots in terms of a model deficiency.
This paper is organized as follows. Section 2 gives the context and objectives of this paper. Section 3 describes the setup of our numerical model. Section 4 presents and discusses the results from the simulations, while Sect. 5 concludes the paper.
2. Context and objectives
2.1. Idealized models of magnetoconvection with an imposed vertical field
In the past decades, theoretical studies of magnetoconvection in a strong, vertically oriented magnetic field have gone a long way, starting from nonlinear extensions of the classic linear stability analysis of Chandrasekhar (1961). These studies have proceeded through the modelling of convecting magnetized Boussinesq and anelastic fluids in two and three dimensions to modelling 3D compressible magnetoconvection in a polytropically stratified atmosphere (Weiss et al. 1996, 2002). One important realization in these studies was that the spatial separation, a mechanism that separates strongly magnetized nonconvecting plasma from weakly magnetized convective flows, allows for overturning convection to also take place in situations where linear theory would predict convective stability or overstability.
As the field strength is decreased from a value high enough to completely suppress motions, convection first sets in in the form of narrow, needlelike vertical convective columns. The horizontal scale of these structures is smaller than the natural scale of unmagnetized convection in the same layer, so this stage is usually referred to as “smallscale convection” (SSC). While these convective columns can even exist as solitary “convecton” solutions (Blanchflower 1999; Blanchflower & Weiss 2002; Houghton & Bushby 2011), they typically appear in large numbers. The “convectons” of earlier Boussinesq models (Blanchflower 1999; Blanchflower & Weiss 2002) are vertically elongated convective cells, with upflows and adjacent downflows embedded in stagnant fluid. In later compressible experiments (Weiss et al. 1996; Houghton & Bushby 2011), however, SSC takes the form of isolated, convective upflow plumes, their mass flux being compensated for by a slow, uniform sinking motion in the surrounding magnetized fluid. The plumes are persistent and immobile, though they pulsate irregularly for a weaker imposed field .
With further reduction of the imposed magnetic field, a phenomenon known as flux separation (FXS) takes place: large isolated patches of unmagnetized, vigorously convecting fluid appear inside the strongly magnetized component where SSC is still going on. The character of the vigorous convection in fieldfree patches is similar to the normal granular convection (GRC) in the unmagnetic case. Flux separation was first detected in a simulation by Tao et al. (1998); however, it seems to be a more fragile mode of magnetoconvective energy transport than the SSC as in simulations without full compressibility or with lower aspect ratios it is not present. Also, magnetoconvective systems show a hysteresislike behaviour where the presence or absence of FXS is not uniquely determined by the physical parameters of the system but also depend on its history (Weiss et al. 2002).
For even lower field strengths, the strongly magnetized component becomes confined to an intermittent network of channels separating the vigorously convecting patches. Finally, with the further decrease in the imposed field the topology changes, giving way to the well known pattern of isolated magnetic flux concentrations arranged in a magnetic networklike pattern between granules (GRC).
The sequence described above, characterized as a competition between two modes of convection (columnal SSC vs. vigorous GRC), suggests an interpretation of UDs as the photospheric manifestation of SSC. At the same time, there are no obvious examples of FXS in the solar photosphere, and the short lifetimes of observed UDs stand in stark contrast to the persistent character of the SSC columns described above.
2.2. Numerical simulations of magnetoconvection in sunspot umbrae
Models like those described above are usually called “idealized” because they aim at a systematic study of increasingly general magnetoconvective systems that are still simple enough to allow parameter studies and to be understood in relatively simple physical terms. In consequence, they do not aspire to faithfully reproducing all the physical conditions relevant to the solar photosphere; instead, they consider convection in a polytropically stratified perfect compressible fluid, with closed boundaries and fixed diffusivities.
Another approach to the problem consists in attempting to simulate conditions prevailing in the solar photosphere as realistically as possible, including radiative transfer and a realistic equation of state in a fully compressible equilibrium model (Moradi et al. 2010). An obvious drawback to this approach is that, when going to the very limits of the available computing power, it is left with too few resources for systematic parametric studies. Nevertheless these ambitious models have been quite successful in e.g. reproducing the observed penumbral fine structure. Structures readily identified with the observed umbral dots are also seen in such simulations, as first demonstrated by Schüssler & Vögler (2006). In contrast, with the SSC columns of idealized models, these simulated UDs have finite lifetimes, and they show a characteristic “coffee bean” structure, being crossed by a dark lane along the long axis of their oval shape. Downflows are present near the end points of the dark lane. The dark lane itself is due to an optical depth effect due to the pileup of upflowing material below a magnetic cusp above the UD. Similar UDlike structures are also seen in the fullspot simulations of Rempel et al. (2009a, 2009b). Observationally, similar structures were found in the cases studied by Bharti et al. (2007b), Rimmele (2008), and Ortiz et al. (2010).
2.3. Objectives
Despite the impressive achievements of these stateoftheart simulations, some differences between the detailed properties of observed and simulated umbral dots remain (cf. Bharti et al. 2010). Even more important differences persist, however, between the nature of magnetoconvective structures in idealized models, on the one hand, and in observed or simulated sunspots, on the other. These differences include the contrast between the observed localized patchy downflows next to UDs vs. the lack thereof around SSC columns in experiments and the persistent character of convective columns vs. the limited lifetimes of UDs.
We cannot claim a satisfactory understanding of magnetoconvection in sunspots until the origins of these differences are properly elucidated and a convincing link between convective columns and UDs is forged. For this, there is a need to bridge the gap between the idealized models and the realistic simulations. As a first step in this process, in this paper we constructed an idealized magnetoconvection model where, however, we relaxed certain constraints employed in earlier idealized models such as Weiss et al. (2002) to bring the model closer to solar conditions. These were:

a twolayer model wherein the convectively unstable fluid isbounded by a stable layer from above;

open lower boundary conditions;

a subgrid closure scheme to account for the effects of smallscale turbulent transport;

overall stratification, while polytropic for simplicity, to be reasonably similar to the actual solar case;

a novel numerical scheme, the BhatnagarGrossKrookMHD (BGKMHD) scheme that is based on gaskinetic theory. We implemented a numerical solver for the induction equation in addition to the BGK NavierStokes equations (BGKNS) solver developed by Tian et al. (2007). The BGK gaskinetic scheme is particularly well suited to the capture of neardiscontinuities in the plasma flows. This admits a satisfactory stabilization of transition from initial perturbation to fully developed turbulent convection.
In the simulations we first relax a convectively unstable layer of perfect gas polytropically stratified under gravity. Then, a vertical uniform magnetic field is imposed on a snapshot of the fully developed, statistically steady turbulent convection. After a statistically steady state of magnetoconvection is attained, we interpolate the results on a higher resolution mesh and resume the calculations. Upon reaching steady state again, the properties of this system are examined. We construct a series of models with different strengths of the imposed magnetic field.
3. Numerical model
In the current study, the 3D BGKNS solver developed by Tian et al. (2007) is coupled with an induction equation solver and applied to solving the following resistive MHD equations under a gravitational field:
where p_{tot} = p_{g} + p_{m}, E = E_{i} + E_{m} + E_{k}, p_{g} is the gas pressure, is the magnetic pressure, E_{i} is the internal energy, is the magnetic energy, and is the kinetic energy. Here, g is the gravitational acceleration, F_{d} the diffusive heat flux, Σ the viscous stress tensor, and η the magnetic resistivity. All the other symbols have their standard meanings.
3.1. Boundary conditions
Side boundaries are periodic. The upper boundary is impenetrable, and the density and internal energy are fixed: (5)The lower boundary is transmitting, so here only the horizontal averages of density and internal energy are fixed at their initial values (a slightly superadiabatic state), while the vertical derivatives of their fluctuations are set to zero: (6)Throughout this paper, overbar and prime denote the average and fluctuating parts, i.e. for a variable a we have . Unless explicitly stated otherwise, overline denotes a horizontal average.
The magnetic field is constrained to become vertical, (7)at the top boundary. The lower boundary is transmitting for magnetic field (8)In both cases, we have ∇·B = 0 at the boundaries.
3.2. Initial setup
Initially, the system consists of two layers: an upper stable radiative layer (1 > z > z_{1}) and a lower unstable convective layer (z_{1} > z > 0). The upper layer is subadiabatically stratified, the lower layer is slightly superadiabatic:
where Z_{1} and Z_{2} are the dimensionless temperature gradients of the subadiabatic and superadiabatic layers, respectively; m_{1} and m_{2} are the polytropic indices for radiative and convective layer, respectively; T_{0}, p_{0} are the temperature and pressure at the top lid; T_{1}, p_{1} are the temperature and pressure at the interface between the initial stable and unstable layers, i.e., In hydrostatic equilibrium we have (13)The choice of Z_{1} and Z_{2} determines the temperature contrast between top and bottom, by fitting to a standard stellar model, the vertical extent of our computational domain. For this reason Z is occasionally called a depth parameter in the literature. For m_{2} we choose a value below the adiabatic value that reproduces the right value of gravity acceleration according to Eq. (13). For Z_{1} a lower value is adopted to make the top layer subadiabatic, thus convectively stable. Then m_{1} can also be determined from Eq. (13). The specific values of these parameters are listed in Table 2.
3.3. Treatment of radiation and subgrid scale motions
Viscosity and diffusivity due to subgrid scale turbulent eddies are calculated according to the Smagorinsky model. The dynamical viscosity is given by (14)where c_{μ} is an adjustable constant, usually chosen from the range 0.1–0.2. The filter width Δ is taken to be the local resolution, a colon stands for tensor contract, and σ = ∂_{i}ν_{j} + ∂_{j}ν_{i}. The turbulent heat transfer coefficient C_{S} is calculated as (15)where the Prandtl number δ is assumed to be constant.
For the current study, radiative transfer is considered in the diffusion approximation. The total (radiative+turbulent) diffusive heat flux is then (16)where S = C_{p}(lnT − ∇_{a}lnp) is the specific entropy, C_{p} the specific heat at constant pressure, and ∇_{a} the adiabatic gradient. In the stable layer, C_{T} is set so that radiation carries out the input energy flux. In the convection zone, C_{T} is very close to zero. As C_{S} = μ/δ represents turbulent diffusion, it is set to zero in the stable region.
Dimensional analysis.
Summary of numerical parameters.
3.4. Dimensional analysis and model parameters
We nondimensionalize all variables by choosing the following basic scales. The scale of the length, l_{scl}, is the depth of the computational domain. For the density, gas pressure, and temperature, their respective values at the top boundary are chosen as the basic scales, i.e., ρ_{scl}, p_{scl}, and T_{scl}. The scales for the rest of the quantities can be deduced according to dimensional analysis. For instance, the scale of velocity and time are (17)respectively. (Note that ν_{scl} is the isothermal sound speed at the top.) In the case of magnetic field strength, a different definition is used in order to give it a more physical meaning (of thermal equipartition field strength at the top of the convective layer): , p_{∗} is the pressure at the top of the convection zone.
Imposed magnetic field strength in the numerical runs.
Of particular importance is the scale of the input energy flux, (18)where ν_{∗} is the sound speed at the convective top. As the solar radiation flux is fixed in physical units, its value in nondimensional units will depend on the value of F_{scl}, i.e. on the position of the upper boundary of our convection zone inside the solar convective envelope. Thus, by specifying the input flux, we can approximately control the radial location of our numerical computational domain within the standard solar model.
In summary, the input energy flux, the aspect ratio and the depth of the domain determine the “geographic” parameters of the model, i.e., size and location. The scales of density, gas pressure, and temperature determine the state of matter. To perform the dimensional analysis, we need a standard solar model. In the current study, we adopt the combined model calculated by Guenther et al. (1992, Table 3B therein).
In our thermally relaxed unmagnetized convection model, the extent of the top conductive layer (z > 0.9) is ~1.87PSHs (pressure scale heights) and the fully convective zone (z < 0.78) ~3.61 PSHs. Between them, there is a transition layer of around ~1.26 PSHs, including an upwardly overshooting zone. The solartype granulation forms at around z = 0.8, the lower part of the transition zone. We define z = 0.8 as the effective top boundary of the convection zone. Thus the upper conductive zone can be regarded as an artificial layer. The total energy flux transported by the system is ~0.1, implying that the top of our convective layer is at a depth of ~520 km from the solar surface (from unit continuum optical depth). This distance roughly corresponds to the depth of the Wilson depression. In such a case, the top of our box is already above the solar photosphere, approaching the height of the temperature minimum.
The scales of our models are summarized in Table 1. Owing to the rapid change of pressure near the solar surface, different standard solar models yield quite different surface pressures, so our current dimensional analysis is preliminary, especially the magnetic strength, which is determined directly by the unit of gas pressure.
For magnetoconvection it is in fact more meaningful to express magnetic field strength in terms of the turbulent equipartition field rather than in Alfvénic units (or thermal equipartition field). The turbulent equipartition magnetic field B_{e} is defined by . We evaluate B_{e} at z = 0.8, the effective top of convection zone. We find B_{e} = 0.51B_{scl} = 1484.89 G. Table 3 presents the strength of the imposed vertical magnetic field in the individual runs in units of B_{e}.
Fig. 1 Statistical properties of steady unmagnetized convection, used as initial condition for magnetoconvection. Shown are the averaged energy fluxes, in units of 6.34 × 10^{11}erg/s/cm^{2}. Averages are taken in time and in the horizontal plane. z is the height from bottom, in units of 1.72 Mm. 
4. Numerical results and discussion
Upon starting the calculation with no magnetic field and after long term relaxation, convection approaches a nearly steady state where the temporally and horizontally averaged total energy flux is nearly constant. The various contributions to the energy flux are shown in Fig. 1. The strong downwardly directed kinetic energy flux is a well known hallmark of steady compressible convection, as expected from convection theory (Unno & Kondo 1989; Petrovay 1990) and as is also consistent with numerical simulations (Chan & Sofia 1986). From the constancy of the total energy flux, it is clear that the unmagnetized convection is completely relaxed.
Next, we superpose a vertical uniform magnetic field on an instantaneous snapshot of the relaxed unmagnetized convection. Figure 2 shows the time history of magnetic energy for the cases listed in Table 3. It is apparent that after a short while, the magnetic energy contained in the system becomes statistically steady. For the cases with weak initial fields, the magnetic fields are amplified by flow motions and eventually reach a nearly steady value. For the strong initial field cases, the magnetic energy grows up to a maximum and then decreases to a nearly constant value.
4.1. Flow morphology
The effects of a vertical magnetic field imposed on the granulation pattern are summarized in Figs. 3 and 4. As revealed by previous studies, e.g., Weiss et al. (1996) and Stein & Nordlund (2000), magnetic fields are intermittent for low field strengths, concentrated into flux tubes confined to the narrow lanes of the intergranular network (case H). The strongest concentration occurs at the junctions of the intergranular lanes of larger granules in the deeper layer, forming porelike structures. This is because in such cases, the strong magnetic fields are generated by amplification of flow motions. The downstream in the intergranular network of deeper bigger granules has a very long extent and persists for a long time, so it can continuously amplify the magnetic fields.
In the opposite limit of a very strong imposed magnetic field, convective flows are significantly suppressed (cases A and B). The remaining flows and magnetic fields show a brainpattern configuration. At first glance, it looks like a numerical effect. Since this kind of structure is quite steady, we cut a piece of the computational domain and refined the grids, then resumed the calculations. We repeated this procedure until the horizontal resolution was improved by a factor of 8 in each direction. Numerical tests show that the brainpattern is an unresolved weak convection pattern on a scale even smaller than the size of SSC. After doubling the grid points per unit, this kind of structure is already resolved. Their shape persists when the resolution is increased further. These high resolution results indicate they are indeed a very weak, nonefficient, laminar, steady tiny scale convection pattern. The fine unresolved steady convection is evident in the vertical cuts in the upper panels of Figs. 5 and 6. We expect that implementing more realistic radiation transfer and boundary conditions would suppress this fine structure, making it undetectable in intensity maps and dopplergrams.
Fig. 2 Histories of magnetic energy for various runs, as indicated by the labels. 
Fig. 3 Dopplergrams at various depths for selected cases. Corresponding magnetic fields are shown in Fig. 4. Velocities are given in units of 6.03 km s^{1}. 
Fig. 4 Contour plots of vertical magnetic fields at various depth for selected cases. Corresponding vertical flow fields are shown in Fig. 3. Magnetic fields are given in units of 2911.55 G. 
Case C acts like a critical situation, where a few small bright convective cells are embedded in the brainpattern background. Their size is very close to the brainpattern. The number and size of such bright dots increase when reducing the imposed magnetic fields, cases D and E show. The size of these SSC elements is smaller than normal granules, inviting comparison to the umbral dots observed in sunspots. This correspondence was indeed suggested by Weiss et al. (2002); however, the SSC elements in their experiments only contained upflows, their upwards directed mass flux being compensated for by a slow sinking motion in the strongly magnetized bulk fluid. In contrast, we find strong localized downflows adjacent to the upflows, within the unmagnetized convective columns. This property makes our SSC elements a close relative of the convectons of Blanchflower (1999) and Blanchflower & Weiss (2002). So for brevity, in the following we refer to them by this name.
In cases F and G we recover the typical FXS phenomenon, first discovered by Tao et al. (1998), where SSC and vigorous GRC coexist, with the latter present in large field free patches (in what follows: inclusions, for brevity) embedded in a strongly magnetized fluid displaying SSC in the form of scattered convectons. Each inclusion includes several granules. The size of the inclusions is bigger than a normal granule and smaller than the granules in the deep convective zone.
Finally, in some regions of the snapshot for case E in Fig. 3 displays multiply segmented convective plumes of a size midway between inclusions and convectons. The lanes segmenting these elements are cool and coincide with downflows. This already suggests that our case E represents a heretofore unknown transitional regime of magnetoconvection between SSC and fully developed FXS, a finding that will be corroborated below.
Fig. 5 Vertical cuts of the computational box along x = 3 at an arbitrary instant. The colour scale represents temperature fluctuation from its horizontal mean, i.e., , in units of 5063.52 K. Rows correspond to different runs as labeled on the left. 
Fig. 7 Root mean square (rms) velocity and magnetic fields from their horizontal means for case F. z is the height from bottom, in units of 1.72 Mm. Units for velocity and magnetic field strength are 6.02 km s^{1} and 2911.55 G, respectively. 
Figures 5 and 6 present vertical cuts of our box. As expected, a stronger magnetic field is generally associated with cooler gas, usually even cooler than in the less magnetized downflows. In the weak field cases, the downflows show typical von Kármán vortices. The magnetic field is amplified and confined by the downward turbulent flows. In the strong field case, on the other hand, the turbulent character of the flow is significantly suppressed, field lines remain close to vertical, and the columnar shape of convectons is apparent. From Fig. 5, we can see that the maximum of temperature fluctuation from its horizontal mean occurs near the convective top. Although for stronger fields, the fluctuation is less turbulent, its amplitude is higher.
Both vertical cuts (Figs. 5 and 6) and the rms of flow and magnetic fields (Fig. 7) show that the most efficient place where magnetic fields are amplified by flow motions is near the top of the convective zone. The enhancement of magnetic field is not completely coherent with flow motions. The maximum fluctuation of magnetic fields is higher than the velocity. This kind of enhancement is significant. For case A, the imposed initial fields are 3B_{e} ~ 4455 G, and the amplified strong fields can be 2B_{scl} ~ 6000 G. In the intermittent case, H, occasionally, the magnetic fields confined in pores can also be very strong (bottom panel in Fig. 6).
4.2. Time evolution
Figure 8 illustrates the time evolution of vertical velocity near the top of the convective layer. The evolution is more clearly seen in the corresponding animations^{1}.
Cases C and D represent SSC. The evolution of convectons in these cases is relatively blurred. Case C seems to be a critical situation where brightening of some unresolved tiny convection cells can be seen occasionally. In case D, the convectons are clearer than in case C. Most of them are steady, evolving slowly.
The time variation of the convective elements in case E reveals a new type of behaviour apparently not seen in previous simulations. During sizeincrease phase of some convectons, one or multiple downdrafts appear inside them; these downdrafts merge into a network of internal downflow lanes, giving the convective elements a segmented appearance reminiscent of granular inclusions. (Again, note the analogy with exploding granules, cf. Hirzberger et al. 1999.) The inclusions have a lifetime of ~70 to ~160 min, depending on their size. The single UD between two inclusions lasts around 30 min. What we see here is essentially convectons developing into inclusions and back, so it represents a hitherto unseen transitional state between SSC and FXS.
In the FXS cases, F and G, we have normal GRC going on inside the inclusions while small, relatively steady convectons are occurring in the magnetized fluid in between. These two kinds of structures evolve separately at different locations. There seems to be no more evolution of convectons to or from inclusions. Although there are smaller granules growing and decaying near the edge of inclusions, they should not be regarded as convectons.
Compared with granules, porelike structures in Case H have a much longer lifetime. In fact, in our simulations, the pores persist ever since they formed. They are just stretched and moved by convective flows from one place to another. The reason might be that pores are confined by big granules in the deep convection zone and they have longer lifetimes than near surface granules. Currently we do not have enough computational resources to let our models run long enough to allow obvious deformation and reformation of these deep big granules.
4.3. Flow statistics
For more quantitative analysis of the properties of turbulent magnetoconvection we plotted the PDFs (probability distribution functions) of the vertical components of the magnetic field strength (in Alfvénic units) and of the flow velocity in Figs. 9 and 10, respectively. Only typical cases were plotted, i.e., completely suppressed (case A), SSC (case C), FXS (case F), and GRC with intermittent fields (case H). In Fig. 10 the unmagnetized convection case is also shown. The data were sampled at two horizontal layers: near the upper convective boundary (z ~ 0.8) and deep inside the convective part (z ~ 0.4).
The shape of PDFs for the distribution of vertical Alfvénic speed is different in the two layers sampled. The width of PDFs is broader and sharper in the top of the convection zone than in the middle, where the magnetic fields are more intermittent. Our PDFs differ from those obtained by Weiss et al. (2002) in the strong field region. We do not have a secondary hump near the maximum. This might be because our models have transmitting lower boundary and some of the strong magnetic fields can be transported out of the box freely, and thus there is no concentration of magnetic fields near the maximum.
When the turbulent convection is completely suppressed, the minimum of Alfvénic speed is far from zero, and its PDF has a narrow width. In such a case, the flow can affect the magnetic fields but cannot bend them downwards.
Fig. 8 Time evolution of vertical velocity at z = 0.8, near the top of the convective layer. The time intervals between snapshots are ~71 min for cases D, E, F, G; ~86 min for case H. Velocities are given in units of 6.03 km s^{1}. 
For SSC (case C), the PDFs of vertical Alfvénic speed are similar to case A, but with a wider width. A minimum of the magnetic fields are approaching zero in both layers. As we mentioned before, this is highly likely to be a transition between SSC and a completely suppressed convection case. An SSC cell is nearly indistinguishable in this case, so we may conclude that when the imposed vertical field lines are so strong that their orientation cannot be turned downwards, the convection starts to be completely suppressed.
For the intermittent field case with GRC (case H), the forms of the PDF of vertical Alfvénic speed are similar at different depths. Both have a symmetric peak around the origin and a broad slow slope at higher strengths. The symmetric peak corresponds to the weakly magnetized, vigorously convecting fluid, while the slope is due to the magnetic flux concentrations. This form of PDF is familiar both from previous simulations and from quiet Sun observations (de Wijn et al. 2009).
For FXS (case F), the shape of the PDF of vertical Alfvénic speed looks like a combination of GRC and SSC, especially near the convective top. It has a sharp symmetric peak, a broad flat region, and a quasi cutoff near the maximum.
The PDFs of the flow velocity sampled in different layers show more obvious discrepancies (Fig. 10). For GRC with intermittent magnetic field (case H), the PDFs of vertical velocity are almost identical to those of unmagnetized convection. Their change in shape between the top and the middle of the convecting layer is related to the well known change in morphology from the observed granular pattern near the surface to isolated fast downdrafts embedded in a slowly rising, warm bulk fluid in the bulk of the convective zone. The lack of a peak at negative velocities in the top panel may be due to the low filling factor of the downdrafts and possibly also to a lack of a characteristic scale among them.
As the imposed magnetic field is increased, the main peak of the PDFs of the vertical velocity shifts towards zero. This indicates a general slow sinking of the magnetized fluid component. An extended shoulder on this peak on the negative side corresponds to the faster downdrafts localized next to the upflows. Both this shoulder and the one on the positive side representing the upflows are more distinct near the top of the convecting layer.
The anisotropy of turbulence is shown in Fig. 11. As expected, the ratio of vertical component to total field strength approaches 1 as we increase the imposed magnetic fields (bottom panel). For the weakest field case in our study, this ratio is around 0.78 near the bottom. The merit of implementing the transmitting lower boundary condition is evident. The strong anisotropy also occurs near the top of transition region, at z = 0.9.
In the upper panel, the anisotropy of velocity for unmagnetized convection is also plotted. It is clear that the magnetic fields suppress the flow motions. For the strong magnetic fields, the velocity field is completely vertical.
4.4. Comparison with observations
In the weakly magnetic cases, the morphology, time dependence, and statistics of our magnetoconvective flow are in good accordance with observations. This is so despite the uncertainty of our dimensional analysis, which makes it inadequate to directly compare the size and lifetime of all kinds of structures with observations. We attribute this difference mainly to the lack of a freely radiating photosphere in our experiments, which is known to play a key role in determining granular scales (Rast 2003). The overall qualitative and quantitative similarity of the magnetoconvective structures in our weak field models to observations encourages us to look for similar observational parallels with the structures seen in our strong field experiments.
The analysis of high resolution Doppergrams and continuum images of sunspot umbrae by Bharti et al. (2007a) show that the umbral dots are often situated adjacent to downwflows. Localized downflows near UDs have also been observed by Ortiz et al. (2010). These observed umbral structures are quite similar in appearance to the convectons surrounded by a patchy ring of downflows, as seen in our experiments in the SSC domain.
Furthermore, dark lanes are often observed in larger umbral dots. These UDs rarely have the regular “coffee bean” structure suggested by the simulations of Rempel et al. (2009b); instead, dark lanes are often asymmetrically positioned inside them, and multiple dark lanes are common, lending a segmented appearance to these UDs. Umbral dots with up to six dark lanes have been observed by Bharti et al. (2007b). Observations also indicate that the UDs evolve rapidly compared to sunspots. Bharti et al. (2007b) show how a UD with multiple dark lanes changed its appearance over the course of about an hour. (However, some of the apparent temporal variations in UDs may be due to seeing effects, cf. Hamedivafa 2011.) In the study of Schüssler & Vögler (2006), small UDs primarily show a more symmetric “bean shape”, while in their Fig. 1 there are also a few examples of UDs with more complicated dark lanes. The size of the UD typically depends on how complicated the substructure is.
These properties of UDs are more reminiscent of the transitional convective structures seen in our case E. Indeed, case E represents a vertical magnetic field strength of B_{0} ~ 2.2 kG in physical units, considering the uncertainty of dimensional analysis, it is just in the range expected for sunspot umbrae. We therefore suggest that conditions in sunspot umbræ correspond to the FXSSSC (F/S for short) transitional regime identified in this paper.
Fig. 9 PDFs for the distribution of vertical Alfvénic speed, , at z = 0.8 and z = 0.4 for various imposed vertical magnetic field strengths as indicated inside the labels, corresponding to, totally suppressed (black solid), SSC (blue dotted), FXS (magenta dashed) and GRC (red dash dot). 
Fig. 10 PDFs for the distribution of vertical velocity, i.e.,v_{z}, at z = 0.8 and z = 0.4 for various imposed vertical magnetic field strengths as indicated inside the labels: totally suppressed field (thick black solid), SSC (blue dotted), FXS (magenta dashed), and GRC (red dash dot). The blue solid thin is the PDF of the unmagnetized case. 
Fig. 11 Vertical anisotropy of the velocity field (upper panel) and magnetic field (in Alfvénic units, lower panel), for cases A (thin black solid), B (green dotted), D (cyan dashed), E (blue dashdot), F (purple dash–3dot), G (magenta long dashes), and H (red solid). z is the height from bottom, and ⟨ ... ⟩ represents horizontal average. 
5. Conclusion
We have constructed a series of numerical experiments of magnetoconvection in a layer of compressible plasma with an open lower boundary and bounded by a stable layer from above, using a BGKMHD scheme. Our model setup differs from that of previous idealized magnetoconvection models (Weiss et al. 2002) in several important aspects by being closer to realistic solar conditions. Despite this, the general behaviour of the solutions as a function of the imposed field strength basically agrees with those models. In particular, FXS still prevails in a wide parameter range, and the SSC dominates the energy transport in the strong field case.
An important difference compared to the results of Weiss et al. (2002) is that in the SSC regime strong, narrow, patchy downflows appear in a ring around the narrow upflows in our experiments. Thus, our “convective columns” may be more closely related to the solitary “convectons” of Blanchflower & Weiss (2002) than to the localized plumes in previous magnetoconvection experiments. This property of our experiments agrees with recent high resolution observations of sunspot umbrae (Bharti et al. 2007a; Bharti et al. 2009; Ortiz et al. 2010).
Furthermore, we have identified a hitherto unknown transitional regime between FXS and SSC. This F/S regime is very dynamic, with a continuous evolution of convectons into granular inclusion and back. Intermediate states of the convective elements are often reminiscent of larger, multisegmented umbral dots, also known to be subject to dynamic changes. This suggests that conditions in sunspot umbræ correspond to the F/S transitional regime identified in this paper and that UDs essentially represent convectons and/or convecton/inclusion transitional phases as seen in our run E.
It should be noted that the realistic numerical simulations of Schüssler & Vögler (2006) and Rempel et al. (2009b) suggest a somewhat different interpretation of the substructure and time dependence seen in UDs. The dark lanes in these simulations are due to an optical depth effect and they do not coincide with downflows. Downflows, however, occur near the end points of the lanes. High resolution observational studies of the correlation and relative position of downflows and dark lanes may help clarify this issue in the future.
It should be added that while observations suggest that the properties of central umbral dots are different from those of peripheric umbral dots (Sobotka & Jurčák 2009), it is hard to draw a sharp line between the two groups, and the finite tilt of the magnetic field may play an important role in determining the properties of central UDs as well. Nonvertical fields, e.g. those due to a fanning out or twisting of the flux rope, may also be instrumental in inhibiting the appearance of FXS inside sunspot umbrae – though this may also be explained simply by the strong fields in sunspot umbrae or by the large depth (and consequently low aspect ratio) of the umbra.
These considerations suggest that one possible way to bring idealized models closer to realistic numerical simulations could be to extend these results to the case where the imposed magnetic field is tilted.
Available at http://astro.elte.hu/~kris/magconvmovies/animations
Acknowledgments
We thank the anonymous referee for constructive comments. The computations were tested and performed on Eötvös University’s 756core, 3.7 Tflop HPC cluster Atlasz and on the supercomputing network of the NIIF Hungarian National Supercomputing Centre (project ID: 1117 fragment). Support by the Hungarian Science Research Fund (OTKA grants no. K81421 and K83133), by the European Commission’s 6th Framework Programme (SOLAIRE Network, MTRNCT2006035484), as well as by the European Union with the cofinancing of the European Social Fund (grant no. TÁMOP4.2.1/B 09/1/KMR20100003) is gratefully acknowledged. CT acknowledges the financial support from University of Helsinki research project “Active Suns”.
References
 Bharti, L., Jain, R., & Jaaffrey, S. N. A. 2007a, ApJ, 665, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Bharti, L., Joshi, C., & Jaaffrey, S. N. A. 2007b, ApJ, 669, L57 [NASA ADS] [CrossRef] [Google Scholar]
 Bharti, L., Joshi, C., Jaaffrey, S. N. A., & Jain, R. 2009, MNRAS, 393, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Bharti, L., Beeck, B., & Schüssler, M. 2010, A&A, 510, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blanchflower, S. 1999, Phys. Lett. A, 261, 74 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Blanchflower, S., & Weiss, N. 2002, Phys. Lett. A, 294, 297 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chan, K. L., & Sofia, S. 1986, ApJ, 307, 222 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford: Clarendon Press) [Google Scholar]
 de Wijn, A. G., Stenflo, J. O., Solanki, S. K., & Tsuneta, S. 2009, Space Sci. Rev., 144, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Guenther, D. B., Demarque, P., Kim, Y., & Pinsonneault, M. H. 1992, ApJ, 387, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Hamedivafa, H. 2011, Sol. Phys., 45 [Google Scholar]
 Hirzberger, J., Bonet, J. A., Vázquez, M., & Hanslmeier, A. 1999, ApJ, 527, 405 [NASA ADS] [CrossRef] [Google Scholar]
 Houghton, S. M., & Bushby, P. J. 2011, MNRAS, 412, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Kilcik, A., Yurchyshyn, V. B., Rempel, M., et al. 2012, ApJ, 745, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Kitai, R., Watanabe, H., Nakamura, T., et al. 2007, PASJ, 59, 585 [Google Scholar]
 Moradi, H., Baldner, C., Birch, A. C., et al. 2010, Sol. Phys., 267, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Ortiz, A., Bellot Rubio, L. R., & Rouppe van der Voort, L. 2010, ApJ, 713, 1282 [NASA ADS] [CrossRef] [Google Scholar]
 Petrovay, K. 1990, ApJ, 362, 722 [NASA ADS] [CrossRef] [Google Scholar]
 Rast, M. P. 2003, ApJ, 597, 1200 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2012, ApJ, 750, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M., Schüssler, M., Cameron, R. H., & Knölker, M. 2009a, Science, 325, 171 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Rempel, M., Schüssler, M., & Knölker, M. 2009b, ApJ, 691, 640 [NASA ADS] [CrossRef] [Google Scholar]
 Riethmüller, T. L., Solanki, S. K., & Lagg, A. 2008a, ApJ, 678, L157 [NASA ADS] [CrossRef] [Google Scholar]
 Riethmüller, T. L., Solanki, S. K., Zakharov, V., & Gandorfer, A. 2008b, A&A, 492, 233 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rimmele, T. 2008, ApJ, 672, 684 [NASA ADS] [CrossRef] [Google Scholar]
 Schüssler, M., & Vögler, A. 2006, ApJ, 641, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Sobotka, M., & Hanslmeier, A. 2005, A&A, 442, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sobotka, M., & Jurčák, J. 2009, ApJ, 694, 1080 [NASA ADS] [CrossRef] [Google Scholar]
 Sobotka, M., & Puschmann, K. G. 2009, A&A, 504, 575 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 2000, Sol. Phys., 192, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Tao, L., Proctor, M. R. E., & Weiss, N. O. 1998, MNRAS, 300, 907 [NASA ADS] [CrossRef] [Google Scholar]
 Tian, C. L., Xu, K., Chan, K. L., & Deng, L. C. 2007, J. Comput. Phys., 226, 2003 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., & Kondo, M. 1989, PASJ, 41, 197 [NASA ADS] [Google Scholar]
 Watanabe, H., Bellot Rubio, L. R., de la Cruz Rodríguez, J., & Rouppe van der Voort, L. 2012, ApJ, 757, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Weiss, N. O., Brownjohn, D. P., Matthews, P. C., & Proctor, M. R. E. 1996, MNRAS, 283, 1153 [NASA ADS] [CrossRef] [Google Scholar]
 Weiss, N. O., Proctor, M. R. E., & Brownjohn, D. P. 2002, MNRAS, 337, 293 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Statistical properties of steady unmagnetized convection, used as initial condition for magnetoconvection. Shown are the averaged energy fluxes, in units of 6.34 × 10^{11}erg/s/cm^{2}. Averages are taken in time and in the horizontal plane. z is the height from bottom, in units of 1.72 Mm. 

In the text 
Fig. 2 Histories of magnetic energy for various runs, as indicated by the labels. 

In the text 
Fig. 3 Dopplergrams at various depths for selected cases. Corresponding magnetic fields are shown in Fig. 4. Velocities are given in units of 6.03 km s^{1}. 

In the text 
Fig. 4 Contour plots of vertical magnetic fields at various depth for selected cases. Corresponding vertical flow fields are shown in Fig. 3. Magnetic fields are given in units of 2911.55 G. 

In the text 
Fig. 5 Vertical cuts of the computational box along x = 3 at an arbitrary instant. The colour scale represents temperature fluctuation from its horizontal mean, i.e., , in units of 5063.52 K. Rows correspond to different runs as labeled on the left. 

In the text 
Fig. 6 Magnetic field strength  B  corresponding to the cuts in Fig. 5. The unit is 2911.55 G. 

In the text 
Fig. 7 Root mean square (rms) velocity and magnetic fields from their horizontal means for case F. z is the height from bottom, in units of 1.72 Mm. Units for velocity and magnetic field strength are 6.02 km s^{1} and 2911.55 G, respectively. 

In the text 
Fig. 8 Time evolution of vertical velocity at z = 0.8, near the top of the convective layer. The time intervals between snapshots are ~71 min for cases D, E, F, G; ~86 min for case H. Velocities are given in units of 6.03 km s^{1}. 

In the text 
Fig. 9 PDFs for the distribution of vertical Alfvénic speed, , at z = 0.8 and z = 0.4 for various imposed vertical magnetic field strengths as indicated inside the labels, corresponding to, totally suppressed (black solid), SSC (blue dotted), FXS (magenta dashed) and GRC (red dash dot). 

In the text 
Fig. 10 PDFs for the distribution of vertical velocity, i.e.,v_{z}, at z = 0.8 and z = 0.4 for various imposed vertical magnetic field strengths as indicated inside the labels: totally suppressed field (thick black solid), SSC (blue dotted), FXS (magenta dashed), and GRC (red dash dot). The blue solid thin is the PDF of the unmagnetized case. 

In the text 
Fig. 11 Vertical anisotropy of the velocity field (upper panel) and magnetic field (in Alfvénic units, lower panel), for cases A (thin black solid), B (green dotted), D (cyan dashed), E (blue dashdot), F (purple dash–3dot), G (magenta long dashes), and H (red solid). z is the height from bottom, and ⟨ ... ⟩ represents horizontal average. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.