EDP Sciences
Free Access
Issue
A&A
Volume 534, October 2011
Article Number A11
Number of page(s) 10
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201117023
Published online 21 September 2011

© 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 missions1. 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 11-year 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 urms 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/vA, 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 dynamo-generated field in the convection zone with a nearly force-free 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 cs. 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 large-scale 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 dynamo-generated 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 self-consistent 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 stress-free 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 bi-layer 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   Mm2, 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 rate-of-strain tensor, semicolons denote covariant differentiation, is the specific pseudo-enthalpy, cs = const. is the isothermal sound speed, and g =  −GMr/r3 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 kfR ≈ 63. This value should be compared with the normalized wavenumber k1R, corresponding to the thickness of the shell ΔR, which yields k1R = 2πRR ≈ 21, so the effective scale separation ratio, kf/k1, 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 Af 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 = cs(|k|cs/δ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 kf. 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 non-helical forcing function, and is an arbitrary unit vector not aligned with k; note that .

The pseudo-enthalpy 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 small-scale fields is expected to be independent of the gauge (Subramanian & Brandenburg 2006; Hubbard & Brandenburg 2010), while that of the large-scale fields is not.

Our wedge is periodic in the azimuthal direction. For the velocity, we use stress-free 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 urms 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, Beq, 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 kf 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 Code2, which is a modular high-order code (sixth order in space and third-order 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 large-scale 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.

thumbnail 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

thumbnail Fig. 2

Equatorward migration, as seen in visualizations of Br for Run D at r = R over a horizontal extent Δθ = 58° and Δφ = 17°.

Open with DEXTER

thumbnail 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

Table 1

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 large-scale 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 non-dimensional 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 km is the relevant wavenumber of the mean field. Using ηt ≈ ηt0 ≡ urms/3kf (Sur et al. 2008), we find τωcyc ≈ 0.2(km/kf)2 ≈ 0.02, where we have assumed km ≈ 2π/0.3R ≈ 20k1 and kf ≈ 60k1. 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 non-monotonous dependence on the magnetic Reynolds number Rm. For larger values of Rm, the magnetic field strength decreases slightly with increasing Rm, 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 Rm is in our simulations approximately 15 times the critical Rm. 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 Rm dependency at Rm = 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 line-of-sight 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 out-of-phase (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 near-surface 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.

thumbnail Fig. 4

Time evolution of the radial magnetic field Br (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.

thumbnail 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 cs, 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.

thumbnail 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.

thumbnail Fig. 7

Radial dependence of the mean squared magnetic field, (solid line), compared with those of Br (dotted), Bθ (dashed), and Bφ (dash-dotted). 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

thumbnail Fig. 8

Time series of formation of plasmoid ejections in spherical coordinates. Contours of are shown together with a color-scale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the 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 close-up 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 counter-clockwise 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 right-hand side, a plasmoid has formed, which is eventually ejected. This plasmoid appears as a CME-like 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 gauge-independent. Current helicity has also been observed in the Sun (Seehafer 1990) and it has been obtained from mean-field 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 r-dependent time-averaged 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 Vej ≈ 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 Af and increases weakly with magnetic Reynolds number, but it seems to be still comparable to half the dynamo cycle period, i.e., Δτ ≈ Tcyc/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.

thumbnail Fig. 9

Time series of a reconnection event in an X-point as a close-up view. Contours of are shown together with a color-scale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the plane, where solid lines represent counter-clockwise 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 three-part 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 ball-shaped 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 “three-part 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 three-part 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 loop-like structures. An animation of the detailed time evolution of the CME-like structures emerging recurrently into the solar corona is available in the on-line 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 large-scale field.

thumbnail Fig. 10

Magnetic field configuration at the time of a ejection. Contours of are shown together with a color-scale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the plane, where the solid lines represent counter-clockwise 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.

thumbnail Fig. 11

Time series of coronal ejections in spherical coordinates. The normalized current helicity, , is shown in a color-scale 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

thumbnail 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

thumbnail 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 time-averaged 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.

thumbnail 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 small-scale magnetic field. Under the assumption of scale separation, this quantity is gauge-independent (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 small-scale 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 small-scale field was computed in Hubbard & Brandenburg (2010), where the magnetic helicity flux from the total and large-scale fields was found to be gauge-dependent, but that from the small-scale field was not. In Fig. 16 we compare the flux of magnetic helicity of the small-scale 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 Beq. 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 × 1047   Mx2 if we use ηt = 3 × 1012   cm2   s-1. This value exceeds the estimated upper limit for the solar dynamo of about 1046   Mx2 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 large-scale field may not a priori be gauge-invariant. 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 gauge-invariant. Furthermore, the shell-integrated 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 small-scale 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 large-scale 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.

thumbnail Fig. 15

Time evolution of the magnetic helicity flux of the large-scale 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

thumbnail Fig. 16

Cumulative mean of the time evolution of the magnetic helicity flux of the small-scale field, , normalized by , where ηt ≈ ηt0 ≡ urms/3kf 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

thumbnail Fig. 17

Dependence of the latitudinal component of the magnetic helicity flux, , compared with the latitudinal gradient of the magnetic helicity density of the small-scale field, , at r/R = 0.85. The large-scale 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 small-scale 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 CME-like 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 Vej (see Table 1) for translating temporal variations (in t) into spatial ones (in r) via r = r0 − Vejt. 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μ0EM/k, where EM 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).

thumbnail 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

thumbnail 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 CME-like phenomena are ubiquitous in simulations that include both a helicity-driven dynamo and a nearly force-free 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 self-consistent dynamo with a corona-like 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.


3

The movie is also available at http://www.youtube.com/watch?v=aR-PgxQyP24.

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. 621-2007-4064, and the National Science Foundation under Grant No. NSF PHY05-51164.

References

Online material

Movie of Fig. 11 (Access here)

All Tables

Table 1

Summary of runs discussed in this paper.

All Figures

thumbnail 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
thumbnail Fig. 2

Equatorward migration, as seen in visualizations of Br for Run D at r = R over a horizontal extent Δθ = 58° and Δφ = 17°.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 4

Time evolution of the radial magnetic field Br (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
thumbnail 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 cs, 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
thumbnail 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
thumbnail Fig. 7

Radial dependence of the mean squared magnetic field, (solid line), compared with those of Br (dotted), Bθ (dashed), and Bφ (dash-dotted). 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
thumbnail Fig. 8

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

Open with DEXTER
In the text
thumbnail Fig. 9

Time series of a reconnection event in an X-point as a close-up view. Contours of are shown together with a color-scale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the plane, where solid lines represent counter-clockwise 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
thumbnail Fig. 10

Magnetic field configuration at the time of a ejection. Contours of are shown together with a color-scale representation of ; dark blue stands for negative and light yellow for positive values. The contours of correspond to field lines of in the plane, where the solid lines represent counter-clockwise 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
thumbnail Fig. 11

Time series of coronal ejections in spherical coordinates. The normalized current helicity, , is shown in a color-scale 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. 15

Time evolution of the magnetic helicity flux of the large-scale 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
thumbnail Fig. 16

Cumulative mean of the time evolution of the magnetic helicity flux of the small-scale field, , normalized by , where ηt ≈ ηt0 ≡ urms/3kf 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
thumbnail Fig. 17

Dependence of the latitudinal component of the magnetic helicity flux, , compared with the latitudinal gradient of the magnetic helicity density of the small-scale field, , at r/R = 0.85. The large-scale 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
thumbnail 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
thumbnail 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.