Solar wind variations over a cycle with asymmetry

Simulations of solar wind variations during an 11-year cycle and the influence of north-south asymmetry

B. Perri\aff1    A. S. Brun\aff1    V. Réville\aff2,1       A. Strugarek\aff1 \aff1 AIM, CEA, CNRS, Université Paris-Saclay, Université Paris-Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France \aff2EPSS, University of California, Los Angeles, CA, United States

We want to study the connections between the magnetic field generated inside the Sun and the solar wind impacting Earth, especially the influence of north-south asymmetry on the magnetic and velocity fields. We study a solar-like 11-year cycle in a quasi-static way : an asymmetric dynamo field is generated through a 2.5D flux-transport model with Babcock-Leighton mechanism, and then is used as bottom boundary condition for compressible 2.5D simulations of the solar wind. We recover solar values for the mass loss rate, the spin-down timescale and the Alfvén radius, and are able to reproduce the observed delay in latitudinal variations of the wind and the general wind structure observed for the Sun. We show that the phase-lag between the energy of the dipole component and the total surface magnetic energy has a strong influence on the amplitude of the variations of global quantities. We show in particular that the magnetic torque variations can be linked to topological variations during a magnetic cycle, while variations in the mass loss rate appeared to be driven by variations of the magnetic energy.

1 Introduction

From a space weather perspective, one of the main challenges is to model accurately the solar wind, for it has a profound effect on the Earth’s space environment. Various in-situ observation missions, like Ulysses for example have shown that the magnetic activity has a strong influence on the wind structure and velocity, depending on the phase of the 11-year cycle (see McComas et al. (2008) and Smith (2011) for summary of the mission highlights, and Issautier et al. (2008) for discussion about north-south asymmetry). At minimum of activity, the magnetic field is mostly dipolar and the wind is slower at the equator (around 400 km/s) and faster at the poles (around 800 km/s); at maximum of activity, the magnetic field is multipolar and the wind distribution is bimodal at all latitudes. This is why there is more and more effort from an instrumental and theoretical point of view, to link in-situ space measurements and remote sensing solar surface observations, one of the goals of ESA’s Solar Orbiter mission.

Fortunately, we have a lot of observational data when it comes to the Sun magnetic activity : sunspots and magnetic field observations from several observatories and several time periods (Royal Greenwich Observatory from 1874 to 1954 in Newton & Milsom (1955) and from 1874 to 1976 in Vizoso & Ballester (1990), Wilcox Solar Observatory from 1976 to 2009 in Hoeksema (2010)) have shown a north-south asymmetry in the sunspots distribution. The southern hemisphere led by 18 months over cycle 19, and since then the northern hemisphere leads by a year on average (Svalgaard & Kamide, 2013). The maximum delay between the two hemispheres measured so far is 2 years and it is worth noting that no systematic pattern has been found for this polar reversal phase delay in the surface magnetic activity (Temmer et al., 2006). We can note however that a systematic asymmetry can be observed for the average heliospheric current sheet at the Earth’s orbit which appears to be systematically shifted southwards since at least cycle 16 (Mursula & Hiltula (2003)).

A possible explanation of this asymmetry has been proposed : it can be due to the interference between the dipolar and quadrupolar components of the magnetic field can lead in extreme cases to the vanishing of the magnetic field in one hemisphere (Gallet & Pétrélis, 2009). See also Tobias (1997), Ossendrijver (2003), DeRosa et al. (2012), Shukuya & Kusano (2017) and references therein for a detailed discussion on dynamo symmetry properties. The antisymmetric family corresponds to odd degrees when projecting the magnetic field onto the spherical harmonics functions under the assumption of axisymmetry (with order ). It is overall dominant in the Sun, which explains the apparent dipolar structure over most of the cycle. However the symmetric family, which corresponds to even degrees (when ), is not negligible; it reaches on average 25% of the amplitude of the antisymmetric family, and becomes dominant during polarity reversals. Such asymmetry has also an impact on the wind structure : in Sokół et al. (2015), the reconstruction of the solar wind using Interplanetary Scintillations (IPS) observations summed up in Tokumaru et al. (2010), shows clearly an asymmetry in the latitudinal distribution of the wind.

For the large-scale magnetic field generation, the general theoretical framework is the dynamo theory, especially the interface dynamo (Parker, 1993) : due to strong shears in the solar differential rotation (Schou et al., 1998), a toroidal field is generated in the tachocline; the poloidal field is regenerated by induction thanks to the turbulent fluid motions in the convection zone, thus sustaining it against ohmic dissipation (see Miesch (2005) for review on the subject). The most realistic models are associated with a flux-transport mechanism and a Babcock-Leighton term to take into account the role of the meridional circulation, linking the surface of the Sun and the tachocline at the poles by redistributing the magnetic field (see Babcock (1961) and Leighton (1969), but also Wang & Sheeley (1991) for the link with observations and Dikpati & Charbonneau (1999) for simulations).

To simulate the whole Sun, magnetohydrodynamic (or MHD) simulations have been used successfully since the 1970s (for a full review on the subject, see Brun & Browning (2017)).

Large Eddy Simulations (LES) are global-scale simulations where a filter is applied to remove the small scales which are treated as dissipation terms. They have been first computed in Gilman (1983) and in Glatzmaier (1985), and first adapted for high resolutions in Brun et al. (2004); more realistic set-ups with high Reynolds numbers can be found in Hotta et al. (2016). LES models are now able to display large-scale magnetic cycles : for the Sun see Ghizaru et al. (2010), for young convective stars see Brown et al. (2011), for exceptional cycles such as grand minima see Augustson et al. (2015) and for generalization to other solar-like stars see Strugarek et al. (2017).

Mean Field Simulations (MFS) are simulations performed under the assumptions of mean-field theory and axisymmetry. They have the advantage of being computationally cheap compared to the other fully 3D MHD methods. A presentation of kinematic dynamo models can be found in Roberts (1972), for more focus on the role of rotation see Stix (1976); for general reviews see Krause & Raedler (1980) and Charbonneau (2010). Though they rely on a simplified physical model, MFS reproduce better the magnetic field at the surface of the Sun for now. New models are even able to produce sunspots (Kumar et al., 2018).

For the solar wind, the current framework has been initiated by the work of Parker (1958) : the observed gap of pressure and temperature between the solar corona and the interstellar medium led to the emerging idea of a dynamic corona, expanding to supersonic speed. This hydrodynamical approach was then expanded to take into account the magnetic field, which yields better computations of the angular momentum loss (Schatzman (1962), Weber & Davis (1967)). It was also realized that, to explain the fast distribution of the wind, we certainly needed to take into account subtle magnetic effects and to adopt a multi-fluid approach (Hollweg & Isenberg, 2002). These 1D analytical equatorial models were then expanded to 2D (Sakurai, 1985), leveraging the conservation of various physical quantities along the poloidal streamlines. There are still a lot of mechanisms that are not fully understood and yet to be modeled. For instance, the coronal heating is still an open problem : we do not know the exact mechanism, although there are very promising studies around magnetic reconnection (Parker, 1988) and MHD turbulence (Cranmer et al., 2007). The open flux problem is also a rising challenge, with the question of observational magnetic maps underestimating the open magnetic flux of the Sun (Linker et al., 2017).

One way to reproduce the wind dynamics is via compressible MHD simulations (Keppens & Goedbloed, 1999) and there are basically two approaches : study of one or several flux tubes starting from the chromosphere to focus on the energy transfer between the surface and the corona in one (Suzuki & Inutsuka, 2005) or either in two dimensions (Matsumoto & Suzuki, 2012), or more global simulations including the whole star either in 2.5D in Matt et al. (2012), Tóth et al. (2012), Réville et al. (2015a) or in 3D in Riley et al. (2015) and Réville & Brun (2017).

Though connected in the real Sun, the solar interior and atmosphere are very difficult to couple numerically. Different plasma parameters, the wide range of time and length scales, the variety of physical processes involved and the stiffness of the MHD equations all together make the modeling of the complete Sun an obvious numerical challenge. Remarkable attempts to deal with this problem are numerical studies made in small Cartesian domains (few tens of Megameters) including the different photospheric layers (Vögler et al. (2005), Martínez-Sykora et al. (2008) for the addition of conduction, Rempel et al. (2009) for focus on sunspots among others ; see also review by Wedemeyer-Böhm et al. (2009)).

In a previous study, published in Pinto et al. (2011), the influence of the magnetic field geometry and amplitude on the solar wind has been studied in a quasi-static way in 2.5D. A first numerical code named STELEM (STellar ELEMents, cf. Jouve & Brun (2007)) computes an 11-year magnetic cycle using a kinematic mean-field dynamo model. The latter is used as a bottom boundary condition for a second numerical code (DIP, cf. Grappin et al. (2010)) which computes a sequence of relaxed steady states of an isothermal wind. Quantities such as the wind speed distribution, the mass loss and the angular momentum loss were thus computed over an activity cycle, showing considerable variations over time.

In this paper, we use a similar approach but with a more realistic dynamo model as described in Jouve & Brun (2007) and emphasizing north-south asymmetry as in DeRosa et al. (2012) (see their Appendix A). We will first present the numerical codes and the physical ingredients in our models in Section 2, then we will present our results as followed : we will begin with a description of the magnetic dynamo cycle in Section 3 (time-latitude diagram, topology of the field and its evolution), then we will focus on the variation of the solar wind speed and its spatial distribution in Section 4 and finally we will describe the evolution of global quantities such as mass and angular momentum loss over the cycle under the influence of asymmetry in Section 5. Discussion and conclusion are made in section 6.

2 Numerical set-up

We used two different MHD codes for our 2.5D axisymmetric simulations : STELEM for the dynamo field generated inside the star with a flux-transport model (Emonet & Charbonneau 1998, private communication; Jouve & Brun (2007)), and PLUTO for the solar corona and the associated wind (Mignone et al., 2007). The coupling between the two codes is made through the magnetic field properties : the dynamo field is matching a potential field at the surface at the star, which is then used as a bottom boundary condition for the corona and wind, and initial background field over the whole computational domain.

In Section 2.1 we describe physical properties and the numerical methods used in the STELEM code; the same is done for the PLUTO code in Section 2.2. Finally we give further details about the coupling process in Section 2.3.

2.1 Dynamo simulations with STELEM

The simulations performed are the same as the one described in DeRosa et al. (2012) in their appendix.

Working in spherical coordinates and under the assumption of axisymmetry, we perform a poloidal-toroidal decomposition and write the mean magnetic and velocity field (respectively and ) as follows :


is decomposed with the poloidal stream function and toroidal field . The velocity field is time-independent and prescribed by profiles of the meridional circulation and differential rotation .

We rewrite the induction equation in the framework of mean-field theory in terms of and and we introduce a poloidal term based on the Babcock-Leighton mechanism. We finally normalize the equations by using the solar radius as length scale and the diffusion time as the time scale, and obtain the following two coupled partial differential equations :


where is the effective magnetic diffusivity, is the turbulent diffusivity in the convection zone and . These equations are controlled by three dimensionless parameters : which characterizes the shear by differential rotation (i.e. the omega effect) ; which characterizes the Babcock-Leighton source term ; (the Reynolds number) which characterizes the intensity of the meridional circulation. , and are respectively the rotation rate, the typical amplitude of the surface source term and the amplitude of the meridional flow. In this study, we have , and , which corresponds to the parameters in DeRosa et al. (2012). The rotation rate thus correspond to the one measured by helioseismology.

For simplicity, we assume that the meridional circulation is equatorially symmetric, having one large single cell in each hemisphere. Flows are directed poleward at the surface and equatorward at depths, vanishing at the bottom radial boundary. The equatorward flow penetrates slightly beneath the tachocline (Dikpati & Charbonneau, 1999).

The rotation profile is inspired by solar angular velocity profile deduced from helioseismic inversions (Thompson et al., 2003). We assume solid-body rotation below and a differential rotation above :


with the parameters : , , , , and .

We assume different diffusivities in the envelope and in the stable interior, smoothly matching the two :


with the parameters : , and .

In Babcock-Leighton flux-transport dynamo models, the poloidal field owes its origin to the tilt of magnetic loops emerging at the solar surface. Thus the source has to be confined to a thin layer just below the surface and, since the process is fundamentally non-local, the source term depends on the variation of at the base of the convection zone. To create asymmetry between the two hemispheres we introduce a modified source term modulated by the asymmetry parameter :


with the parameters : , , and .

The STELEM code is used to solve equations (3) and (4) : it uses a finite element method in space and a third-order scheme in time. The temporal scheme used is adapted from Spalart et al. (1991) and is similar to a Runge-Kutta 3 method. The STELEM code has been thoroughly tested and validated via an international mean-field dynamo benchmark involving eight different codes (Jouve et al., 2008).

The numerical domain is an annular meridional cut of the Sun with the colatitude and the radius (ie. from slightly below the tachocline located at up to the solar surface). We use a uniform grid with 64 points in radius and cosine of the latitude. At the latitudinal boundaries ( and ) and at the bottom radial boundary (), and are set to 0. At the upper radial boundary (), the solution is matched to an external potential field. Usual initial conditions involve setting a confined dipole field configuration, i.e. is set to in the convection zone and to 0 below the tachocline ; the toroidal field is initialized to 0 everywhere.

2.2 Wind simulations with PLUTO

This set-up is inspired by the one described in Réville et al. (2015a), but adapted to spherical coordinates.

We solve the set of the conservative ideal MHD equations composed of the continuity equation for the density , the momentum equation for the velocity field with its momentum written , the energy equation which is noted and the induction equation for the magnetic field :


where is the total pressure (thermal and magnetic), is the identity matrix and is a source term (gravitational acceleration in our case). We use the ideal equation of state :


where is the thermal pressure, is the internal energy per mass and is the adiabatic exponent. This gives for the energy : .

PLUTO solves normalized equations, using three variables to set all the others: length, density and speed. If we note with the parameters related to the star and with the parameters related to the normalization, we have , and , where is the Keplerian speed at the stellar surface and the gravitational constant. By choosing the physical values of , and , one can deduce all of the other values given by the code in physical units. In our set-up, we choose cm, (which corresponds to the density in the solar corona above 2.5 km, cf. Vernazza et al. (1981)) and km/s. Our wind simulations are then controlled by three parameters : the adiabatic exponent for the polytropic wind, the rotation of the star normalized by the escape velocity and the speed of sound normalized also by the escape velocity . Note that the escape velocity is defined as . For the rotation speed, we take the solar value, which gives . We choose to fix , which corresponds to a K hot corona for solar parameters and . This choice of is dictated by the need to maintain an almost constant temperature as the wind expands, which is what is observed in the solar wind. Hence, choosing is a simplified way of taking into account heating, which is not modeled here.

We assume axisymmetry and use the spherical coordinates . Since PLUTO is a multi-physics and multi-solver code, we choose a finite-volume method using an approximate Riemann Solver (here the HLL solver, cf. Einfeldt (1988)). PLUTO uses a reconstruct-solve-average approach using a set of primitive variables to solve the Riemann problem corresponding to the previous set of equations.

The numerical domain is an annular meridional cut with the colatitude and the radius . We use an uniform grid in latitude with 512 points, and a stretched grid in radius with 512 points; the grid spacing is geometrically increasing from at the surface of the star to at the outer boundary. At the latitudinal boundaries ( and ), we set axisymmetric boundary conditions. At the top radial boundary (), we set an outflow boundary condition which corresponds to for all variables, except for the radial magnetic field where we enforce . Because the wind has opened the field lines and under the assumption of axisymmetry, this ensures the divergence-free property of the field. At the bottom radial boundary (), we set a condition similar to the one described in Zanni & Ferreira (2009) to be as close as possible to a perfect rotating conductor. We also implement the same differential rotation as described in (5). We initialize the velocity field with a polytropic wind solution and the magnetic field with a potential extrapolation of the field produced by STELEM at a given time. Please note that there is a splitting in PLUTO between the background field (which is curl-free and provided by the dynamo model) and the deviation field (which is a perturbation of the background field and carries the magnetic energy). Please note also that, to enforce the divergence-free property of the field, we use a hyperbolic divergence cleaning, which means that the induction equation is coupled to a generalized Lagrange multiplier in order to compensate the deviations from a divergence-free field (Dedner et al., 2002).

2.3 Coupling method

To couple the two codes described above, we use the projection of the surface magnetic field produced by STELEM on spherical harmonics, which PLUTO can read. To proceed so, we select a given time , then we project the surface radial magnetic field at that time on spherical harmonics up to a degree :


Note that there is no dependency in due to axisymmetry (). In this study, we choose , as it provides the best compromise between accuracy and costs for the reconstruction. Then inside PLUTO, we reconstruct the field by reading these coefficients and combining them with spherical harmonics, and finally we extrapolate the field inside the whole wind computation domain. For the extrapolation, we chose a Potential-Field Source Surface (PFSS) method with a source surface radius equal to . For more information about the PFSS method, for original computation see Schatten et al. (1969) and Altschuler & Newkirk (1969), for more explicit calculation see Schrijver & De Rosa (2003). We obtain:


In Réville et al. (2015b), it was explained that the fiducial value of 2.5 usually found in the literature for was an approximation, and that you can define an optimal source surface radius that allows you to recover the best value for the open flux in your simulation. This optimal value increases with the magnetic field strength and decreases with high order magnetic topology and rotation rate. It was suggested that a different value of should thus be used at different times over the cycle, but we did not investigate such parametrization as the PFSS is only used to initialize our wind model. Our value is an estimate of the larger value needed in our simulations, which correspond to a strong dipole at low rotation rate. We recall also that the PFSS extrapolation is only used to initialize the simulation, and the final configuration does not depend much on the initial extrapolation.

This yields the following field reconstruction :


The magnetic field provided by STELEM is dimensionless (cf. (3) and (4)). To recover the proper physical units, we calibrate the radial magnetic field amplitude so that the mass loss rate is as close as possible to estimations of the global mass loss rate deduced from Ulysses data (McComas et al., 2008) ; this yields values between and (Réville & Brun, 2017). The magnetic field is then normalized by the PLUTO units described in section 2.2.

From our dynamo run, we obtain thus a times series of 54 snapshots , sampled over an 11-year cycle with a time step of about 2 months. For each input magnetic field, we let the wind relax and reach a steady state ; this takes around 500 characteristic times (defined as ), which translates to 9 days using solar units. The result is then a sequence of steady-state solutions, providing general properties of the coronal magnetic field and the wind between solar minimum and maximum of activity, but without taking into account the back-reaction of the wind on the dynamo. All simulations are performed from scratch, so there is no dependence for any state of relaxation on the previous states.

3 Description of the cycle

3.1 Properties of the dynamo field

Figure 1: Time-latitude diagrams of the radial surface and the azimuthal tachocline magnetic fields for the dynamo model produced by STELEM. The time is shown in years and the vertical axis shows the cosine of the angle associated to the latitude. The left panel shows the general aspect of the magnetic field over 3 cycles. The right panel is a zoom on the specific cycle we studied, with the black dashed lines indicating remarkable times of the simulation.

First we will begin by presenting the main features of the solar dynamo solution generated by STELEM.

Figure 2: Evolution of the proxy for the sunspot number (SSN) over time. The SSN for the northern (southern) hemisphere is in blue (red), and the total number is in black. The black dashed lines indicate the minima of the studied cycle as defined by topology.

The left panel of Figure 1 displays the time-latitude diagrams for the radial magnetic field at the surface of the star and the toroidal magnetic field at the base of the convective zone over several dynamo cycles. The cycle period is approximately 12.0 years, which is close to the observational mean Sun value of approximately 11 years (Clette & Lefèvre, 2012). The asymmetry of the model results in a delay of 9 months of the magnetic configuration between the two hemispheres (for example the southern one reverses first, like in cycle 18 of the Sun). This is in qualitative agreement with the current observed delay of 1 year (Temmer et al., 2006). This diagram for the surface radial magnetic field is typical of a flux-transport dynamo model, and exhibits features similar to the solar field : at each instant of the cycle there are two branches, one equatorward and the other one poleward. The right panel of Figure 1 shows more precisely which cycle we decided to study, with black dashed lines indicating some remarkable times at which we computed the associated wind solution with PLUTO. These times are 0.0, 2.9, 4.2, 4.9, 5.8 and 12.0 years after the cycle minimum. We will explain in section 3.2 why we chose to show these specific moments among the 54 we computed. The most right and most left lines correspond to the cycle minima.

Note that there are several ways to define a cycle minimum in our numerical simulations when comparing with real sunspot time series. In Figure 1, we fixed the minimum as the time of the cycle when the magnetic field configuration is more dipolar, which means the time when the ratio of the dipole energy over the quadrupole energy is maximum. Likewise, the maximum of the cycle is defined as the maximum ratio of energy between the quadrupole and dipole. Another way to define a minimum of activity is to say that it is the time when there are the lowest number of sunspots on the solar surface. Since our model does not generate sunspots, we use a proxy to determine an equivalent sunspot number (or SSN) based on the generation of toroidal magnetic field at the tachocline, given by:


where is a scaling constant to find values for our SSN proxy of the same order of magnitude as the Sun (Hung et al., 2017).

Separating the SSN from the northern and southern hemisphere in Figure 2 allows us to see clearly the north-south asymmetry. For the Sun, these two definitions of minimum or maximum of the solar cycle give the same time. However, we can notice that in our model the minimum of sunspot activity is delayed by 2.5 years compared to our minimum of quadrupolar energy (cf. Figure 2). This is because we did not fine-tune the model to match an exact solar cycle, what we seek foremost is to understand the link between the dynamo field and the wind on a fundamental level. However it allows us to isolate the effects of those different contributions and understand more precisely the underlying physics.

Figure 3: Time evolution of the coefficients of the surface magnetic field projection on spherical harmonics over 2.5 dynamo cycles, grouped as equatorially symmetric and antisymmetric families. Only the modes are displayed due to the assumption of axisymmetry.

Another impact of the introduced asymmetry is the ability to couple the equatorially symmetric and antisymmetric family modes for the magnetic field. To demonstrate this point, we display in Figure 3 the time evolution of the coefficients of the projection of the surface radial magnetic field on the spherical harmonics. They are gathered as equatorially symmetric and antisymmetric families, which correspond to the even and odd degrees for the projection on the spherical harmonics when considering only modes. First thing we can notice is that the amplitude of the antisymmetric family modes is similar to the one of the Sun (between -4 and 4G, as shown in DeRosa et al. (2012)). The component is stronger than what is observed in the Sun : this seems to be induced by the flux-transport model, which tends to accumulate magnetic flux at the poles, thus favoring higher modes than the dipole and octupole. Despite the fact that the simulation was initialized with a dipole field, we can see that the symmetric family modes develop in a significant way, reaching on average about 10% of the antisymmetric family mean amplitude. This is different from simple 2.5D mean-field dynamo models : usually the symmetric family modes are unable to develop when the model is initialized with a dipole and is symmetric, whereas in the Sun the quadrupole amplitude is measured to be around 25% that of the dipole most of the time. We can also note that the amplitude of the symmetric family modes (between -0.5 and 0.5G) is close to what has been observed since cycle 21 (see for instance DeRosa et al. (2012)).

Figure 4: Comparison of the evolution of the surface magnetic energy and the energy of 2 mode components over several cycles : for the left panel, for the right panel.

The final property of our dynamo model we wanted to highlight is the time evolution of the total radial surface magnetic energy versus the energy of the dipole component. With our decomposition in spherical harmonics (cf. (16)), we can define the total radial surface magnetic energy as :


Then we can define the energy of a specific harmonic component of the radial surface field as :


Using data from the Wilcox Observatory it can be shown that and were anti-correlated from cycle 20 to 23 (Réville & Brun, 2017), which confirms that the Sun is mostly dipolar during minimum of activity and multipolar near maximum. We plot the time-evolution of these two quantities over 3 cycles in the left panel of Figure 4. In our case, and have a phase delay of one quarter of a period, so not fully correlated nor anti-correlated but in phase quadrature. This is a direct consequence of what we have just highlighted concerning the modes : since the dipole is not here the strongest magnetic mode, it is not the relevant one to study from an energetic point of view. In the right panel of Figure 4, we show the same comparison but for the mode, and here we have a phase delay of one third of a period, which is closer to anti-correlation. This shows that, although Babcock-Leighton models seem to allow higher modes to reach a bigger amplitude than in the Sun, they capture most of the interplay between topology and energy.

3.2 Coronal magnetic field

The time evolution of the coronal magnetic field is displayed in Figure 5 at remarkable times during the cycle. If we consider the initial minimum of activity at the beginning of our simulation as instant , the different panels a, b, c, d, e and f correspond respectively to times 0.0, 2.9, 3.8, 4.5, 6.0 and 11.9 years. These times were indicated as black dashed lines in the right panel of Figure 1. The first 4 close to the star are visible for both hemispheres in order to see the effect of asymmetry. The color scale represents the following quantity : , which is the solar wind velocity projected on the magnetic field in units of Mach number. Since the wind always flows outwards the star, this quantity allows us to track the changes of polarity of the magnetic field in the open field regions ; with this colortable, yellow corresponds to positive polarity, and dark blue to negative polarity. The transitions between polarities are in red and correspond to current sheets. The poloidal field lines are plotted in white, in solid (dashed) for positive (negative) polarity, corresponding to a positive (negative) potential vector .

Figure 5: Snapshots of the evolution of the coronal magnetic field at different times : if we set years at the first minimum, then we have 0, 2.9, 4.2, 4.9, 5.8 and 12.0 years. The color scale represents the following quantity : , which is the solar wind velocity projected on the magnetic field in units of Mach number. White lines correspond to the poloidal magnetic field lines of positive polarity in solid and negative polarity in dashed lines. We represent only the 4 first solar radii.

We now analyze the changes of topology in the coronal magnetic field observed over one dynamo cycle and its interaction with the wind. To describe properly the structures, we will differentiate the helmet streamers, which separate coronal holes of opposite magnetic polarities, from the pseudo-streamers, which overlie twin loop arcades and separate holes of the same polarity, as defined in Wang et al. (2007).

We start at a minimum of activity, go through a maximum and then return to minimum. At minimum of activity (panel a), the magnetic field is dominated by the dipole starting from 2-3 solar radii, but is mostly octupolar at the surface with 3 arcades of closed magnetic field loops. This is consistent with the spectrum analysis displayed in Figure 3. We have a central streamer similar to what is observed in the Sun : it extends up to 4 solar radii, where most coronograph pictures show a streamer between 2 and 4 solar radii. The magnetic field is of positive polarity in the northern coronal hole and of negative polarity in the southern coronal hole.

We can then clearly see the reversal of the two hemispheres happening one after the other. In panel b, we can see that the equatorial streamer has been disrupted at high latitudes in the southern hemisphere, leading to the appearance of pseudo-streamers. In the northern hemisphere, the main streamer is still intact and connected transequatorially with the southern hemisphere just under the equator. In panel c, the first coronal hole of opposite polarity opens in the southern hemisphere, creating a streamer at mid-latitude. The equatorial streamer is completely disrupted at this point, pseudo-streamers have been formed as well in the northern hemisphere. In panel d, another coronal hole of opposite polarity forms this time in the northern hemisphere. In the southern hemisphere, the coronal hole is spreading due to the wind opening up the coronal field lines. The magnetic configuration is very close to quadrupolar. In panel e, the southern hemisphere has now reversed, since the overall magnetic field, especially at the poles, is now the opposite of the one in panel a. In the northern hemisphere, the reversal is not complete, the pole is the last region where the coronal hole of opposite polarity has not spread to yet. In panel f, we jump from the middle of the cycle to the end of the cycle: the field has returned to a dipolar configuration, but with the exact opposite polarity compared to panel a, hence completing an activity cycle.

Figure 6: Radial cuts at minimum and maximum of activity, for , and between 2 and 10 solar radii.

We then describe further our model by showing radial cuts of various quantities in the solar corona. We want to check the radial dependency of , and to see if we recover expected behaviors. In Figure 6, we plot these quantities in the equatorial plane for the velocity and at mid-latitude for the magnetic field to avoid the current sheets, at minimum and maximum of activity. For the velocity profile, we compare it at minimum with predictions of the Weber and Davis wind, as shown in Keppens & Goedbloed (1999). The radial velocity has the profile of an accelerated transonic solution, which at the end is well fitted by a logarithmic function of the radius. The longitudinal velocity is slightly increasing in the first solar radii before slowly decreasing. Its initial slope is equal to the rotation rate of the star, because of the co-rotation of the corona with the star. Further from the star, this becomes less and less true and the surrounding corona rotates slower; we tend to have a dependency. Far from the star, displays an expected dependency of , typical of a purely radial field stretched by the wind to satisfy the divergence-free property (in spherical coordinates : ). This is true both at minimum and maximum of activity.

4 Solar wind speed

We will now focus on the variation of the wind speed over a solar cycle. This is one of the most important features because it can be compared with in situ measurements and directly affect planets or satellites. Figure 7 shows the evolution of the radial wind velocity at in km/s over an 11-year cycle in our simulation. We show only the radial velocity because so far from the star, the wind has become mostly radial, so it is which contains most of the wind total velocity. We chose because it is our domain’s outer boundary, which corresponds to almost 0.1 AU; most wind in situ measurements are made at 1 AU, so we have to extrapolate our numerical values to compare them with observational values by using a logarithmic fit.

Figure 7: Time-latitude diagram of the wind speed over an 11-year cycle in km/s at .

First of all, we can see that the amplitude of the wind is different from what is measured : the solar wind flow amplitude is usually between 400 and 800 km/s, while here we reach between 200 and 235 km/s at 0.1 AU, which gives between 420 and 470 km/s at 1 AU. This difference in amplitude is due to our polytropic corona that prevents us from recovering the fast wind, because we have a uniform heating. To recover fast winds numerically, we have to make sure a large part of the energy is deposited beyond the critical point in the supersonic region (Leer & Holzer, 1980). Such complex developments will be undertaken in a later study.

We will focus then on the latitudinal distribution of the wind. At minimum of activity, the wind is very organized : at the equator, the magnetic structures being mainly closed loops, the wind tends to be trapped in the corona and is slower ; at higher latitudes, the magnetic field is more open, which allows the wind to reach its maximum speed. As the star approaches maximum of activity, the topology becomes more complex, with more closed loops forming at the surface of the star, which reduces the wind velocity. What is very interesting in our model is that, due to asymmetry, we can clearly see the wind slowing down first in the southern hemisphere only around year 3 after the first minimum, and then in the northern hemisphere around year 4 ; then it speeds up again at high latitudes around year 6 for the southern hemisphere and year 7 for the northern hemisphere. At , the time delay of the wind velocity between the northern and southern poles is 5 months. Compared to 9 months for the delay in the magnetic field asymmetry, we can see that the wind tends to smooth the asymmetry. It is difficult to say if this tendency is recovered in observations: in Tokumaru et al. (2015), the asymmetry in the wind speed is found to be at most a year, with the possibility of being less; for the corresponding cycles, the asymmetry in the sunspot numbers is between one and two years (Svalgaard & Kamide, 2013). So, like in the model, there seems to be a reduction of the north-south asymmetry between the surface and far away in the wind.

When we compare this result with the same figure from Pinto et al. (2011), we can notice a few differences. Of course, the asymmetry was absent from their model, which means that their latitudinal slow-downs are perfectly synchronized. Also, in their model, the wind seems to return almost immediately to a minimum-like state after the maximum, while in ours, it really takes most of the second half of the cycle to reach an equatorial zone of slow-down which is as narrow as in the minimum state of activity. This feature in our model is actually closer to what is observed in the Sun : solar wind reconstruction using IPS by Tokumaru et al. (2010) and Sokół et al. (2015) display the same behavior, with some kind of relaxing time between the maximum and minimum of activity of 7.5 years, compared to 6 in our model. We can also point out that, as asymmetry is present in the Sun, we can also notice a delay between the latitudinal variations of the wind in Sokół et al. (2015).

Finally, we can notice that during maximum, the wind velocity is dropping at the poles in the two hemispheres. This has been observed in Ulysses data and has been interpreted as the presence of high-latitude current sheets (Khabarova et al., 2017).

Figure 8 displays the Alfvén surface profile at minimum and maximum (respectively left and right panel). The Alfvén surface corresponds to the surface at which the wind speed becomes faster that the Alfvén speed . At minimum, the Alfvén surface is very regular. The Alfvén radius is bigger at the poles than at the equator, due to the magnetic field being stronger. At maximum, we can see the Alfvén surface becoming more irregular due to the change of topology. However at the poles, the Alfvén radius remains pretty large, due to both the magnetic field still being stronger and the drop in the wind velocity observed in Figure 7 between year 2.5 and 5.5.

We see that our model is qualitatively closer to observations thanks to the asymmetric magnetic field produced by the underlying dynamo model.

Figure 8: Meridional cuts of the Alfvén surface (black line) at minimum and maximum of activity (respectively first and second panel). The color in the background shows the velocity projected on the magnetic field in unit of Mach number. The white lines are the poloidal magnetic field lines of positive polarity in solid and negative polarity in dashed lines. We represent only the 10 first solar radii.

5 Mass and momentum flux

In this section, we will focus on quantities of interest and their time evolution over the cycle in our simulations. The three quantities we are interested in are the average Alfvén radius, the mass loss rate and the angular momentum loss rate.

Figure 9: Time evolution of several global quantities of interest over an 11-year cycle : average Alfvén radius (first panel), mass loss (second panel) and angular momentum loss (third panel). The purple (orange) ticks indicate the beginning and end of the field reversal in the southern (northern) hemisphere. The grey line is either the parity factor or the surface magnetic energy as indicated in the right y-axis.

The Alfvén radius is defined as the radius at which the wind velocity is equal to the Alfvén velocity. We then define an average Alfvén radius (as in Pinto et al. (2011)) as the average cylindrical radii of the Alfvén surface weighted by the local mass flux crossing the surface of a sphere :


The mass loss rate is defined as :


where is the radius at which you want to measure the mass loss rate (ie. the radius of the spherical integration surface). In theory, the mass loss rate should be independent of . In our case, the mass loss is evaluated at the outer boundary of the domain.

In the same way, we define the angular momentum loss rate as :


where and are respectively the poloidal magnetic field and velocity flow.

Figure 9 shows the evolution of these three quantities over a solar cycle, starting at minimum of activity. Since our wind solution is being relaxed for a given time and a given magnetic field configuration, it means that we can make temporal averages. Note then that the curves we show in Figure 9 are averaged in time over 25 reference times, which means that the variations observed are significant compared to temporal variations. We overplotted two relevant quantities to understand our variations : the magnetic surface energy, already defined in equation (18), and the parity factor defined as follows:


where is the energy of the quadrupole and is the energy of the dipole.

The average Alfvén radius goes from 5.5 solar radii at minimum to 3.0 solar radii at maximum. These values are slightly smaller than what was estimated from the angular momentum loss from the Helios mission in Pizzo et al. (1983) and Marsch & Richter (1984), where its largest value was estimated at 12-14 and 13.6-16.6 solar radii respectively. However, when we compare our values to numerical models, the range varies between 2.5 and 60 : isothermal and polytropic models tend to give a lower estimate for the Alfvén radius (Pneuman & Kopp, 1971). When we compare with Pinto et al. (2011) to see the impact of the Babcock-Leighton model versus the alpha-omega dynamo model, we see that their Alfvén radius goes from 9 at minimum to 2.2 at maximum. The values are quite similar, it is more the amplitude of the variation over the cycle which has been modified : there is a factor 4 in their case, a factor 2 in our case. When we compare with Réville & Brun (2017) for the impact of 2D versus 3D and theoretical dynamo versus observational maps, we see that their Alfvén radius goes from 5 at minimum to 7 at maximum. The amplitude of the variation is slightly smaller. An explanation proposed for these discrepancies lies in the variations of the dipole component and the correlation with the total magnetic energy, as shown in Figure 4 : in our case, contrary to Pinto et al. (2011), we do have an anti-correlation, which is not perfect but does capture most of this effect, which could explain our intermediate results.

Figure 10: Mass flux in the meridional plane at different times in the cycle : the first panel is at year from the first minimum of activity, the second at years, the third at years and the last at years. Bright shades indicate higher mass loss rate. The white lines are magnetic field lines of positive (negative) polarity in solid (dashed) lines. We represent only the 4 first solar radii.

The mass loss we compute varies from 2.35 to 2.70 at maximum. From the data in McComas et al. (2008) and the calculations from Réville & Brun (2017), the mass loss inferred by Ulysses reached a minimum value of 2.3 in 1991 and a maximum value of 3.1 in 1992-1993. Our minimum mass loss is correct, but the amplitude of the variations is smaller. Wang (1998) indicates that the mass loss rate is minimum at minimum of activity, then rises during the maximum, and then returns to its original value as the star goes back to minimum. Our mass loss has a different behavior : it decreases till it reaches it minimum value around year 3, then rises and reaches its maximum around year 8, before decreasing again until the end of the cycle. This behavior can be explained by the fact that our maximum of polarity and maximum of magnetic energy do not happen at the same time : our mass loss is mostly anti-correlated with the evolution of the magnetic surface energy, which is also the case in the model of Réville & Brun (2017). The variation between the mass loss at minimum and maximum is smaller than in most simulations (13% in our case, factor 2 in Pinto et al. (2011), 20% in Réville & Brun (2017)). This small variation of the mass loss may well be a consequence of the delay between the minimum of activity and the minimum defined by topology. We recall also that the geometry of the magnetic field has not been fine-tuned to represent any solar cycle in particular.

Figure 10 shows a 2D map of the mass flux in the meridional plane ; this allows us to determine the latitude at which the mass loss is more important. We plotted different relevant moments of the evolution of the mass loss rate : at minimum (year 0), during the local minimal (year 2.2), during the rising phase (year 5.1) and at the maximal value (year 7.6). At minimum of activity, we have a repartition very similar to that of Pinto et al. (2011) : the mass loss comes mainly from the coronal holes located at the poles near the equatorial streamer. However, since the reversals of the two hemisphere happen at different times but still close to one another, we can see that there are less coronal holes than expected that form at year 5.1 (which is almost the multipolar maximum). The minimal value of the mass loss corresponds to a configuration where the polar coronal holes are closing and yet no coronal hole associated to a pseudo-streamer has opened. The maximal value of the mass loss corresponds to a configuration where the polar coronal holes are bigger than at the minimum of activity due to the equatorial streamer which is not completely formed yet.

Finally, our angular momentum loss varies from g cm2 s-1 at minimum of activity to g cm2 s-1 at maximum of activity, hence a variation of 70%. This is the same trend as found in Pinto et al. (2011) or Réville et al. (2015a) : the star loses less angular momentum when it is more multipolar, as shown by the anti-correlation with the parity factor. It can be used to compute the magnetic spin-down timescale, defined as :


where is the Sun’s angular momentum, estimated at g cm2 s-1 in Pinto et al. (2011) from a 1D seismic solar model (Brun et al., 2002). This yields a magnetic spin-down timescale varying from years to years, which is of the same order of magnitude of what was found in Pinto et al. (2011).

Figure 11: Time evolution of the north-south asymmetry over an 11-year cycle for the average Alfvén radius (first panel), mass loss (second panel) and angular momentum loss (third panel). The purple (orange) lines indicate the beginning and end of the field reversal in the southern (northern) hemisphere.

To finish this section, we study the impact of north-south asymmetry on the variations of these global quantities. In Figure 11, we plot the time evolution of the north-south relative difference for the average Alfvén radius, the mass loss and the angular momentum loss in percentage. To compute these quantities for the northern or southern hemisphere only, we use the same definitions as in (20), (21) and (22), except that the integral goes from 0 to or from to . For the average Alfvén radius, the asymmetry goes up to 8%. It gets bigger around year 2.5 and 6.5, which correspond to right before opening and right after closing of coronal holes of opposite polarities, as shown by the purple (orange) lines for the southern (northern) hemisphere. For the mass loss, due to the fact that , we would not expect any asymmetry between the two hemispheres. However, since the wind velocity is influenced by the magnetic field, the asymmetry propagates and is enough to change the mass flux of the hemispheres. The asymmetry is nonetheless less significant that for other quantities, reaching only 3% at most around year 4, which correspond to opening of coronal holes of opposite polarity. Finally it is for the angular momentum loss that the asymmetry is actually more visible, reaching up to a difference of 20% between the two hemispheres. Here, the asymmetry is stronger around year 2.5 and year 4.5. To conclude, in these simulations, the asymmetry is also present in global quantities and is more visible at opening and closing of coronal holes of opposite polarity for each hemisphere, which correspond to the reversal of the field. The fact that the asymmetry is stronger during maximum of activity has been observed for the Sun (Tokumaru et al., 2015).

6 Conclusion and Perspectives

We have studied in this work how the properties of the solar dynamo cycle influence its wind and corona. We followed a similar approach as in Pinto et al. (2011), but went beyond by using more realistic models and considering the effect of north-south asymmetry in the dynamo field. We have used a magnetic field cycle produced by a dynamo code based on a Babcock-Leighton flux-transport model described in Jouve & Brun (2007) and including an asymmetry parameter as in DeRosa et al. (2012). This allowed us to recover solar-like features such as the phase-lag between the poloidal and toroidal field or the relative extent of the equatorward and poleward branches seen in the butterfly diagram (Figure 1). The north-south asymmetry in the magnetic field, that we took into account by introducing an asymmetry in the Babcock-Leighton term, couples the symmetric and antisymmetric family modes of the dynamo, resulting in a delay between the polarity reversals of the northern and southern hemisphere of the star as observed in the Sun. We recall that this model was not parameterized to reproduce a specific solar cycle, our aim was to understand from a theoretical point of view which physical ingredients allow us to recover features similar to observations of the solar corona. We then coupled the output of the dynamo model to a polytropic wind model (note that in Pinto et al. 2011 the wind was considered isothermal) in spherical geometry (adapted from Réville et al., 2015a). We computed 54 states of relaxed wind over an 11-year dynamo cycle : we studied the influence of a complex magnetic topology with north-south asymmetry on the wind.

We find that an anti-correlation between the energy of the dominant magnetic mode and the total magnetic energy at the surface tends to reduce the variations of the global quantities, as in Réville & Brun (2017). In Pinto et al. (2011), the Alfvén radius varies by a factor 4 and the mass loss varies by 60% over the cycle ; in our study, we find that the Alfvén radius varies by a factor 2 and the mass loss by 13%, which is closer to what was found in Réville & Brun (2017) based on magnetograms of the Sun, where the Alfvén radius varies by 30% and the mass loss by 20%.

The variations of the different magnetic modes over the cycle induce a complex corona with very different features at minimum and maximum. At minimum we have faster wind at the poles and slower wind at the equator ; the wind is very organized with a main streamer at the equator separating the positive from the negative polarity regions and coronal holes are located only at the poles. At maximum the distribution of the wind speed loses a clear latitudinal dependency, similarly to the Sun as observations with Ulysses clearly demonstrated (McComas et al., 2008). Along the cycle more and more pseudo-streamers and streamers emerge, with coronal holes of opposite polarity opening first at mid-latitudes, first in the southern hemisphere and then in the northern, leading to a complex latitudinal polarity repartition with up to 6 different regions of different polarity next to one another at maximum. The mass loss originates from coronal holes and close to the streamers : as the corona evolves over the dynamo cycle, so does the mass loss in latitude and amplitude. All these elements show that at a given latitude, the dynamics of the corona change drastically over a dynamo cycle.

In the dynamo model studied here, the dipole and quadrupole were not necessarily the strongest mode in their mode family. This effect is particularly visible at the poles, where higher modes tend to dominate the magnetic field. This influences the delay between the minimum of sunspot activity and dipolar activity, resulting in a different mass loss profile as what is usually computed (see, e.g. Pinto et al., 2011). This however allows us to really pinpoint the interplay between the different modes and the solar wind. As shown in Figure 9, the Alfvén radius and the angular momentum loss rate are anti-correlated with the parity factor, because they are highly sensitive to topology, and which also indicates that only the large-scale modes have a significant impact on the lever arm (e.g. Alfvén radius) for the braking of the star. On the other hand in our model, the mass loss rate is much less sensitive to topology, it mainly follows the global magnetic surface energy. We recall that this result may depend on the fact that we used basic heating for the corona, this needs to be confirmed for more realistic solar atmospheres. Nevertheless we stress here the importance of the detailed energy repartition between the scales and the relative phase of the two symmetry families of the dynamo to assess the quality of a given dynamo solution when compared to the Sun.

The asymmetry in the dynamo model also allowed us to recover more realistic profiles for the wind distribution along the cycle : in Figure 7 we obtain qualitatively the same structure as displayed in Sokół et al. (2015) using IPS data from Tokumaru et al. (2010), with a time-latitude variation correlated with the dynamo cycle. The asymmetry is clearly visible, with the southern hemisphere wind speed slowing down first in this particular example. The transition between maximum and minimum of activity is as smooth as in the solar observations. It is also worth noting that the north-south asymmetry was a 9-month delay for the dynamo generated field, and translated into a 5-month delay for the wind speed at 0.1 AU. The wind seems thus to smooth the asymmetry. Some studies indicate even that, given the proper range of parameters, the asymmetry can slightly appear naturally in flux-transport models if there is a non-linear coupling (Shukuya & Kusano, 2017). It could be interesting to do a more systematic study of the influence of the delay between the two hemisphere reversals. In our case, the delay is shorter than in the Sun. The asymmetry is also present in global quantities such as the Alfvén radius, the mass loss and the angular momentum. If this result can be applied to the Sun, it means that we have to be take it into account when inferring theses global quantities from measures taken at a single latitude, especially for the angular momentum loss.

We have to bear in mind the limitations of our model. The fact that the dynamo was not calibrated to be exactly solar-like has for direct consequence that all physical quantities may not necessarily exhibit a solar-like behavior. For instance the polytropic approach is a simple approximation of the heating of the corona, which leads to e.g. a smaller Alfvén radius than the expected solar Alfvén radius. We could have calibrated our wind model in order to increase it, but this would have also modified other physical quantities. Instead, we chose to calibrate our model to reproduce the mean mass loss rate of the Sun for it is a better known value deduced from in situ observations. Finally we recall that with such a quasi-static approach, we can understand some of the influence of the dynamo on the wind, but we cannot include the influence of the wind on the dynamo, if any. Still our results show interesting relation between dynamo and wind properties.

For the future, several aspects remain to be explored. A more realistic coronal heating such as the ones prescribed in Cranmer et al. (2007) or Riley et al. (2015) should be use to improve the realism of the wind solutions. Réville & Brun (2017) have shown that considering the full 3D topology of the solar magnetic field also brings a lot more information and can help recovering more realistic features of the solar environment. For example, an important feature when comparing to observational data is the tilt of the heliospheric current sheet, which is omitted here due to the assumption of axisymmetry. We are currently developing a 3D set-up up to 1 AU to look at such tendencies. The north-south asymmetry we studied was also not variable in time, which should be the case in the Sun: we could add stochastic noise in the future to make it even more realistic on longer timescales and thus obtain different delays depending on the cycle. Please note that the two codes used can also be used to model other stars than the Sun, a similar study could be performed for any star on the main sequence. Finally, the best way to understand the full interplay between the dynamo and the wind would be to consider a dynamical coupling, to avoid any non-causal perturbations, which we intend to explore in a future work, along with some of the improvements listed above.


We thank R. Pinto for useful discussions. This work was supported by a CEA "Thèse Phare" grant, by CNRS and INSU/PNST program and by CNES SHM funds. Computations were carried out using CEA CCRT and CNRS IDRIS facilities within the GENCI 20410133 allocation.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

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

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