Dynamodriven plasmoid ejections above a spherical surface^{⋆}
^{1}
Nordita, AlbaNova University Center, Roslagstullsbacken 23, 10691 Stockholm, Sweden
email: joern@nordita.org
^{2}
Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden
Received: 4 April 2011
Accepted: 26 July 2011
Aims. We extend earlier models of turbulent dynamos with an upper, nearly forcefree exterior to spherical geometry, and study how flux emerges from lower layers to the upper ones without being driven by magnetic buoyancy. We also study how this affects the possibility of plasmoid ejection.
Methods. A spherical wedge is used that includes northern and southern hemispheres up to midlatitudes and a certain range in longitude of the Sun. In radius, we cover both the region that corresponds to the convection zone in the Sun and the immediate exterior up to twice the radius of the Sun. Turbulence is driven with a helical forcing function in the interior, where the sign changes at the equator between the two hemispheres.
Results. An oscillatory largescale dynamo with equatorward migration is found to operate in the turbulence zone. Plasmoid ejections occur in regular intervals, similar to what is seen in earlier Cartesian models. These plasmoid ejections are tentatively associated with coronal mass ejections (CMEs). The magnetic helicity is found to change sign outside the turbulence zone, which is in agreement with recent findings for the solar wind.
Key words: magnetohydrodynamics (MHD) / turbulence / Sun: dynamo / Sun: coronal mass ejections (CMEs) / stars: magnetic field
Movie is available in electronic form at http://www.aanda.org
© ESO, 2011
1. Introduction
Observations show that the Sun sheds mass through twisted magnetic flux configurations (Démoulin et al. 2002). Remark able examples of such helical ejections can be seen in the movies produced by the SOHO and SDO missions^{1}. Such events may be important for the solar dynamo (Blackman & Brandenburg 2003). They are generally referred to as coronal mass ejections (CMEs). Conventionally, CMEs are modeled by adopting a given distribution of magnetic flux at the solar surface and letting it evolve by shearing and twisting the magnetic field at its footpoints at the surface (Antiochos et al. 1999; Török & Kliem 2003). This approach is also used to model coronal heating (Gudiksen & Nordlund 2005; Bingert & Peter 2011). The success of this method depends crucially on the ability to synthesize the velocity and magnetic field patterns at the surface. Of course, ultimately such velocity and magnetic field patterns must come from a realistic simulation of the Sun’s convection zone, where the field is generated by dynamo action. In other words, we need a unified treatment of the convection zone and the CMEs. The difficulty here is the large range of time scales, from the 11year dynamo cycle to the time scales of hours and even minutes on which CMEs develop. Such a large range of time scales is related to the strong density stratification in the Sun, as can be seen from the following argument. In the bulk of the convection zone, the dynamo is controlled by rather slow motions with turnover times of days and months. The typical velocity depends on the convective flux via , where is the mean density and u_{rms} is the rms velocity of the turbulent convection. The dynamo cycle time can even be several hundred times the turnover time. In the corona, on the other hand, the typical time scale depends on the Alfvén time, L/v_{A}, where L is the typical scale of magnetic structures and is the Alfvén speed for a given magnetic field strength B. Here, μ_{0} is the vacuum permeability.
In a recent paper, Warnecke & Brandenburg (2010) attempted a new approach of a unified treatment by combining a dynamogenerated field in the convection zone with a nearly forcefree coronal part, albeit in a local Cartesian geometry. In this paper, we go a step further by performing direct numerical simulations (DNS) in spherical geometry. We also include density stratification due to gravity, but with a density contrast between the dynamo interior and the outer parts of the simulation domain that is much less (about 20) than in the Sun (around 14 orders of magnitude). This low density contrast is achieved by using an isothermal configuration with constant sound speed c_{s}. Hence, the average density depends only on the gravitational potential and is given by , where G is Newton’s gravitational constant, M is the central mass, and r is the distance from the center. As convection is not possible in such an isothermal setup, we drive turbulence by an imposed helical forcing that vanishes outside the convection zone. This also helps achieving a strong largescale magnetic field. The helicity of the forcing is negative (positive) in the northern (southern) hemisphere and smoothly changes sign across the equator. Such a forcing gives rise to an α^{2} dynamo with periodic oscillations and equatorward migration of magnetic activity (Mitra et al. 2010a). We ignore differential rotation, so there is no systematic shearing in latitude. The only twisting comes then from the same motions that also sustain the dynamogenerated magnetic field. Note that in our model the mechanism of transport of the magnetic field to the surface is not magnetic buoyancy. Instead, we expect that, twisted magnetic fields will expel themselves to the outer regions by the Lorentz force.
Our aim in this paper is not to provide a model as close to reality as possible, but to show that it is possible to capture the phenomenon of CMEs (or, more generally, plasmoid ejections) within a minimalistic model that treats the convection zone and the outer parts of the Sun in a selfconsistent manner. That is, the magnetic field in the convection zone is dynamically generated by dynamo action and the motions are not prescribed by hand, but they emerge as a solution of the momentum equation and include magnetic backreaction from the Lorentz force.
Given that gravity decreases with radius, there is in principle the possibility of a radial wind with a critical point at (Choudhuri 1998). However, as we use stressfree boundary conditions with no mass flux in the radial direction, no such wind can be generated in our simulations. Nevertheless, we observe radially outward propagation of helical magnetic field structures without mass flux. Furthermore, our results for the flux of magnetic helicity compare well with recent measurements of the same in the solar wind (Brandenburg et al. 2011). Our approach might therefore provide new insights not only for CMEs and dynamo theory, but also for solar wind turbulence.
2. The model
We use spherical polar coordinates, (r,θ,φ). As in earlier work of Mitra et al. (2009, 2010a), our simulation domain is a spherical wedge. The results of Mitra et al. (2009) for such a wedge are consistent with those of Livermore et al. (2010) for a full spherical shell, both of which ignored the effects of an equator, which was included in the work of Mitra et al. (2010a). Our model is a bilayer in the radial direction. The inner layer is forced with random helical forcing functions which have different signs of helicity in the two hemispheres. This models the helical aspects of convection in the Sun. We shall often call the inner layer “turbulence zone”. The radius separating the two layers corresponds to the solar radius, r = R. This length scale is used as our unit of length. The inner layer models some aspects of the convection zone (0.7R ≤ r ≤ R) without however having any real convection, and the outer layer (R ≤ r ≤ 2R) models the solar corona. We consider the range π/3 ≤ θ ≤ 2π/3, corresponding to ± 30° latitude, and 0 < φ < 0.3, corresponding to a longitudinal extent of 17°. Here, θ is the polar angle and φ the azimuth. At the solar surface at R = 700 Mm, this would correspond to an area of about 730 × 210 Mm^{2}, which could encompass several active regions in the Sun. In our model the momentum equation is (1)where is the viscous force, ν is the kinematic viscosity, is the traceless rateofstrain tensor, semicolons denote covariant differentiation, is the specific pseudoenthalpy, c_{s} = const. is the isothermal sound speed, and g = −GMr/r^{3} is the gravitational acceleration. We choose , so r_{∗} = 1.5R lies within our domain. This value is rather close to the surface and would lead to significant mass loss if there was a wind, but this is suppressed by using impenetrative outer boundaries.
The forcing function is given as the product of two parts, (2)where is a profile function connecting the two layers and w is the width of the transition at the border between the two layers (r = R). In other words, we choose the external force to be zero in the outer layer, r > R. The function f consists of random plane helical transverse waves with relative helicity and wavenumbers that lie in a band around an average forcing wavenumber of k_{f}R ≈ 63. This value should be compared with the normalized wavenumber k_{1}R, corresponding to the thickness of the shell ΔR, which yields k_{1}R = 2πR/ΔR ≈ 21, so the effective scale separation ratio, k_{f}/k_{1}, is about 3.
In Eq. (2) the last argument of f(r,θ,φ,t;σ) denotes a parametric dependence on the helicity which is here chosen to be σ = −cosθ such that the kinetic helicity of the turbulence is negative in the northern hemisphere and positive in the southern. Specifically, f is given by (Haugen et al. 2003) (3)where A_{f} is a nondimensional forcing amplitude, and x is the position vector. The wavevector k(t) and the random phase − π < φ(t) ≤ π change at every time step, so f(x,t) is δcorrelated in time. The normalization factor N is chosen on dimensional grounds to be N = c_{s}(kc_{s}/δt)^{1/2}. At each timestep we select randomly one of many possible wavevectors in a certain range around a given forcing wavenumber. The average wavenumber is referred to as k_{f}. We ignore curvature effects in the expression for the forcing function and thus force the system with what would correspond to transverse helical waves in a Cartesian coordinate system, i.e., (4)where − 1 ≤ σ ≤ 1 is the helicity parameter of the forcing function, (5)is a nonhelical forcing function, and is an arbitrary unit vector not aligned with k; note that .
The pseudoenthalpy term in Eq. (1) emerges from the fact that for an isothermal equation of state the pressure is given by , so the pressure gradient force is given by . The continuity equation is then written in terms of h as (6)Equations (1) and (6) are solved together with the uncurled induction equation for the vector potential A in the resistive gauge (Candelaresi et al. 2011), (7)where η is the (constant) magnetic diffusivity, so the magnetic field is given by B = ∇ × A and thus obeys ∇·B = 0 at all times. The gauge can in principle become important when calculating the magnetic helicity density A·B, although the part resulting from the smallscale fields is expected to be independent of the gauge (Subramanian & Brandenburg 2006; Hubbard & Brandenburg 2010), while that of the largescale fields is not.
Our wedge is periodic in the azimuthal direction. For the velocity, we use stressfree boundary conditions on all other boundaries. For the magnetic field we employ vertical field conditions on r = 2R and perfect conductor conditions on both r = 0.7R and the two θ boundaries. Time is expressed in units of , which is the eddy turnover time in the turbulence zone, and u_{rms} is the rms velocity in r < R. Density is given in units of the mean density in the turbulence zone, . The magnetic field is expressed in units of the mean equipartition value, B_{eq}, defined via . We use a magnetic diffusivity that is constant in space and time and its value is given in terms of the magnetic Reynolds number, defined as (8)where k_{f} is the characteristic scale of the external force, defined above. This also turns out to be the energy containing scale of the fluid. In the following analysis, we use φ averages, defined as =^{∫}F(r,θ,φ,t)dφ/2π. Occasionally we also use time averages denoted by ⟨ . ⟩ _{t}. We perform DNS with the Pencil Code^{2}, which is a modular highorder code (sixth order in space and thirdorder in time, by default) for solving a large range of partial differential equations, including the ones relevant in the present context.
3. Results
3.1. Dynamo in the turbulence zone
We start our DNS with seed magnetic field everywhere in the domain. Owing to the helical forcing in the turbulent layer, a largescale magnetic field is produced by dynamo action. The dynamo is cyclic with equatorward migration of magnetic fields. This dynamo was studied by DNS in Mitra et al. (2010a) and has been interpreted as an α^{2} dynamo. The possibility of oscillating α^{2} dynamos was known since the early papers of Baryshnikova & Shukurov (1987) and Rädler & Bräuer (1987), who showed that a necessary condition for oscillations is that the α effect must change sign in the domain.
Fig. 1 Initial exponential growth and subsequent saturation behavior of the magnetic field in the interior for forced turbulence with dynamo action for Run A. The magnetic field strength is oscillating with twice the dynamo frequency 2ω_{cyc}. 

Open with DEXTER 
Fig. 2 Equatorward migration, as seen in visualizations of B_{r} for Run D at r = R over a horizontal extent Δθ = 58° and Δφ = 17°. 

Open with DEXTER 
Fig. 3 Periodic variation of and in the turbulence zone. Dark blue stands for negative and light yellow for positive values. The dotted horizontal lines show the location of the equator at θ = π/2. The magnetic field is normalized by the equipartition value. Taken from Run A. 

Open with DEXTER 
Summary of runs discussed in this paper.
The dynamo first grows exponentially and then saturates after around 300 turnover times, see Fig. 1. After saturation the dynamo produces a largescale magnetic field with opposite polarities in the northern and southern hemispheres. In Fig. 2 we plot the radial magnetic field at the surface of the dynamo region at r = R. This layer would correspond to the solar photosphere if we had a more realistic solar model, which would include higher stratification as well as cooling and radiation effects. The six wedges represent different times and show clearly an equatorward migration of the radial magnetic field with polarity reversal every half cycle. The other components of the magnetic field (not plotted) also shows the same behavior. Comparing the first (t/τ = 3028) and the last (t/τ = 3101) panel, the polarity has changed sign in a time interval Δt/τ ≈ 100. The oscillatory and migratory properties of the dynamo is also seen in the butterfly diagram of Fig. 3 for and . In Fig. 1 one can also verify that the oscillation period is around 200 turnover times, corresponding to a nondimensional dynamo cycle frequency of τω_{cyc} = 0.032 and the field strength in the turbulent layer varies between 1.2 and 1.6 of the equipartition field strength. This value of the cycle frequency is roughly consistent with an estimate of Mitra et al. (2010b) that , where k_{m} is the relevant wavenumber of the mean field. Using η_{t} ≈ η_{t0} ≡ u_{rms}/3k_{f} (Sur et al. 2008), we find τω_{cyc} ≈ 0.2(k_{m}/k_{f})^{2} ≈ 0.02, where we have assumed k_{m} ≈ 2π/0.3R ≈ 20k_{1} and k_{f} ≈ 60k_{1}. The estimate of Mitra et al. (2010b) applies to perfectly conducting outer boundary conditions, which might explain the remaining discrepancy.
A summary of all runs is given in Table 1, where the amplitudes of the magnetic field show a weak nonmonotonous dependence on the magnetic Reynolds number R_{m}. For larger values of R_{m}, the magnetic field strength decreases slightly with increasing R_{m}, but it is weaker than in some earlier α^{2} dynamos with open boundaries (Brandenburg & Dobler 2001). This could be due to two reasons. Firstly, our simulations are far from the asymptotic limit of large magnetic Reynolds numbers, in which the results of Brandenburg & Dobler (2001) are applicable. The maximum value R_{m} is in our simulations approximately 15 times the critical R_{m}. The second reason could be that we have expulsion of magnetic helicity from our domain which was not present in (Brandenburg & Dobler 2001). We find the peak of the R_{m} dependency at R_{m} = 9, corresponding to Run C. The dynamo cycle frequency shows a weak decrease (by a factor of 2) as the magnetic Reynolds number increases (by a factor of 30).
3.2. Phase relation between radial and azimuthal fields
Although our dynamo model does not include important features of the Sun such as differential rotation, some comparison may still be appropriate. For the Sun, one measures the mean radial field by averaging the lineofsight magnetic field from synoptic magnetograms. The azimuthal field is not directly observed, but its sign can normally be read off by looking at the magnetic field orientation of sunspot pairs. Existing data suggest that radial and azimuthal fields are approximately in outofphase (Yoshimura 1976). This is reasonably well reproduced by αΩ dynamos models, where the radial field lags behind the azimuthal one by 0.75π (Stix 1976). However, in the present work, radial and azimuthal fields are approximately in phase with a phase difference of 0.3π inside the dynamo region; see Fig. 4. Future studies will include the nearsurface shear layer, which has been suspected to play an important role in producing equatorward migration (Brandenburg 2005). This would also help reproducing the observed phase relation.
Fig. 4 Time evolution of the radial magnetic field B_{r} (solid line) and the azimuthal magnetic field B_{φ} (dotted line) in the dynamo region averaged over the radius r and azimuth at θ = ± 7°. To improve the statistics, we calculate the components of the magnetic field as the antisymmetric part in latitude, i.e., for i = r,φ. 

Open with DEXTER 
3.3. Relation between kinetic and magnetic energies
Next we investigate the relation between the rms values of the magnetic field and the velocity. Both quantities are oscillating in time with a typical period of 200 turnover times. In Fig. 5 we compare the time evolution of the magnetic field strength and the rms velocity. The magnetic field is calculated in the dynamo region and normalized to the thermal equipartition field strength. The phase difference between the two is slightly less than π within the dynamo region. This basically shows that the magnetic field quenches the turbulence.
Fig. 5 Phase relations between the magnetic field and the velocity in the dynamo region. The magnetic field is plotted as , normalized with the equipartition field of the sound speed, (=) as a solid and black line. The rms velocity, normalized by the sound speed c_{s}, is plotted as a dashed red line, and has been smoothed over 5 neighboring data points to make it more legible. Taken from Run A. 

Open with DEXTER 
3.4. Density variations
The density is stratified in radius and varies by over an order of magnitude. For all the runs listed in Table 1 the density fluctuates about the hydrostatic equilibrium value, . The relative fluctuations are of comparable strength at all radii; see Fig. 6.
Fig. 6 Radial dependence of density overplotted at different times. In the inset is the linear behavior of the logarithmic density log ρ/ρ_{0} to the inverse of the radius shown. Taken from Run A. 

Open with DEXTER 
3.5. Magnetic field outside the turbulence zone
The magnetic field averaged over the entire domain is more than 5 times smaller than in the turbulence zone. In Fig. 7 we show that the magnetic field is concentrated to the turbulence zone and drops approximately exponentially with a scale height of about 0.23R in the outer parts for r > R. The toroidal component of the magnetic field is dominant in the turbulence layer, but does not play a significant role in the outer part. By contrast, the radial field is weak in the inner parts and dominates in the outer.
Magnetic structures emerge through the surface and create field line concentrations that reconnect, separate, and rise to the outer boundary of the simulation domain. This dynamical evolution is seen in a sequence of field line images in Fig. 8, where field lines of the mean field are shown as contours of and colors represent . Prior to a plasmoid ejection we see a convergence of antiparallel radial field lines, which then reconnect such that the newly reconnected field lines move away from the reconnection site. The actual reconnection seems to happen much faster than the subsequent ejection.
Fig. 7 Radial dependence of the mean squared magnetic field, (solid line), compared with those of B_{r} (dotted), B_{θ} (dashed), and B_{φ} (dashdotted). All quantities are averaged over 13 dynamo cycles. The inset shows the same quantities in a logarithmic representation. Taken from Run A. 

Open with DEXTER 
Fig. 8 Time series of formation of plasmoid ejections in spherical coordinates. Contours of are shown together with a colorscale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the rθ plane. The dashed horizontal lines show the location of the surface at r = R. Taken from Run D. 

Open with DEXTER 
In the outer layers, the magnetic field emerges as large structures that correlate with reconnection events of magnetic fields. In Rust (1994) such phenomena have been described as magnetic clouds. We find recurrent ejections of magnetic field lines with concentrations and reconnection events, but the occurrence of structures such as magnetic clouds does not happen completely regularly, i.e., these structured events would be difficult to predict.
In Fig. 9 we show a closeup of the magnetic field. A configuration resembling a reconnection event is clearly seen. Here, the contours represent magnetic field lines with solid and dashed lines denoting counterclockwise and clockwise orientations, respectively. The solid antiparallel field lines reconnect around r = 0.9R and separate to the left and to the right. On the righthand side, a plasmoid has formed, which is eventually ejected. This plasmoid appears as a CMElike ejection in the first panel of Fig. 11. These plasmoids, as seen more clearly in Figs. 8 and 10, appear as a concentration of field lines that propagate outwards. The fact that reconnection happens predominantly in the upper parts of the turbulence zone suggests that turbulence is needed to enable fast reconnection (Lazarian & Vishniac 1999).
Additionally, we find reconnection as a result of the interaction between ejections. As plotted in Fig. 10, the field lines of two subsequent events have the opposite field line direction, which can then interact in the outer layers. Comparison with the first panel of Fig. 11 shows that the current helicity has a correlation with the separatrices of the two polarities of the field lines. We also find that in the interaction region the field lines have high density and are more strongly concentrated.
3.6. Current helicity
The current helicity (J·B) is often used as a useful proxy for the magnetic helicity (A·B) at small scales, because, unlike magnetic helicity, it is gaugeindependent. Current helicity has also been observed in the Sun (Seehafer 1990) and it has been obtained from meanfield dynamo models (Dikpati & Gilman 2001). In the present paper we are particularly interested in the current helicity outside the Sun. We normalize it by the rdependent timeaveraged mean squared field to compensate for the radial decrease of . In Fig. 12, we have also averaged in latitude from 20° to 28°. In the turbulence zone the sign of is the same as that of kinetic helicity which, in turn, has the same sign as the helicity of the external forcing, i.e. of σ; see Figs. 11 and 12.
However, to our surprise, above the surface, and separately for each hemisphere, the signs of current helicity tend to be opposite to those in the turbulence zone; see Fig. 11 for the panels of t/τ = 1669 and t/τ = 1740. To demonstrate that plasmoid ejections are recurrent phenomena, we look at the evolution of as a function of t and r. This is done in Fig. 12 for Run A. It turns out that the speed of plasmoid ejecta is about 0.45 of the rms velocity of the turbulence in the interior region for Reynolds numbers up to 15 and about 0.3 up to 30. To compare with the Sun, we estimate the rms velocity of the turbulence in terms of the convective energy flux via . The density of the corona is ρ_{cor} ≈ 10^{12} kg m^{3}, so our estimate would suggest V_{ej} ≈ 0.3(F/ρ_{cor})^{1/3} ≈ 1200 km s^{1}, which is somewhat above the observed speeds of 400–1000 km s^{1}. The time interval between subsequent ejections is around 100 τ for Run A. As seen from Table 1, the ejection interval is independent of the forcing amplitude A_{f} and increases weakly with magnetic Reynolds number, but it seems to be still comparable to half the dynamo cycle period, i.e., Δτ ≈ T_{cyc}/2. This means that plasmoid ejections happen about twice each cycle. It is therefore not clear whether such a result can be extrapolated to the real Sun.
Fig. 9 Time series of a reconnection event in an Xpoint as a closeup view. Contours of are shown together with a colorscale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the rθ plane, where solid lines represent counterclockwise magnetic field lines and the dash ones clockwise. The dashed vertical lines show the location of the surface at r = R. Taken from Run D. 

Open with DEXTER 
In our simulations, we find the ejections to have the shape of the characteristic threepart structure that is observed in real CMEs. This is particularly clear in Fig. 11, where the ejections seem to contain three different parts. In the center we find a ballshaped structure consisting of one polarity of current helicity, where at the front of the ejection a bow of opposite polarity had formed. In between these two structures the current helicity is close to zero and appears as a cavity. These three parts (prominence, cavity, and front) are described by Low (1996) for CMEs and are generally referred to as “threepart structure”. The basic shape of the ejection is independent of the used forcing amplitudes and the kinetic and magnetic Reynolds numbers. It should, however, be noted that the threepart structure of the ejections becomes clearer at magnetic higher Reynolds numbers (e.g., for Runs D and G compared with Run A, for example). In the Sun, the plasma is confined to loops of magnetic field with flows along field lines due to the low plasma beta in the solar corona. This is also seen in our simulations displayed in Fig. 11, where the ejections follow field lines and appear to create looplike structures. An animation of the detailed time evolution of the CMElike structures emerging recurrently into the solar corona is available in the online edition (see Fig. 11)^{3}. However, since our choice of boundary conditions does not allow mass flux at the outer boundary, no plasma can actually leave the domain. The recurrent nature of the plasmoid ejections found here and in Warnecke & Brandenburg (2010) is not yet well understood. In contrast to Warnecke & Brandenburg (2010), where there are no strong oscillations present, here the ejections seem to correlate with the dynamo cycle. In each hemisphere of the turbulence zone a magnetic field is created with different polarity. After they have migrated to the equator, they merge and produce an ejection. However, comparing with the results of Warnecke & Brandenburg (2010), which are similar to those in the present paper, the oscillation cannot be the only explanation for the recurrence of the ejections. As we have seen in Fig. 12, these events export magnetic helicity out of the domain. For the dynamo, on the other hand, magnetic helicity losses play a role only in the nonlinear stage. It is therefore conceivable that the regular occurrence of plasmoid ejections is connected with nonlinear relaxation oscillations rather than with the dynamo cycle which is essentially a linear phenomenon. This is also suggested by the results of Warnecke & Brandenburg (2010), where recurrent ejections occur without oscillations of the largescale field.
Fig. 10 Magnetic field configuration at the time of a ejection. Contours of are shown together with a colorscale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the rθ plane, where the solid lines represent counterclockwise magnetic field lines and the dash ones clockwise. The dashed vertical lines show the location of the surface at r = R. Taken from Run D. 

Open with DEXTER 
From Figs. 11 and 12 we conclude that in each hemisphere the sign of current helicity outside the turbulence zone is mostly opposite to that inside the turbulence zone. A stronger trend is shown in the cumulative mean of current helicity over time. This is shown in Fig. 13, where we plot the time evolution of the φ averaged current helicity at r = 1.5 R and 28° latitude, which is a safe distance away from the outer r and θ boundaries so as not to perturb our results, which should thus give a reasonable representation of the outer layers. For the northern hemisphere the current helicity (solid black line) and the cumulative mean (solid red line) show positive values and for the southern hemisphere (dotted lines) negative values. This agrees with results of Brandenburg et al. (2009), where the magnetic helicity of the field in the exterior has the opposite sign than in the interior.
Fig. 11 Time series of coronal ejections in spherical coordinates. The normalized current helicity, , is shown in a colorscale representation for different times; dark blue stands for negative and light yellow for positive values. The dashed horizontal lines show the location of the surface at r = R. Taken from Run D. The temporal evolution is shown in a movie available as online material. 

Open with DEXTER 
Fig. 12 Dependence of the dimensionless ratio on time t/τ and radius r in terms of the solar radius. The top panel shows a narrow band in θ in the northern hemisphere and the bottom one a thin band in the southern hemisphere. In both plots we have also averaged in latitude from 20° to 28°. Dark blue stands for negative and light yellow for positive values. The dotted horizontal lines show the location of the surface at r = R. 

Open with DEXTER 
Fig. 13 Dependence of the dimensionless ratio on time t/τ at radius r = 1.5 R and 28° latitude. The solid line stands for the northern hemisphere and the dotted for the southern hemisphere. The red lines represent the cumulative mean for each hemisphere. 

Open with DEXTER 
To investigate whether the sign of the current helicity in the turbulence zone is different from that and in the outer parts, we show in Fig. 14 the azimuthally and timeaveraged current helicity as a function of radius and colatitude. It turns out that, even though we have averaged the result over several thousand turnover times, the hemispheric sign rule of current helicity is still only approximately obeyed in the outer layers – even though it is nearly perfectly obeyed in the turbulence zone. There remains substantial uncertainty, especially near the equator. This suggests that meaningful statements about magnetic and current helicities in the solar wind can only be made after averaging over sufficiently long stretches of time.
Fig. 14 Current helicity averaged over 3900 turnover times. Legend is the same as in Fig. 11. Dark blue corresponds to negative values, while the light yellow corresponds to positive value. Taken from Run D. 

Open with DEXTER 
3.7. Magnetic helicity fluxes
In view of astrophysical dynamo theory it is important to understand the amount of magnetic helicity that can be exported from the system. Of particular interest here is the magnetic helicity associated with the smallscale magnetic field. Under the assumption of scale separation, this quantity is gaugeindependent (Subramanian & Brandenburg 2006), so we can express it in any gauge. This has been verified in simulations both with an equator (Mitra et al. 2010b) and without (Hubbard & Brandenburg 2010). Here, we compute the magnetic helicity flux associated with the smallscale field by subtracting that of the azimuthally averaged field from that of the total field, i.e., (9)where E = μ_{0}ηJ − U × B is the electric field. This is also the way how the magnetic helicity flux from the smallscale field was computed in Hubbard & Brandenburg (2010), where the magnetic helicity flux from the total and largescale fields was found to be gaugedependent, but that from the smallscale field was not. In Fig. 16 we compare the flux of magnetic helicity of the smallscale field across the outer surfaces in the northern and southern hemispheres with that through the equator. It turns out that a major part of the flux goes through the equator. The part of the magnetic helicity flux that goes through the surface is about 20% of . However, the magnetic helicity flux should primarily depend on rather than B_{eq}. In the present simulations, where the dynamo works with a fully helical field, the two are comparable. This is not the case in the Sun, where the estimated field strength is expected to be about 300 G (Brandenburg 2009). Thus, to compare with the Sun, a more reasonable guess for the magnetic helicity flux would be about 20% of . Integrating this over one hemisphere and multiplying this with the 11 year cycle time, we find the total magnetic helicity loss to be , which corresponds to 5 × 10^{47} Mx^{2} if we use η_{t} = 3 × 10^{12} cm^{2} s^{1}. This value exceeds the estimated upper limit for the solar dynamo of about 10^{46} Mx^{2} per cycle given by Brandenburg (2009). However the estimate by Brandenburg (2009) is based on an αΩ dynamo model with α effect and shear that yield a period comparable with the 11 year period of the Sun. Therefore, the discrepancy with the present model, where shear plays no role, should not be surprising. Instead, it tells us that a dynamo without shear has to shed even more magnetic helicity than one with shear.
The magnetic helicity flux of the largescale field may not a priori be gaugeinvariant. However, the system is in a statistically steady state and, in addition, the magnetic helicity integrated over each cycle is constant. In that case the divergence of the magnetic helicity flux is also gaugeinvariant. Furthermore, the shellintegrated magnetic helicity cannot have a rotational component and is therefore uniquely defined. In Fig. 15 we plot this flux and see that its maxima tend to occur about 50 turnover times after magnetic field maxima; see Figs. 4 and 5. The helicity flux of the smallscale field does not show a clear behavior. Since the ejections appear to be related to the magnetic field strength in this way, one might conclude that the magnetic helicity flux of the largescale field is transported through these ejections. This result is somewhat unexpected and deserves to be reexamined more thoroughly in future simulations where cycle and ejection frequencies are clearly different from each other.
Fig. 15 Time evolution of the magnetic helicity flux of the largescale field, smooth over two data points. Here, the mean of magnetic helicity flux out through the surface of the northern hemisphere (black) is shown, together with that through the southern hemisphere (dotted red), and the equator (dashed blue). 

Open with DEXTER 
Fig. 16 Cumulative mean of the time evolution of the magnetic helicity flux of the smallscale field, , normalized by , where η_{t} ≈ η_{t0} ≡ u_{rms}/3k_{f} was defined in Sect. 3.1. Here, the mean of magnetic helicity flux out through the surface of the northern hemisphere (black) is shown, together with that through the southern hemisphere (dotted red), and the equator (dashed blue). 

Open with DEXTER 
Fig. 17 Dependence of the latitudinal component of the magnetic helicity flux, , compared with the latitudinal gradient of the magnetic helicity density of the smallscale field, , at r/R = 0.85. The largescale variation of the latter agrees with that of the former if the gradient is multiplied by an effective diffusion coefficient for magnetic helicity of κ_{t} ≈ 3η_{t0}. 

Open with DEXTER 
Next, let us look at the magnetic helicity flux of the smallscale field. On earlier occasions, Mitra et al. (2010b) and Hubbard & Brandenburg (2010) have been able to describe the resulting magnetic helicity flux by a Fickian diffusion ansatz of the form , where κ_{h}/η_{t0} was found to be 0.3 and 0.1, respectively. In Fig. 17 we show that the present data allow a similar representation, although the uncertainty is large. It turns out that κ_{h}/η_{t0} is about 3, suggesting thus that turbulent magnetic helicity exchange across the equator can be rather efficient. Such an efficient transport of magnetic helicity out of the dynamo region is known to be beneficial for the dynamo in that it alleviates catastrophic quenching (Blackman & Brandenburg 2003). In this sense, the inclusion of CMElike phenomena is not only interesting in its own right, but it has important beneficial consequences for the dynamo itself in that it models a more realistic outer boundary condition.
3.8. Comparison with solar wind data
Our results suggest a reversal of the sign of magnetic helicity between the inner and outer parts of the computational domain. This is in fact in agreement with recent attempts to measure magnetic helicity in the solar wind (Brandenburg et al. 2011). They used the Taylor hypothesis to relate temporal fluctuations of the magnetic field to spatial variations by using the fact that the turbulence is swept past the space craft with the mean solar wind. This idea can in principle also be applied to the present simulations, provided we use the obtained mean ejection speed V_{ej} (see Table 1) for translating temporal variations (in t) into spatial ones (in r) via r = r_{0} − V_{ej}t. Under the assumption of homogeneity, one can then estimate the magnetic helicity spectrum as ; see Matthaeus et al. (1982) and Eq. (9) of Brandenburg et al. (2011). Here, hats indicate Fourier transforms and the asterisk denotes complex conjugation.
In Figs. 18 and 19 we show the results for the northern and southern hemispheres, as well as time series of the two relevant components B_{θ} and B_{φ}. The resulting magnetic helicity spectra, normalized by 2μ_{0}E_{M}/k, where E_{M} is the magnetic energy spectrum, give a quantity that is between − 1 and + 1. Note that the time traces are governed by a low frequency component of fairly large amplitude. In addition, there are also other components of higher frequency, but they are harder to see. The results suggest positive magnetic helicity in the north and negative in the south, which would be indicative of the helicities of the solar wind at smaller length scale. It also agrees with the current helicities determined using explicit evaluation in real space. On the other hand, the Parker spiral (Parker 1958) might be responsible for the magnetic helicity at large scales (Bieber et al. 1987a,b).
Fig. 18 Helicity in the northern outer atmosphere. The values are written out at the point, r = 1.5 R, 90° − θ = 17°, and φ = 9°. Top panel: phase relation between the toroidal B_{φ} and poloidal B_{θ} field, plotted over time t/τ. Bottom panel: helicity H(k) is plotted over normalized wave number kR. The helicity is calculated with the Taylor hypothesis using the Fourier transformation of the poloidal and toroidal field. 

Open with DEXTER 
Fig. 19 Helicity in the southern outer atmosphere. The values are written out at the point, r = 1.5 R, 90° − θ = −17° and φ = 8.6°. Top panel: phase relation between the toroidal B_{φ} and poloidal B_{θ} field, plotted over time t/τ. Bottom panel: helicity H(k) is plotted over normalized wave number kR. The helicity is calculated with the Taylor hypothesis using the Fourier transformation of the poloidal and toroidal field. 

Open with DEXTER 
4. Conclusions
In the present work we have demonstrated that CMElike phenomena are ubiquitous in simulations that include both a helicitydriven dynamo and a nearly forcefree exterior above it. This was first shown in Cartesian geometry (Warnecke & Brandenburg 2010) and is now also verified for spherical geometry. A feature common to both models is that the helical driving is confined to what we call the turbulence zone, which would correspond to the convection zone in the Sun. In contrast to the earlier work, we have now used a helical forcing for which the kinetic helicity changes sign across the equator. This makes the dynamo oscillatory and displays equatorward migration of magnetic field (Mitra et al. 2010a). More importantly, unlike our earlier work where the gas pressure was neglected in the outer parts, it is fully retained here, because it does automatically become small away from the surface due to the effect of gravity that is here included too, but was neglected in the earlier Cartesian model. The solutions shown here and those of Warnecke & Brandenburg (2010) demonstrate that this new approach of combining a selfconsistent dynamo with a coronalike exterior is a viable one and can model successfully features that are similar to those in the Sun. However, our model is still not sophisticated enough for direct quantitative comparisons.
Of particular interest is the sign change of magnetic and current helicities with radius. Although similar behavior has also been seen in other Cartesian models of Brandenburg et al. (2009), its relevance for the Sun was unknown until evidence for similar sign properties emerged from solar wind data (Brandenburg et al. 2011). In the present case we were also able to corroborate similar findings by using the Taylor hypothesis based on the plasmoid ejection speed. It is remarkable that this appears to be sufficient for relating spatial and temporal fluctuations to each other.
There are many ways in which the present model can be extended and made more realistic. On the one hand, the assumption of isothermal stratification could be relaxed and the increase of temperature in the corona together with the solar wind could be modeled in a reasonably realistic way. On the other hand, the dynamo model could be modified to include the effects of convection and of latitudinal differential rotation. Among other things, differential rotation would lead to the Parker spiral (Parker 1958), which is known to produce magnetic helicity of its own (Bieber et al. 1987a,b). It would then be interesting to see how this affects the magnetic helicity distribution seen in the present model.
The movie is also available at http://www.youtube.com/watch?v=aRPgxQyP24.
Acknowledgments
We thank the anonymous referee for useful suggestions and Matthias Rheinhardt for help with the implementation. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm, the National Supercomputer Centers in Linköping and the High Performance Computing Center North. This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952 and the Swedish Research Council Grant No. 62120074064, and the National Science Foundation under Grant No. NSF PHY0551164.
References
 Antiochos, S. K., De Vore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485 [NASA ADS] [CrossRef] (In the text)
 Baryshnikova, Y., & Shukurov, A. M. 1987, Astron. Nachr., 308, 89 [NASA ADS] [CrossRef] (In the text)
 Bieber, J. W., Evenson, P. A., & Matthaeus, W. H. 1987a, ApJ, 315, 700 [NASA ADS] [CrossRef] (In the text)
 Bieber, J. W., Evenson, P. A., & Matthaeus, W. H. 1987b, J. Geophys. Res., 14, 864 (In the text)
 Bingert, S., & Peter, H. 2011, A&A, 530, A112 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Blackman, E. G., & Brandenburg, A. 2003, ApJ, 584, L99 [NASA ADS] [CrossRef] (In the text)
 Brandenburg, A. 2005, ApJ, 625, 539 [NASA ADS] [CrossRef] (In the text)
 Brandenburg, A. 2009, Plasma Phys. Control. Fusion, 51, 124043 [NASA ADS] [CrossRef] (In the text)
 Brandenburg, A., & Dobler, W. 2001, A&A, 369, 329 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Brandenburg, A., Candelaresi, S., & Chatterjee, P. 2009, MNRAS, 398, 1414 [NASA ADS] [CrossRef] (In the text)
 Brandenburg, A., Subramanian, K., Balogh, A., & Goldstein, M. L. 2011, ApJ, 734, 9 [NASA ADS] [CrossRef] (In the text)
 Candelaresi, S., Hubbard, A., Brandenburg, A., & Mitra, D. 2011, Phys. Plasmas, 18, 012903 [NASA ADS] [CrossRef] (In the text)
 Choudhuri, A. R. 1998, The Physics of Fluids and Plasmas (Cambridge University Press) (In the text)
 Démoulin, P., Mandrini, C. H., van DrielGesztelyi, L., et al. 2002, ApJ, 382, 650 (In the text)
 Dikpati, M., & Gilman, P. A. 2001, ApJ, 559, 428 [NASA ADS] [CrossRef] (In the text)
 Gudiksen, B. V., & Nordlund, Å. 2005, ApJ, 618, 1031 [NASA ADS] [CrossRef] (In the text)
 Haugen, N. E. L., Brandenburg, A., & Dobler, W. 2004, Phys. Rev. E, 70, 016308 [NASA ADS] [CrossRef] (In the text)
 Hubbard, A., & Brandenburg, A. 2010, Geophys. Astrophys. Fluid Dyn., 104, 577 [NASA ADS] [CrossRef] (In the text)
 Käpylä, P. J., & Brandenburg, A. 2009, ApJ, 699, 1059 [NASA ADS] [CrossRef] (In the text)
 Krause, F., & Rädler, K.H. 1980, Meanfield magnetohydrodynamics and dynamo theory (Oxford: Pergamon Press) (In the text)
 Lazarian, A., & Vishniac, E. T. 1999, ApJ, 517, 700 [NASA ADS] [CrossRef] (In the text)
 Low, B. C. 1996, Sol. Phys., 167, 217 [NASA ADS] [CrossRef] (In the text)
 Livermore, P. W., Hughes, D. W., & Tobias, S. M. 2010, Phys. Fluids, 22, 037101 [NASA ADS] [CrossRef] (In the text)
 Matthaeus, W. H., Goldstein, M. L., & Smith, C. 1982, Phys. Rev. Lett., 48, 1256 [NASA ADS] [CrossRef] (In the text)
 Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D. 2009, ApJ, 697, 923 [NASA ADS] [CrossRef] (In the text)
 Mitra, D., Tavakol, R., Käpylä, P. J., & Brandenburg, A. 2010a, ApJ, 719, L1 [NASA ADS] [CrossRef] (In the text)
 Mitra, D., Candelaresi, S., Chatterjee, P., Tavakol, R., & Brandenburg, A. 2010b, Astron. Nachr., 331, 130 [NASA ADS] [CrossRef] (In the text)
 Parker, E. N. 1958, ApJ, 128, 664 [NASA ADS] [CrossRef] (In the text)
 Rädler, K.H., & Bräuer, H.J. 1987, Astron. Nachr., 308, 101 [NASA ADS] [CrossRef] (In the text)
 Rust, D. M. 1994, Geophys. Res. Lett., 21, 241 [NASA ADS] [CrossRef] (In the text)
 Seehafer, N. 1990, Sol. Phys., 125, 219 [NASA ADS] [CrossRef] (In the text)
 Stix, M. 1976, A&A, 47, 243 [NASA ADS] (In the text)
 Subramanian, K., & Brandenburg, A. 2006, ApJ, 648, L71 [NASA ADS] [CrossRef] (In the text)
 Sur, S., Subramanian, K., & Brandenburg, A. 2007, MNRAS, 376, 1238 [NASA ADS] [CrossRef] (In the text)
 Török, T., & Kliem, B. 2003, A&A, 406, 1043 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Warnecke, J., & Brandenburg, A. 2010, A&A, 523, A19 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Yoshimura, H. 1976, Sol. Phys., 50, 3 [NASA ADS] [CrossRef] (In the text)
Online material
Movie of Fig. 11 (Access here)
All Tables
All Figures
Fig. 1 Initial exponential growth and subsequent saturation behavior of the magnetic field in the interior for forced turbulence with dynamo action for Run A. The magnetic field strength is oscillating with twice the dynamo frequency 2ω_{cyc}. 

Open with DEXTER  
In the text 
Fig. 2 Equatorward migration, as seen in visualizations of B_{r} for Run D at r = R over a horizontal extent Δθ = 58° and Δφ = 17°. 

Open with DEXTER  
In the text 
Fig. 3 Periodic variation of and in the turbulence zone. Dark blue stands for negative and light yellow for positive values. The dotted horizontal lines show the location of the equator at θ = π/2. The magnetic field is normalized by the equipartition value. Taken from Run A. 

Open with DEXTER  
In the text 
Fig. 4 Time evolution of the radial magnetic field B_{r} (solid line) and the azimuthal magnetic field B_{φ} (dotted line) in the dynamo region averaged over the radius r and azimuth at θ = ± 7°. To improve the statistics, we calculate the components of the magnetic field as the antisymmetric part in latitude, i.e., for i = r,φ. 

Open with DEXTER  
In the text 
Fig. 5 Phase relations between the magnetic field and the velocity in the dynamo region. The magnetic field is plotted as , normalized with the equipartition field of the sound speed, (=) as a solid and black line. The rms velocity, normalized by the sound speed c_{s}, is plotted as a dashed red line, and has been smoothed over 5 neighboring data points to make it more legible. Taken from Run A. 

Open with DEXTER  
In the text 
Fig. 6 Radial dependence of density overplotted at different times. In the inset is the linear behavior of the logarithmic density log ρ/ρ_{0} to the inverse of the radius shown. Taken from Run A. 

Open with DEXTER  
In the text 
Fig. 7 Radial dependence of the mean squared magnetic field, (solid line), compared with those of B_{r} (dotted), B_{θ} (dashed), and B_{φ} (dashdotted). All quantities are averaged over 13 dynamo cycles. The inset shows the same quantities in a logarithmic representation. Taken from Run A. 

Open with DEXTER  
In the text 
Fig. 8 Time series of formation of plasmoid ejections in spherical coordinates. Contours of are shown together with a colorscale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the rθ plane. The dashed horizontal lines show the location of the surface at r = R. Taken from Run D. 

Open with DEXTER  
In the text 
Fig. 9 Time series of a reconnection event in an Xpoint as a closeup view. Contours of are shown together with a colorscale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the rθ plane, where solid lines represent counterclockwise magnetic field lines and the dash ones clockwise. The dashed vertical lines show the location of the surface at r = R. Taken from Run D. 

Open with DEXTER  
In the text 
Fig. 10 Magnetic field configuration at the time of a ejection. Contours of are shown together with a colorscale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the rθ plane, where the solid lines represent counterclockwise magnetic field lines and the dash ones clockwise. The dashed vertical lines show the location of the surface at r = R. Taken from Run D. 

Open with DEXTER  
In the text 
Fig. 11 Time series of coronal ejections in spherical coordinates. The normalized current helicity, , is shown in a colorscale representation for different times; dark blue stands for negative and light yellow for positive values. The dashed horizontal lines show the location of the surface at r = R. Taken from Run D. The temporal evolution is shown in a movie available as online material. 

Open with DEXTER  
In the text 
Fig. 12 Dependence of the dimensionless ratio on time t/τ and radius r in terms of the solar radius. The top panel shows a narrow band in θ in the northern hemisphere and the bottom one a thin band in the southern hemisphere. In both plots we have also averaged in latitude from 20° to 28°. Dark blue stands for negative and light yellow for positive values. The dotted horizontal lines show the location of the surface at r = R. 

Open with DEXTER  
In the text 
Fig. 13 Dependence of the dimensionless ratio on time t/τ at radius r = 1.5 R and 28° latitude. The solid line stands for the northern hemisphere and the dotted for the southern hemisphere. The red lines represent the cumulative mean for each hemisphere. 

Open with DEXTER  
In the text 
Fig. 14 Current helicity averaged over 3900 turnover times. Legend is the same as in Fig. 11. Dark blue corresponds to negative values, while the light yellow corresponds to positive value. Taken from Run D. 

Open with DEXTER  
In the text 
Fig. 15 Time evolution of the magnetic helicity flux of the largescale field, smooth over two data points. Here, the mean of magnetic helicity flux out through the surface of the northern hemisphere (black) is shown, together with that through the southern hemisphere (dotted red), and the equator (dashed blue). 

Open with DEXTER  
In the text 
Fig. 16 Cumulative mean of the time evolution of the magnetic helicity flux of the smallscale field, , normalized by , where η_{t} ≈ η_{t0} ≡ u_{rms}/3k_{f} was defined in Sect. 3.1. Here, the mean of magnetic helicity flux out through the surface of the northern hemisphere (black) is shown, together with that through the southern hemisphere (dotted red), and the equator (dashed blue). 

Open with DEXTER  
In the text 
Fig. 17 Dependence of the latitudinal component of the magnetic helicity flux, , compared with the latitudinal gradient of the magnetic helicity density of the smallscale field, , at r/R = 0.85. The largescale variation of the latter agrees with that of the former if the gradient is multiplied by an effective diffusion coefficient for magnetic helicity of κ_{t} ≈ 3η_{t0}. 

Open with DEXTER  
In the text 
Fig. 18 Helicity in the northern outer atmosphere. The values are written out at the point, r = 1.5 R, 90° − θ = 17°, and φ = 9°. Top panel: phase relation between the toroidal B_{φ} and poloidal B_{θ} field, plotted over time t/τ. Bottom panel: helicity H(k) is plotted over normalized wave number kR. The helicity is calculated with the Taylor hypothesis using the Fourier transformation of the poloidal and toroidal field. 

Open with DEXTER  
In the text 
Fig. 19 Helicity in the southern outer atmosphere. The values are written out at the point, r = 1.5 R, 90° − θ = −17° and φ = 8.6°. Top panel: phase relation between the toroidal B_{φ} and poloidal B_{θ} field, plotted over time t/τ. Bottom panel: helicity H(k) is plotted over normalized wave number kR. The helicity is calculated with the Taylor hypothesis using the Fourier transformation of the poloidal and toroidal field. 

Open with DEXTER  
In the text 