Simulation of the Formation of a Solar Active Region

Simulation of the Formation of a Solar Active Region

M. C. M. Cheung Lockheed Martin Solar and Astrophysics Laboratory, Palo Alto, CA 94304, USA. M. Rempel High Altitude Observatory, NCAR, P.O. Box 3000, Boulder, Colorado 80307, USA A. M. Title Lockheed Martin Solar and Astrophysics Laboratory, Palo Alto, CA 94304, USA. M. Schüssler Max Planck Institute for Solar System Research, Katlenburg-Lindau, 37191, Germany.

We present a radiative magnetohydrodynamics simulation of the formation of an Active Region on the solar surface. The simulation models the rise of a buoyant magnetic flux bundle from a depth of 7.5 Mm in the convection zone up into the solar photosphere. The rise of the magnetic plasma in the convection zone is accompanied by predominantly horizontal expansion. Such an expansion leads to a scaling relation between the plasma density and the magnetic field strength such that . The emergence of magnetic flux into the photosphere appears as a complex magnetic pattern, which results from the interaction of the rising magnetic field with the turbulent convective flows. Small-scale magnetic elements at the surface first appear, followed by their gradual coalescence into larger magnetic concentrations, which eventually results in the formation of a pair of opposite polarity spots. Although the mean flow pattern in the vicinity of the developing spots is directed radially outward, correlations between the magnetic field and velocity field fluctuations allow the spots to accumulate flux. Such correlations result from the Lorentz-force driven, counter-streaming motion of opposite-polarity fragments. The formation of the simulated Active Region is accompanied by transient light bridges between umbrae and umbral dots. Together with recent sunspot modeling, this work highlights the common magnetoconvective origin of umbral dots, light bridges and penumbral filaments.

1 Introduction

Active regions (ARs) are the most prominent manifestation of the large scale solar magnetic field in the photosphere of the Sun. Understanding the underlying flux emergence process is a crucial step toward understanding the connection between the dynamo processes in the solar convection zone and magnetic activity in the photosphere and above.

Observations of active regions on the solar surface indicate that, in terms of magnetic flux content, there is a continuous spectrum of active region sizes (Zwaan, 1978, 1987; Harvey & Martin, 1973; Schrijver & Zwaan, 2000; Hagenaar, 2001; Hagenaar et al., 2003). At the small end of the spectrum, ephemeral ARs (which have a flux of Mx) manifest themselves in the solar photosphere in terms of mixed polarity field and magnetic bright points but not necessarily pores (Cheung et al., 2008)Small ARs, which have fluxes of to Mx in each polarity, contain solar pores. Large ARs, which have even more magnetic flux (up to Mx), contain sunspots with penumbrae.

Since the original suggestion by Parker (1955), there has been a large body of theoretical and observational work supporting the scenario that sunspots (and ARs) form as a result of the buoyant rise of magnetic flux bundles from the solar convection zone to the solar atmosphere. On the observational side, the work of Zwaan (1978, 1987) and Strous & Zwaan (1999) showed that ARs form as a consequence of a succession of small-scale flux emergence events on the solar surface. The observational studies in these papers noted that although magnetic flux emerges as small bundles, these bundles subsequently migrate and coalesce in a systematic fashion, eventually leading to formation of pores and sunspots (i.e. patches of predominantly one polarity). Recently, an analysis of spectropolarimetric observations by Schlichenmaier et al. (2010) shows how the accumulative of flux by a solar pore resulted in the formation of penumbral filaments.

Numerical magnetohydrodynamics (MHD) simulations of the rise of buoyant magnetic flux tubes have yielded a wealth of information regarding the flux emergence process. For example, numerical models of the rise of toroidal flux tubes from the bottom of the convection zone based on the thin flux tube approximation (Roberts & Webb, 1978; Spruit, 1981) are able to reproduce global properties of ARs such as their tilt angles (D’Silva & Choudhuri, 1993), their emergence latitudes (Caligari et al., 1995, 1998) and asymmetries between the leading and following polarities (Fan et al., 1993; Moreno-Insertis et al., 1994; Caligari et al., 1995, 1998). Due to the invalidity of the thin flux tube approximation in the uppermost ten or so Mm of the solar convection zone, such studies have limited applicability for examining how ARs form on the solar surface. Complementing this line of work, multidimensional MHD simulations (e.g. Forbes & Priest, 1984; Isobe et al., 2005; Shibata et al., 1989; Archontis et al., 2004; Manchester IV et al., 2004; Magara, 2006; Fang et al., 2010) into the solar atmosphere have greatly improved our understanding of how the buoyant rise of magnetic fields into the solar atmosphere affect local dynamics. For a historical perspective of such work, the reader is referred to the review articles by Moreno-Insertis (1997)Fisher et al. (2000) and Fan (2004, 2009).

Recent models that include the treatment of convective motions and radiative heating/cooling (either by solving the radiative transfer equation or by parameterized terms in the energy equation) have shown these two effects to be very important to the properties of flux emergence (Stein & Nordlund, 2006; Cheung et al., 2007, 2008; Isobe et al., 2008; Abbett, 2007; Martínez-Sykora et al., 2008, 2009; Tortosa-Andreu & Moreno-Insertis, 2009; Yelles Chaouche et al., 2009; Fang et al., 2010). Despite these advances, a comprehensive numerical model of the formation of an AR resulting from flux emergence has been lacking. This is a consequence of the vast range of length and time scales encountered in the convection zone as consequence of a density variation by about six orders of magnitude from the base to the top. To cope with this difficulty numerical models focus either on the deep convection zone leaving out the upper most 10 - 20 Mm or on the upper most 10 Mm including the photosphere.

While the former typically utilize the anelastic approximation (Abbett et al., 2000; Fan et al., 2003; Jouve & Brun, 2009) the latter have to be based on fully compressible MHD. Regardless of the differences a common challenge in both cases is to explain the formation of rather coherent sunspots with several kG field strength from magnetic field which has risen through a highly stratified atmosphere and consequently should have weakened considerably through expansion. Addressing this aspect through a self-consistent numerical simulation capturing the rise of magnetic flux from the base of the convection zone into the photosphere is currently out of reach, however much can be learned from a numerical simulation containing the last Mm beneath the photosphere. While the geometric extent is only of the convection zone depth, the density drop of orders of magnitude (corresponding to ten pressure scale-heights) is comparable to the drop in the remaining of the convection zone. In this paper we investigate the flux emergence process in the upper most Mm of the convection zone resulting from the advection of a semitorus-shaped flux tube across the bottom boundary of our computational domain.

2 Simulation Setup

The radiative MHD simulation was carried out with the MURaM code (Vögler et al., 2005; Rempel et al., 2009b), which takes into account radiative energy transport in the energy equation (assuming local thermodynamic equilibrium) and a realistic equation of state. The MURaM code has been used to model a variety of magnetic structures in the solar photosphere and convection zone. The Cartesian domain spans Mm Mm in the two horizontal directions and Mm in the vertical direction (with horizontal and vertical grid-spacing of and km respectively). The base of the photosphere (, the mean geometrical height where the Rosseland optical depth has a value of unity) is located Mm above the bottom boundary of the simulation domain. Before the introduction of a magnetic flux tube, a purely hydrodynamic simulation was performed to allow the radiatively driven convection to relax to statistical equilibrium.

To mimic the rise of magnetic flux into the topmost layers of the convection zone, an axisymmetric, twisted flux tube with the shape of a semitorus is kinematically advected into the domain through the bottom boundary. Similar time-dependent boundary conditions to kinematically advect magnetic field into the top layers of the convection zone have previously been used for flux emergence simulations. Stein & Nordlund (2006) and Stein et al. (2010) used a boundary condition which advect purely horizontal field (of uniform field strength and orientation) through upflow regions. Martínez-Sykora et al. (2008) used a boundary condition which advected an horizontal twisted flux tube. The boundary condition for the current simulation is implemented by specifying the time-dependent fluxes of vertical mass, momentum, energy and magnetic field (i.e. terms acted upon by the divergence operator in the conservation form of the MHD equations) consistent with the rise of the torus. Following Fan & Gibson (2003), the magnetic field distribution of the flux tube in a spherical coordinate system with the torus centered at the origin is given by


is the distance of a point from the toroidal axis of the tube, and and are the semi-minor and -major radii of the torus, respectively. The dimensionless twist parameter is equivalent to the parameter used by Fan & Gibson (2003). For this simulation, we took , Mm and Mm. kG such that the toroidal field strength at the tube axis is kG (corresponding to a plasma of ). For , is set to zero. The total toroidal flux content of the tube is Mx. Averaged over the cross-section of the tube (i.e. ), the mean toroidal field strength is kG. This is somewhat stronger, but still comparable to field strengths of a few kG expected from thin flux tube simulations that begin with kG at the base of the convection zone (see, e.g., Caligari et al., 1995). Thin flux tube simulations also predict a rise speed of between and km s at this depth. For this simulation, we chose to impose a rise speed of km s.

The axis of symmetry of the torus is parallel to the -axis of the Cartesian domain so that the resulting magnetic loop in the domain has its axis lying in the - plane. At the bottom boundary, the gas pressure within the torus is given by , where dyne cm is the mean gas pressure at the bottom boundary. The inflow specific entropy of the plasma in the torus is set to the same value as that of ambient convective upflows. This choice of the thermodynamic properties results in a relative density deficit of ( is the first adiabatic exponent) at the tube axis.

Beginning at , the magnetic torus is kinematically advected into the domain with a constant rise speed of km s until the top half of the torus has traversed through the bottom boundary. Beyond this point in time (i.e., hrs), the velocity at the bottom boundary within the tube is set to zero. Outside of the tube, the bottom boundary condition allows for smooth inflows and outflows (see Vögler et al. 2005 for details). Mass flows across the top boundary are allowed (upflows are unimpeded and downflows are damped). By virtue of the low densities near the top boundary, the mass flux across the top has a negligible effect on the mass content in the simulation domain. The magnetic field at the top boundary is matched to a potential field. Periodic boundary conditions apply at the side boundaries. Following Rempel, Schüssler, & Knölker (2009b), we limit the strength of the Lorentz force in low- regions () to limit the Alfvén speed to a maximum of km s. This artificial limiting of the Lorentz force mainly acts within the central umbral regions of sunspots in mid-to-high photosphere, where local Alfvén speeds reach up to km s and would impose prohibitively small time-steps on the explicit scheme of the code.

Recently, MacTaggart & Hood (2009) modeled the rise of buoyant toroidal magnetic flux tubes. Their simulation setup is different from ours in a number of ways. First of all, they used a plane-parallel background convection zone and atmospheric model to mimic the stratification in the Sun but do not include convective flows nor radiative transfer. Secondly, the toroidal flux tubes in their simulations were embedded as buoyant, stationary structures as an initial condition. An interesting result from their calculations is that the choice of a toroidal flux tube (as opposed to an originally horizontal tube) facilitated the emergence of the axis of the tube into the photosphere. MacTaggart & Hood (2009) attributed this to the ability (in the toroidal tube case) of plasma to drain down the legs of of the tube to enhance buoyant at the tube apex. In section 3.3, we showed that convective flows acting on the rising magnetic plasma have a similar effect.

3 Simulation Results

3.1 Photospheric Emergence and Active Region Formation

Before proceeding to discuss the physics underlying the flux emergence and AR formation process, we describe how the emerging flux region evolves in time. Magnetic flux begins to emerge into the photosphere beginning at hrs and emergence progresses for several hours. A time sequence of the grey intensity and maps of the vertical field (sampled at ) are shown in Fig. 1. In the intensity map at hr, elongated granules begin to appear. This stretching of the convective cells is due to the horizontal expansion of the rising magnetic plasma. A more detailed discussion of the physics of horizontal expansion is given in section 3.2. Due to the subsurface horizontal expansion and undulation of the field lines by convective flows, the emerging flux region covers an extend surface area within which small-scale bipoles emerge with a systematic orientation (see also Cheung et al., 2008). By hr, small pores begin to appear in multiple locations in the intensity map (c.f. Stein et al., 2010). A corresponding synthetic ‘magnetogram’ ( sampled at ) is shown in Fig. 2. Two hours later, two spots have appeared near the horizontal positions coincident with where the subsurface footpoints of the initial semitorus are rooted. As time progresses, the intensity maps show the spots growing in size. A synthetic magnetogram at a later time shows two large, vertical concentrations of opposite polarity field at the locations of the spot (see Fig. 3). Throughout the simulation, the two spots never attain fully developed penumbrae which is consistent with observations (Zwaan, 1978, 1987) since the two large magnetic concentrations have (over the course of the simulation) at most Mx. Nevertheless, transient penumbral filaments, umbral dots and even light bridges are found in the synthetic intensity maps.

3.2 Subsurface Rise and Expansion of the Magnetic Tube

Figure 4 shows the structure of the torus as it begins to erupt from the convection zone into the overlying photosphere. The subsurface roots of the torus at the bottom of the computational domain are shown in the magnetic map at a depth of 7.5 Mm (left panel). Field lines traced from the opposite polarity roots diverge upwards to covered an extended horizontal region. The right panel shows the associated emerging flux region at the surface, where the magnetic field lines have been undulated by their interaction with the convective downflows as is consistent with observations (see e.g. Pariat et al., 2004; Cheung et al., 2007; Centeno et al., 2007; Martínez González et al., 2007; Martínez González & Bellot Rubio, 2009; Cheung et al., 2008; Ishikawa et al., 2010).

The number of pressure scale-heights spanned between Mm and (i.e. the photospheric base) is . To reach pressure equilibrium with the surroundings, the rising magnetic structure expands strongly in the horizontal directions. This is an important stage in the evolution for reasons of hydrostatic balance and mass conservation. Since the rise time of the flux tube is much longer than both the free-fall and sound-crossing times (both s) across the layers of the convection zone captured in this simulation, its rise does not significantly alter the mean stratification. On the other hand, the injection of mass from the rise of the half-torus is equal to of the original mass in the entire domain, which is roughly equal to the total mass content of the top Mm of the convection zone (as captured in the computational domain). Thus it is not surprising that the rising magnetic plasma displaces a substantial fraction of the original mass in the domain and fills the near-surface layers during the first stage of the active region formation process. As a consequence, the top of the magnetic torus takes on a flattened, pancake-like structure in near-surface layers (see Fig. 5). A similar behavior is also reported in the flux emergence simulations by Toriumi & Yokoyama (2010).

Given the field strength of the original magnetic torus at a depth of Mm, what strengths can we expect for the field emerging into the photosphere? This question can be addressed by considering a number of different simplified scenarios. For a purely horizontal magnetic flux tube expanding in directions transverse to its axis, conservation of mass and conservation of magnetic flux give the linear scaling relation . A different scaling relation can be derived by considering how an upwelling in the stratified convection zone modifies horizontal field threading the rising plasma. Let the initial horizontal field be and for simplicity assume that the upflow is axisymmetric about the (vertical) -axis and centered at the origin (the absolute coordinate is unimportant). Let the rates of expansion in the horizontal and vertical directions be and respectively. Here is a measure of the flow anisotropy ( for isotropic expansion). From the ideal Induction Equation, the rate of change of the magnetic field attached to a Lagrangian fluid element is


which in this case reduces to


The corresponding Lagrangian form of the continuity equation is


Combining equations (5) and (7), we obtain the scaling relation


For isotropic expansion (), . In the case where the horizontal expansion rate is much larger than that of the vertical expansion rate we have , so that .

For a field strength of kG at g cm (density at Mm), the linear scaling relation gives, for g cm, a photospheric field strength of G. This is evidently too weak compared to what is found in the simulation and to what is observed, both of which give emerging horizontal field strengths on the order of a few hundred G (Lites et al., 1998; Kubo et al., 2003).

The second scaling relation yields more consistent values for emerging horizontal field strengths. For (zero expansion in the vertical direction) and (isotropic expansion), a density drop of yields horizontal photospheric field strengths of G and G, respectively. Figure 6 shows a joint probability density function (JPDF) of the horizontal field strength versus mass density sampled at the mid-plane between and hours. Points with low values of specific entropy () have been excluded since they correspond to plasma that has undergone radiative cooling at the surface. The maximum magnetic field strength of kG corresponds to the field strength at the axis of the torus that is introduced through the bottom boundary. The solid and dashed lines correspond to the scaling relations and respectively. The latter scaling relation is clearly too steep and underestimates the field strengths at photospheric densities ( g cm). The square-root scaling relation provides a much better match to the JPDF. This result is consistent with the fact that the plasma experiences predominantly horizontal expansion ( close to zero) during its rise to the surface (see Fig. 5).

Figure 7 shows how the arrival of buoyant, magnetic plasma at the photosphere leads to a transient pressure excess which drives diverging horizontal flows about the flux emergence site. The coincidence (spatially and temporally) of this pressure enhancement with the onset of large-scale diverging flows (beginning at about hrs) is evidence that the flows are driven by the associated horizontal pressure gradient. is of the order , so that the horizontal flows have Mach numbers of similar amplitude. After the initial emergence phase ( hr), the divergent horizontal flows at the periphery of the magnetic complex persist, albeit with lower amplitude ( km s). This raises the question of how magnetic flux is still capable of accumulating in the two spots, which will be addressed in section 3.4.

3.3 Discharge of Mass from the Rising Magnetic Structure

Given that the near-surface field is rather dispersed, the next stage of active region evolution is the coalescence of the dispersed field into coherent concentrations. In order for the fragments to collect together, excess mass must be unloaded from the emerged field lines. To illustrate the amount of mass that is eventually removed, we compare the mass carried within the original half-torus and the mass remaining within the two opposite polarity concentrations at hrs. Let and be azimuthal averages of one such flux concentration about its vertical axis. Using these, we calculated surfaces of constant magnetic flux and the corresponding mass content contained within the flux surfaces. The mass contained within a volume defined by the flux surface Mx (summed over both spots) at hrs is only of the mass contained in a sub-volume of the original half-torus with the same flux content. Thus most of the original injected mass has somehow escaped to the surroundings.

A number of physical mechanisms may be invoked to explain the discharge of mass from the original magnetic structure. First of all, magnetic diffusion of the torus allows mass transfer across magnetic field lines. However, this would result in an increase of mass in the semitorus rather than a decrease and thus can be ruled out as the underlying mechanism for mass discharge. As a second possibility, we consider the mass flux through the top and bottom boundaries. By virtue of the low densities at the top boundary, the vertical mass flux there is negligible for the mass budget of the semitorus. At the bottom boundary, vertical (and horizontal) velocities within the torus () are set to zero after the initial prescribed rise. Downflows are permitted at other magnetic footpoints at the bottom boundary but since the majority of the magnetic flux at that layer remains within the torus, mass flux through the bottom boundary can also be ruled out as the main discharge mechanism. A third possibility is the removal of mass via outflows associated with the horizontal expansion (see section 3.2). While it is true that outflows carry mass away from the center of the emerging region, these flows also advect magnetic field with them and cause a weakening of the field. So outflows alone are insufficient to explain how mass is removed from the field lines.

The responsible mass removal mechanism is illustrated in Fig. 8. The left panel of this figure shows a schematic representation of convective flows undulating a field line which has one end anchored deep below the photosphere. As already reported in Cheung et al. (2008), expulsion of flux from the granular centers lead to the encounter and cancellation of opposite polarity field. Such opposite polarity pairs are connected below the surface in the form of U-loops (see also Stein et al., 2010). The first consequence of this magnetic reconnection between the opposite polarities is the creation of an O-loop, which discharges mass from the original field line. This is akin to the process suggested by Parker (1994) for the discharge of mass from and intensification of magnetic fields at the bottom of the convection zone. Based on observations of emerging flux regions with Spectro-Polarimeter instrument (Lites et al., 2001) of the Solar Optical Telescope (Tsuneta et al., 2008) onboard the Hinode satellite (Kosugi et al., 2007)Lites (2009) put forward this mechanism to explain how mass is discharged from magnetic field lines in emerging flux regions. In the simulation, mass is subducted in a similar fashion though in the 3D case plasmoids take the place of O-loops.

The amount of vertical unsigned flux at the surface may erode depending on the height at which reconnection occurs. Specifically, if reconnection happens exactly at the surface, the unsigned flux will immediately decrease. If reconnection occurs above , the unsigned flux will decrease only when resulting O-loops (plasmoids) submerge. This second scenario has indeed been reported by Iida, Yokoyama, & Ichimoto (2010). If reconnection occurs below the surface, the surface unsigned flux will decrease if the remnant U-loop (counterpart to the O-loop after reconnection) rises above the surface. This is less likely since the U-loop will still be anchored to downflows. We find in our simulation that reconnection occurs both in the photosphere and in the convection zone. However, the unsigned flux (at any one horizontal level) erodes due to retraction of inverse O-loops (plasmoids). The process described has been numerically modeled by Isobe, Tripathi, & Archontis (2007, in 2D) and by Archontis & Hood (2009, in 3D) for explaining the origin of Ellerman bombs resulting from magnetic reconnection occurring in emerging flux regions (Pariat et al., 2004, 2009).

The right panel of Fig. 8 shows an example of a collection of sinking U-loops in the simulation. The downward transport of such U-loops, which allows the discharge of mass from the rising magnetic field, is related to the phenomenon known as turbulent pumping (Tobias et al., 1998; Brummell et al., 2002). Turbulent pumping manifests itself in terms of vertical component the Poynting flux of magnetic energy, which for ideal MHD is


where the first term in the integrand represents the bodily ascent (or descent) of horizontal fields and the second term represents the Poynting flux due to horizontal flows. Figure 9 shows plots of the Poynting flux through two horizontal planes ( and Mm) as functions of time. The individual contributions from vertical and horizontal flows are also plotted. The traversal of the rising magnetic structure through the Mm plane appears as a positive hump in both the emergence term and the combined Poynting flux between and hr. Thereafter the effect of turbulent pumping takes over and the net Poynting flux due to vertical flows is negative. At the surface ( Mm), the convective downflows are much stronger (in terms of Mach number, they approach ) and the effect of turbulent pumping on the Poynting flux is pronounced throughout the emergence and post-emergence phase. It should be noted that although the Poynting flux (a flux of energy) has a negative sign, it does not necessarily imply that net magnetic flux itself is migrating downwards.

3.4 Growth of Magnetic Flux Content in the Developing Spots

As already discussed in the previous section, mean horizontal outflows persist through the simulated duration of AR formation. Given such mean flows, how does the magnetic flux at the photosphere migrate inwards to yield coherent spots? To examine this, the magnetic and velocity fields (for a fixed time) were decomposed into the sum of azimuthal averages and corresponding fluctuating components, namely etc, where the bar and prime denote an azimuthal average and the fluctuation about the average, respectively. For ideal MHD, the mean field induction equation is


where is the (azimuthally-averaged) mean-field electromotive force resulting from correlations between the magnetic field and velocity field fluctuations (Krause & Rädler, 1980; Rädler, 1980). To examine how (vertical) magnetic flux accumulates in the growing spots, consider a circular area of fixed radius , at height and centered at . The time rate of change of the total magnetic flux crossing this circular area is




Equation (14) indicates that the magnetic flux enclosed within changes as a result of (1) advection of the mean magnetic field by the mean flow (first term on the r.h.s.), and (2) correlations between their fluctuating components.

Figure 10 shows plots of the mean magnetic field and of the flux transport rates defined in Eq. (14). To smooth out fluctuations in time, the azimuthal means and azimuthally-fluctuating quantities were also averaged between and hr in time. This remains consistent with the derivation above since the operations of averaging temporally and azimuthally commute. At this stage of the simulation, the mean magnetic field within a radius of Mm has already reached kG values and we focus our attention on how the weaker field at larger radii is transported. Inspection of the mean velocity field reveals that the contribution from the radially directed outflow () dominates such that magnetic flux is transported away from the spot. However, it is found that correlations between the velocity and magnetic fields (i.e. ) provide a means for the inward migration of vertical field. The amplitude of is in fact larger than that of and the net result is an accumulation of flux in the spot.

The streaming of opposite polarity (vertical) field in opposite horizontal directions (Strous et al., 1996; Strous & Zwaan, 1999; Bernasconi et al., 2002; Cheung et al., 2008) provides the systematic correlation needed for accumulation of magnetic flux in the spot. This begs the question of what is the underlying driver of the counter-streaming motion. The left panel of Figure 8 shows how the granular expulsion (Weiss, 1966; Hurlburt & Toomre, 1988) of magnetic flux from serpentine field lines leads to a migration of positive polarity in one horizontal direction and negative polarity field in the opposite direction (i.e. a non-zero ). The presence of sea serpents, however, is an insufficient condition for flux accumulation in spots. For instance, consider a scenario wherein uniform magnetic field (with a net horizontal component) threads a series of granules. Although is non-zero, the translational symmetry of the setup causes the loop integral to vanish so that there is no net change of flux.

The missing ingredient is some underlying large-scale structure which remains coherent over time-scales much greater than the granulation turnover times at the surface. Here it is provided by the coherent subsurface roots of the nascent active region (which is illustrated by the anchored field line in the left panel of Fig. 8). The upper four panels of Fig. 11 show radial profiles of the azimuthally-, temporally- and spatially- (over a height of km about ) averaged forces about the developing spot. For both positive and negative polarities, the gas and magnetic pressure gradient forces roughly (but not exactly) balance each other. The profile for the magnetic tension force () shows a tendency to accelerate positive (negative) polarities inward (outward). This trend is even clearer when one looks at the net force (-) and is consistent with the counter-streaming behavior of the opposite polarities (see bottom panel of Fig. 11). From this analysis, we conclude that the organization of the surface flux fragments into the two spots results from the presence of the coherent surface roots which influence the surface dynamics via the Lorentz force (especially the tension force).

This result closely resembles the tethered-balloon model of Spruit (1981). In this model, emerging magnetic loops consist of buoyant gas in the crests and neutrally-buoyant or anti-buoyant gas at footpoints anchored below the photosphere. Like tethered balloons, the equilibrium state is reached when buoyant elements hover vertically above their footpoints. The evolution of a field-line segment from the oblique to the vertical is mediated by the tension force. In our case, the interaction of granular convective flows with emerging magnetic fields, resulting in the undulation of field lines, flux expulsion to intergranular lanes and subsequent convective intensification (Cheung et al., 2008; Danilovic et al., 2010) add a layer of complexity to the problem. Nevertheless, the present analysis demonstrates that Spruit’s model is qualitatively compatible with our simulation results.

3.5 Light Bridge Formation and Disappearance

As coalescence of flux concentrations progresses, ambient material with relatively weak field can become entrained in between adjacent high field-strength concentrations. An extended lane of entrained material trapped between neighboring regions of strong field appears as a light bridge in intensity maps. Figure 12 shows an example from the simulation. In this case, the horizontal convergence of strong field concentrations (manifested as dark umbrae in continuum intensity) traps ambient material with weak field to form a light bridge. Three vertical cross-sectional cuts oriented perpendicular to the light bridge show the vertical velocity distribution and magnetic field strength along parts of the bridge. At its end, located deep inside the spot ( Mm), the light bridge has a width of a few hundred km. The corresponding cross-sectional plots of and show that it is associated with an upflow of plasma with relatively weak magnetic field (a few hundred G at maximum). The upwelling locally lifts the surface by about km relative to the adjacent umbra. The expansion with height of the neighboring magnetic field pinches the upflow to force a cusp-like structure near  (Nordlund, 2004). This is consistent with the observational study of light bridges by Lites et al. (2004). Just below , the pinching accelerates the upflow from km s to a few km s. Above , the upflow speed decreases. This stagnation-like point partially traps low entropy material (radiatively cooled) and leads to a central dark lane within the light bridge. The slower speeds at the surface are consistent with the observational results by Rimmele (1997). Cool material escape by feeding the downflows flanking the upflow. All these processes are similar to those occurring in simulations of umbral dots (Schüssler & Vögler, 2006).

Vertical cross-sections across wider segments of the light bridge (Fig. 12) tell a similar story. The vertical cross-section given in the middle panel with weak field, also flanked by a pair of periphery downflows. At an even wider section of the bridge with a width of a few Mm, the convective flow consists of a few overturning cells instead of a single one.

As the strong field concentrations continue to converge, the light bridge eventually disappears (at around hrs) by means of outflows channeled toward the edge of the spot. However, until the time when we halted the simulation ( hrs), transient light bridges and umbral dots continued to appear in the pair of spots due to the tendency of entrained material to attempt to escape. Towards the later stages of the simulation, however, the frequency (and size) of light bridges diminished together with the amount of entrained material available for light bridge formation.

Together with previous models of umbral convection (Schüssler & Vögler, 2006) and sunspot simulations (Heinemann et al., 2007; Rempel et al., 2009a, b) suggest that umbral dots, penumbral filaments and light bridges are all manifestations of overturning magnetoconvection.

4 Discussion

In this paper, we presented the first radiative MHD simulation of the birth of an AR on the solar surface. To mimic the emergence of magnetic flux from the convection zone to the photosphere, a magnetic semitorus was kinematically advected upward through the bottom boundary, which is located at Mm below the photospheric base. Here is a list of interesting findings regarding the physics of how AR formation occurs in the simulation.

The magnetic field strength of the rising plasma in the emerging flux region roughly follow the scaling relation , where is the mass density. This scaling relation is consistent with the scenario that rising plasma preferentially expand in the horizontal directions. This horizontal expansion disperses field over a large horizontal area. As time progresses, the dispersed magnetic fragments coalesce into gradually larger concentrations (in terms of flux content). When compact magnetic concentrations attain a flux of Mx or more, they appear as dark pores in intensity images.

In order for dispersed field to become compact again, excess mass must be removed from the original emerging structure. This is facilitated by convective downflows, which act on horizontal field to form U-loops anchored below the surface. Magnetic reconnection within these loops transports mass away from rising field lines. This physical mechanism is the same as the one suggested by Lites (2009) for explaining how the mass is discharged from the rising magnetic field in observed emerging flux regions.

Following the removal of mass from the original magnetic structure, magnetic flux must migrate into two concentrated trunks to form spots. It is found that the azimuthally-averaged flows around the developing magnetic concentrations (corresponding to spots) contribute to the erosion of magnetic flux. However, correlations between the fluctuations of velocity and magnetic fields result in net inward migration of magnetic flux with the same sign as that of the spot. Such correlations result from the counter-streaming motion of opposite polarity fragments. This re-organization of the dispersed flux fragments is due to the presence of coherent subsurface magnetic roots, which influence the surface dynamics via the Lorentz force. In this sense, our result is akin to the tethered-balloon model of sunspot formation suggested by Spruit (1981). We emphasize that these streaming magnetic polarities are more akin to the moving dipolar features (MDFs, Bernasconi et al., 2002) found in emerging flux regions (Strous et al., 1996; Schlichenmaier et al., 2010) than the so-called Type I moving magnetic features (MMFs), which likely also have serpentine field line geometries (Sainz Dalda & Bellot Rubio, 2008; Kitiashvili et al., 2010) but which may be associated with the decay of sunspots (Kubo, Shimizu, & Tsuneta, 2007).

As the dispersed magnetic field collect together to form coherent spots, plasma with relatively weak field (at most a few hundred G) and high entropy may be entrained between regions of predominantly vertical field of a few kG strength. These strong field regions correspond to dark umbrae in intensity maps whereas the entrained weak field plasma in between manifest itself as a light bridge. Cross-sectional cuts of the vertical velocity and magnetic field strength distribution perpendicular to a light bridge in the simulation shares strong resemblance with similar cross-sectional cuts through umbral dots and penumbral filaments. This suggests a common convective origin for all three intensity features.

In this paper, we have focused our attention on analyzing how physical processes occurring in our radiative MHD simulation allow for the formation of an AR with certain observational properties (e.g. elongated granules and mixed polarity patterns in the emerging flux region, pore formation, transient light bridges etc). In a follow-on paper, we intend to carry out more detailed comparisons with observations of emerging flux regions. Such an exercise would be especially timely considering the rise of solar activity and the availability of vector magnetograms from both the Helioseismic Magnetic Imager onboard the Solar Dynamics Observatory and from the Hinode Solar Optical Telescope.

This work was made possible by NASA’s High-End Computing Program. The simulation presented in this paper was carried out on the Pleiades cluster at the Ames Research Center. We thank the Advanced Supercomputing Division staff for their technical support. This work was supported by NASA contract NNM07AA01C at LMSAL. M. Rempel is partially supported through NASA grant NNH09AK02I to the National Center for Atmospheric Research. NCAR is sponsored by the National Science Foundation.


  • Abbett (2007) Abbett, W. P. 2007, ApJ, 665, 1469
  • Abbett et al. (2000) Abbett, W. P., Fisher, G. H., & Fan, Y. 2000, ApJ, 540, 548
  • Archontis & Hood (2009) Archontis, V., & Hood, A. W. 2009, A&A, 508, 1469
  • Archontis et al. (2004) Archontis, V., Moreno-Insertis, F., Galsgaard, K., Hood, A., & O’Shea, E. 2004, ApJ, 441, 886
  • Bernasconi et al. (2002) Bernasconi, P. N., Rust, D. M., Georgoulis, M. K., & Labonte, B. J. 2002, Sol. Phys., 209, 119
  • Brummell et al. (2002) Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825
  • Caligari et al. (1995) Caligari, P., Moreno-Insertis, F., & Schüssler, M. 1995, ApJ, 441, 886
  • Caligari et al. (1998) Caligari, P., Schüssler, M., & Moreno-Insertis, F. 1998, ApJ, 502, 481
  • Centeno et al. (2007) Centeno, R., Socas-Navarro, H., Lites, B., Kubo, M., Frank, Z., Shine, R., Tarbell, T., Title, A., Ichimoto, K., Tsuneta, S., Katsukawa, Y., Suematsu, Y., Shimizu, T., & Nagata, S. 2007, ApJ, 666, L137
  • Cheung et al. (2007) Cheung, M. C. M., Schüssler, M., & Moreno-Insertis, F. 2007, A&A, 467, 703
  • Cheung et al. (2008) Cheung, M. C. M., Schüssler, M., Tarbell, T. D., & Title, A. M. 2008, ApJ, 687, 1373
  • D’Silva & Choudhuri (1993) D’Silva, S., & Choudhuri, A. R. 1993, A&A, 272, 621
  • Danilovic et al. (2010) Danilovic, S., Schüssler, M., & Solanki, S. K., 2010, A&A, 509, 76
  • Fan (2004) Fan, Y. 2004, Living Rev. Solar Phys., 1, 1
  • Fan (2009) —. 2009, Living Reviews in Solar Physics, 6, 4
  • Fan et al. (2003) Fan, Y., Abbett, W. P., & Fisher, G. H. 2003, ApJ, 582, 1206
  • Fan et al. (1993) Fan, Y., Fisher, G. H., & Deluca, E. E. 1993, ApJ, 405, 390
  • Fan & Gibson (2003) Fan, Y., & Gibson, S. E. 2003, ApJ, 589, L105
  • Fang et al. (2010) Fang, F., Manchester, IV, W., Abbett, W. P., & van der Holst, B. 2010, ApJ, 714, 1649
  • Fisher et al. (2000) Fisher, G. H., Fan, Y., Longcope, D. W., Linton, M. G., & Pevtsov, A. A. 2000, Sol. Phys., 192, 119
  • Forbes & Priest (1984) Forbes, T. G., & Priest, E. R. 1984, Sol. Phys., 94, 315
  • Hagenaar (2001) Hagenaar, H. J. 2001, ApJ, 555, 448
  • Hagenaar et al. (2003) Hagenaar, H. J., Schrijver, C. J., & Title, A. M. 2003, ApJ, 584, 1107
  • Harvey & Martin (1973) Harvey, K. L., & Martin, S. F. 1973, Sol. Phys., 32, 389
  • Heinemann et al. (2007) Heinemann, T., Nordlund, Å., Scharmer, G. B., & Spruit, H. C. 2007, ApJ, 669, 1390
  • Hurlburt & Toomre (1988) Hurlburt, N. E., & Toomre, J. 1988, ApJ, 327, 920
  • Iida et al. (2010) Iida, Y., Yokoyama, T., & Ichimoto, K. 2010, ApJ, 713, 325
  • Ishikawa et al. (2010) Ishikawa, R., Tsuneta, S., & Jurcak, J. 2010, ApJ, 713, 1310
  • Isobe et al. (2005) Isobe, H., Miyagoshi, T., Shibata, K., & Yokoyama, T. 2005, Nature, 434, 478
  • Isobe et al. (2008) Isobe, H., Proctor, M. R. E., & Weiss, N. O. 2008, ApJ, 679, L57
  • Isobe et al. (2007) Isobe, H., Tripathi, D., & Archontis, V. 2007, ApJ, 657, L53
  • Jouve & Brun (2009) Jouve, L., & Brun, A. S. 2009, ApJ, 701, 1300
  • Kitiashvili et al. (2010) Kitiashvili, I. N., Bellot Rubio, L. R., Kosovichev, A. G., Mansour, N. N., Sainz Dalda, A., & Wray 2010, ApJ, 716, L181
  • Kosugi et al. (2007) Kosugi, T., Matsuzaki, K., Sakao, T., Shimizu, T., Sone, Y., Tachikawa, S., Hashimoto, T., Minesugi, K., Ohnishi, A., Yamada, T., Tsuneta, S., Hara, H., Ichimoto, K., Suematsu, Y., Shimojo, M., Watanabe, T., Shimada, S., Davis, J. M., Hill, L. D., Owens, J. K., Title, A. M., Culhane, J. L., Harra, L. K., Doschek, G. A., & Golub, L. 2007, Sol. Phys., 243, 3
  • Krause & Rädler (1980) Krause, F., & Rädler, K. 1980, Mean-field magnetohydrodynamics and dynamo theory, ed. Goodman, L. J. & Love, R. N.
  • Kubo et al. (2003) Kubo, M., Shimizu, T., & Lites, B. W. 2003, ApJ, 595, 465
  • Kubo et al. (2007) Kubo, M., Shimizu, T., & Tsuneta, S. 2007, ApJ, 659, 812
  • Lites (2009) Lites, B. W. 2009, Space Science Reviews, 144, 197
  • Lites et al. (2001) Lites, B. W., Elmore, D. F., & Streander, K. V. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 236, Advanced Solar Polarimetry – Theory, Observation, and Instrumentation, ed. M. Sigwarth, 33
  • Lites et al. (2004) Lites, B. W., Scharmer, G. B., Berger, T. E., & Title, A. M. 2004, Sol. Phys., 221, 65
  • Lites et al. (1998) Lites, B. W., Skumanich, A., & Martinez Pillet, V. 1998, A&A, 333, 1053
  • MacTaggart & Hood (2009) MacTaggart, D., & Hood, A. W. 2009, A&A, 507, 995
  • Magara (2006) Magara, T. 2006, ApJ, 653, 1499
  • Manchester IV et al. (2004) Manchester IV, W., Gombosi, T., DeZeeuw, D., & Fan, Y. 2004, ApJ, 610, 588
  • Martínez González et al. (2007) Martínez González, M. J., Collados, M., Ruiz Cobo, B., & Solanki, S. K. 2007, A&A, 469, L39
  • Martínez González & Bellot Rubio (2009) Martínez González, & Bellot Rubio, L. R. 2009, ApJ, 700, 1391
  • Martínez-Sykora et al. (2008) Martínez-Sykora, J., Hansteen, V., & Carlsson, M. 2008, ApJ, 679, 871
  • Martínez-Sykora et al. (2009) Martínez-Sykora, J., Hansteen, V., & Carlsson, M. 2009, ApJ, 702, 129
  • Moreno-Insertis (1997) Moreno-Insertis, F. 1997, Memorie della Societa Astronomica Italiana, 68, 429
  • Moreno-Insertis et al. (1994) Moreno-Insertis, F., Caligari, P., & Schüssler, M. 1994, Sol. Phys., 153, 449
  • Nordlund (2004) Nordlund, A. 2004, unpublished, presentation at
  • Pariat et al. (2004) Pariat, E., Aulanier, G., Schmieder, B., Georgoulis, M. K., Rust, D. M., & Bernasconi, P. N. 2004, ApJ, 614, 1099
  • Pariat et al. (2009) Pariat, E., Masson, S., & Aulanier, G. 2009, ApJ, 701, 1911
  • Parker (1955) Parker, E. N. 1955, ApJ, 121, 491
  • Parker (1994) —. 1994, ApJ, 433, 867
  • Rädler (1980) Rädler, K. 1980, Astronomische Nachrichten, 301, 101
  • Rempel et al. (2009a) Rempel, M., Schüssler, M., Cameron, R. H., & Knölker, M. 2009a, Science, 325, 171
  • Rempel et al. (2009b) Rempel, M., Schüssler, M., & Knölker, M. 2009b, ApJ, 691, 640
  • Rimmele (1997) Rimmele, T. R. 1997, ApJ, 490, 458
  • Roberts & Webb (1978) Roberts, B., & Webb, A. R. 1978, Sol. Phys., 56, 5
  • Sainz Dalda & Bellot Rubio (2008) Sainz Dalda, A., & Bellot Rubio, L. R. 2008, A&A, 481, L21
  • Schlichenmaier et al. (2010) Schlichenmaier, R., Bello Gonzalez, N., Rezaei, R., & Waldmann, T. A. 2010, A&A, 512, L1
  • Schrijver & Zwaan (2000) Schrijver, C. J., & Zwaan, C. 2000, Solar and Stellar Magnetic Activity, ed. Schrijver, C. J. & Zwaan, C.
  • Schüssler & Vögler (2006) Schüssler, M., & Vögler, A. 2006, ApJ, 641, L73
  • Shibata et al. (1989) Shibata, K., Tajima, T., Steinolfson, R. S., & Matsumoto, R. 1989, ApJ, 345, 584
  • Spruit (1981) Spruit, H. C. 1981, A&A, 98, 155
  • Spruit (1981) Spruit, H. C. 1981, The Physics of Sunspots, 98
  • Spruit et al. (1987) Spruit, H. C., Title, A. M., & van Ballegooijen, A. A. 1987, Sol. Phys., 110, 115
  • Stein et al. (2010) Stein, R. F., Lagerfjärd, A., Nordlund, Å., & Georgobiani, D. 2010, Sol. Phys., 34
  • Stein & Nordlund (2006) Stein, R. F., & Nordlund, Å. 2006, ApJ, 642, 1246
  • Strous et al. (1996) Strous, L. H., Scharmer, G., Tarbell, T. D., Title, A. M., & Zwaan, C. 1996, A&A, 306, 947
  • Strous & Zwaan (1999) Strous, L. H., & Zwaan, C. 1999, ApJ, 527, 435
  • Tobias et al. (1998) Tobias, S. M., Brummell, N. H., Clune, T. L., & Toomre, J. 1998, ApJ, 502, L177
  • Toriumi & Yokoyama (2010) Toriumi, S., & Yokoyama, T. 2010, ApJ, 714, 505
  • Tortosa-Andreu & Moreno-Insertis (2009) Tortosa-Andreu, A., & Moreno-Insertis, F. 2009, A&A, 507, 949
  • Tsuneta et al. (2008) Tsuneta, S., Ichimoto, K., Katsukawa, Y., Nagata, S., Otsubo, M., Shimizu, T., Suematsu, Y., Nakagiri, M., Noguchi, M., Tarbell, T., Title, A., Shine, R., Rosenberg, W., Hoffmann, C., Jurcevich, B., Kushner, G., Levay, M., Lites, B., Elmore, D., Matsushita, T., Kawaguchi, N., Saito, H., Mikami, I., Hill, L. D., & Owens, J. K. 2008, Sol. Phys., 74
  • Vögler et al. (2005) Vögler, A., Shelyag, S., Schüssler, M., Cattaneo, F., Emonet, T., & Linde, T. 2005, A&A, 429, 335
  • Weiss (1966) Weiss, N. O. 1966, Royal Society of London Proceedings Series A, 293, 310
  • Yelles Chaouche et al. (2009) Yelles Chaouche, L., Cheung, M. C. M., Solanki, S. K., Schüssler, M., & Lagg, A. 2009, A&A, 507, L53
  • Zwaan (1978) Zwaan, C. 1978, Sol. Phys., 60, 213
  • Zwaan (1987) —. 1987, ARA&A, 25, 83
Figure 1: Time sequence of continuum intensity images at nm (left) and synthetic longitudinal magnetograms (sampled at , right). The images show the full horizontal extent ( Mm) of the simulation domain.
Figure 2: The upper panel shows a synthetic magnetogram (sampled at ) of the simulated emerging flux region at hrs. The lower panel shows the vertical cross-section of the magnetic field strength () along (delineated by red dashed line).
Figure 3: Same as Fig 2 but at hrs. The two polarities of the active region at this time are much more compact.
Figure 4: Magnetic structure of the (semi)torus as it erupts onto the solar surface ( hrs). The same set of field lines are shown in both panels. The left panel shows the distribution of the vertical component of the magnetic field (scaled between kG) at Mm. At this depth the opposite polarities are in the form of two coherent roots. The panel on the right shows the corresponding magnetic map at (scaled between kG), which highlights the serpentine nature of the magnetic field lines near the surface as a consequence of the interaction between the emerging field and the granular convective flows.
Figure 5: Velocity field below the surface ( Mm) as the emerging magnetic field is traversing this layer ( hrs). Red and blue colors indicate down- and upflows, respectively, and the arrows indicate the horizontal velocity. The outflows at the periphery of the emerging flux region reaches speeds of up to km s whereas upflow velocities are closer to km s.
Figure 6: Joint probability density function of the horizontal field strength. The values were sampled at the vertical plane between and hours. The solid and dashed lines respectively show the scaling relations and with kG and g cm.
Figure 7: The emerging flux tube leads to a enhancement of pressure which drives outflow away from the emergence site. The three panels in this figure show, respectively, the time evolution of the surface (photospheric base, ) magnetic field strength, relative gas pressure perturbation, and component of velocity averaged over the band Mm. In all three panels, the green contours indicate enhancements of the gas pressure (relative to dyne cm) by and (for the purpose of having fewer contours, the pressure values have been smoothed in time with a Gaussian filter with min).
Figure 8: Mechanism for removal of mass and unsigned flux from the surface. (a) Schematic representation of how mass is removed from emerging magnetic field lines in a 2D scenario. In addition to undulating field lines, convective flows expel emerged flux (indicated by ovals labeled with positive and negative signs) from granular upflows. (b) 3D rendering of near-surface field lines in the simulated emerging flux region. Field lines in the upper panel are colored in accordance with the local density perturbation (about horizontal mean) with dark blue indicating density enhancement. Field lines in the lower panel are colored according to the vertical component of the momentum with red indicating downflowing material.
Figure 9: Poynting fluxes of magnetic energy through horizontal planes at and Mm. In both plots, the dashed curves and diamonds indicate the contributions by the first (bodily emergence or submergence of horizontal fields) and second (shearing by horizontal flows) terms in Eq. (9), respectively. The solid curves indicate the sum of the two contributions.
Figure 10: The top panel shows the azimuthally- and temporally-averaged (between and hr, sampled at ) magnetic field strengths (solid and dashed lines indicate the vertical and radial components, respectively) as functions of the radial distance from the axis of the positive-polarity spot in the simulation. The lower panels shows the corresponding quantities as defined in Eq. (14): (dotted), (dashed) and their sum (solid). Positive values of indicate an increase of flux (within the area ) with time.
Figure 11: The panels above show the azimuthally-, temporally- and spatially- (over a height-range of km about ) averaged profiles of the radial components of the magnetic pressure gradient force (), gas pressure gradient force (), magnetic tension (), net force and velocity () for gas with positive polarity field (, thick solid lines) and negative polarity field (, thick dashed lines). The developing spot centered at has positive polarity.
Figure 12: Structure of a light bridge - The left panel shows the continuum brightness of one of the spots in the simulated active region ( hrs). The formation of this light bridge is due to buoyant material entrained between adjacent regions of high magnetic field strength. Vertical cross-sections of the vertical velocity and magnetic field strength are shown for three cuts across segments of the light bridge of varying width. The green lines in the cross-sectional panels indicate the surface.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description