1 Introduction
From a space weather perspective, one of the main challenges is to accurately model 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. (Reference McComas, Ebert, Elliott, Goldstein, Gosling, Schwadron and Skoug2008) and Smith (Reference Smith2011) for a summary of the mission highlights, and Issautier et al. (Reference Issautier, Le Chat, Meyer-Vernet, Moncuquet, Hoang, MacDowall and McComas2008) for a discussion about north-south asymmetry). At the minimum of activity, the magnetic field is mostly dipolar and the wind is slower at the equator (approximately
$400~\text{km}~\text{s}^{-1}$
) and faster at the poles (approximately
$800~\text{km}~\text{s}^{-1}$
); at the maximum of activity, the magnetic field is multipolar and the wind distribution is bimodal at all latitudes. This is why there is an increasing 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’s magnetic activity: sunspots and magnetic field observations from several observatories and several time periods (Royal Greenwich Observatory from 1874 to 1954 in Newton & Milsom (Reference Newton and Milsom1955) and from 1874 to 1976 in Vizoso & Ballester (Reference Vizoso and Ballester1990), Wilcox Solar Observatory from 1976 to 2009 in Hoeksema (Reference Hoeksema, Kosovichev, Andrei and Rozelot2010)) have shown a north-south asymmetry in the sunspot 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 Reference Svalgaard and Kamide2013). 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. Reference Temmer, Rybák, Bendík, Veronig, Vogler, Otruba, Pötzi and Hanslmeier2006). 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 Reference Mursula and Hiltula2003).
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 which can lead in extreme cases to the vanishing of the magnetic field in one hemisphere (Gallet & Pétrélis Reference Gallet and Pétrélis2009). See also Tobias (Reference Tobias1997), Ossendrijver (Reference Ossendrijver2003), DeRosa, Brun & Hoeksema (Reference DeRosa, Brun and Hoeksema2012), Shukuya & Kusano (Reference Shukuya and Kusano2017) and references therein for a detailed discussion on dynamo symmetry properties. The antisymmetric family corresponds to odd
$\ell$
degrees when projecting the magnetic field onto the spherical harmonics functions
$Y_{\ell }^{m}$
under the assumption of axisymmetry (with order
$m=0$
). 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
$\ell$
degrees (when
$m=0$
), is not negligible; it reaches on average 25 % of the amplitude of the antisymmetric family, and becomes dominant during polarity reversals. Such asymmetry also has an impact on the wind structure: in Sokół et al. (Reference Sokół, Swaczyna, Bzowski and Tokumaru2015), the reconstruction of the solar wind using interplanetary scintillation (IPS) observations summed up in Tokumaru, Kojima & Fujiki (Reference Tokumaru, Kojima and Fujiki2010), 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 Reference Parker1993): due to strong shears in the solar differential rotation (Schou et al. Reference Schou, Antia, Basu, Bogart, Bush, Chitre, Christensen-Dalsgaard, Di Mauro, Dziembowski and Eff-Darwich1998), 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 (Reference Miesch2005) for a 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 (Reference Babcock1961) and Leighton (Reference Leighton1969), but also Wang & Sheeley (Reference Wang and Sheeley1991) for the link with observations and Dikpati & Charbonneau (Reference Dikpati and Charbonneau1999) for simulations).
To simulate the whole Sun, magnetohydrodynamic (MHD) simulations have been used successfully since the 1970s (for a full review on the subject, see Brun & Browning (Reference Brun and Browning2017)).
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 were first computed in Gilman (Reference Gilman1983) and in Glatzmaier (Reference Glatzmaier1985), and first adapted for high resolutions in Brun, Miesch & Toomre (Reference Brun, Miesch and Toomre2004); more realistic set-ups with high Reynolds numbers can be found in Hotta, Rempel & Yokoyama (Reference Hotta, Rempel and Yokoyama2016). LES models are now able to display large-scale magnetic cycles: for the Sun see Ghizaru, Charbonneau & Smolarkiewicz (Reference Ghizaru, Charbonneau and Smolarkiewicz2010), for young convective stars see Brown et al. (Reference Brown, Miesch, Browning, Brun and Toomre2011), for exceptional cycles such as grand minima see Augustson et al. (Reference Augustson, Brun, Miesch and Toomre2015) and for generalization to other solar-like stars see Strugarek et al. (Reference Strugarek, Beaudoin, Charbonneau, Brun and do Nascimento2017).
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 three-dimensional (3-D) MHD methods. A presentation of kinematic dynamo models can be found in Roberts (Reference Roberts1972), for more focus on the role of rotation see Stix (Reference Stix1976); for general reviews see Krause & Raedler (Reference Krause and Raedler1980) and Charbonneau (Reference Charbonneau2010). Although 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. Reference Kumar, Jouve, Pinto and Rouillard2018).
For the solar wind, the current framework has been initiated by the work of Parker (Reference Parker1958): 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 Reference Schatzman1962; Weber & Davis Reference Weber and Davis1967). 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 Reference Hollweg and Isenberg2002). These 1-D analytical equatorial models were then expanded to two dimensions (Sakurai Reference Sakurai1985), 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 modelled. 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 Reference Parker1988) and MHD turbulence (Cranmer, van Ballegooijen & Edgar Reference Cranmer, van Ballegooijen and Edgar2007). 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. Reference Linker, Caplan, Downs, Riley, Mikic, Lionello, Henney, Liu, Derosa and Yeates2017).
One way to reproduce the wind dynamics is via compressible MHD simulations (Keppens & Goedbloed Reference Keppens and Goedbloed1999) 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 Reference Suzuki and Inutsuka2005) or two dimensions (Matsumoto & Suzuki Reference Matsumoto and Suzuki2012), or more global simulations including the whole star either in 2.5 dimensions in Matt et al. (Reference Matt, MacGregor, Pinsonneault and Greene2012), Tóth et al. (Reference Tóth, van der Holst, Sokolov, De Zeeuw, Gombosi, Fang, Manchester, Meng, Najib and Powell2012), Réville et al. (Reference Réville, Brun, Matt, Strugarek and Pinto2015a ) or in three dimensions in Riley et al. (Reference Riley, Lionello, Linker, Cliver, Balogh, Beer, Charbonneau, Crooker, DeRosa and Lockwood2015) and Réville & Brun (Reference Réville and Brun2017).
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 altogether make the modelling 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 megametres) including the different photospheric layers (Vögler et al. (Reference Vögler, Shelyag, Schüssler, Cattaneo, Emonet and Linde2005), Martínez-Sykora, Hansteen & Carlsson (Reference Martínez-Sykora, Hansteen and Carlsson2008) for the addition of conduction, Rempel, Schüssler & Knölker (Reference Rempel, Schüssler and Knölker2009) for focus on sunspots among others; see also review by Wedemeyer-Böhm, Lagg & Nordlund (Reference Wedemeyer-Böhm, Lagg and Nordlund2009)).
In a previous study, published in Pinto et al. (Reference Pinto, Brun, Jouve and Grappin2011), the influence of the magnetic field geometry and amplitude on the solar wind has been studied in a quasi-static way in 2.5 dimensions. A first numerical code named STELEM (STellar ELEMents, cf. Jouve & Brun (Reference Jouve and Brun2007)) computes an 11-year magnetic cycle using a kinematic mean-field
$\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FA}$
dynamo model. The latter is used as a bottom boundary condition for a second numerical code (DIP, cf. Grappin et al. (Reference Grappin, Léorat, Leygnac and Pinto2010)) 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 (Reference Jouve and Brun2007), and emphasizing north–south asymmetry as in DeRosa et al. (Reference DeRosa, Brun and Hoeksema2012) (see their appendix A). We will first present the numerical codes and the physical ingredients in our models in § 2, then we will present our results as follows: we will begin with a description of the magnetic dynamo cycle in § 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 § 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 § 5. Discussion and conclusions are given in § 6.
2 Numerical set-up
We used two different MHD codes for our 2.5-D axisymmetric simulations: STELEM for the dynamo field generated inside the star with a flux-transport model (T. Emonet & P. Charbonneau 1998, private communication; Jouve & Brun (Reference Jouve and Brun2007)), and PLUTO for the solar corona and the associated wind (Mignone et al. Reference Mignone, Bodo, Massaglia, Matsakos, Tesileanu, Zanni and Ferrari2007). 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 § 2.1 we describe physical properties and the numerical methods used in the STELEM code; the same is done for the PLUTO code in § 2.2. Finally we give further details about the coupling process in § 2.3.
2.1 Dynamo simulations with STELEM
The simulations performed are the same as that described in DeRosa et al. (Reference DeRosa, Brun and Hoeksema2012) in their appendix.
Working in spherical coordinates
$(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$
and under the assumption of axisymmetry, we perform a poloidal–toroidal decomposition and write the mean magnetic and velocity field (respectively
$\boldsymbol{B}$
and
$\boldsymbol{U}$
) as follows:
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}(r,\unicode[STIX]{x1D703},t)=\unicode[STIX]{x1D735}\times (A_{\unicode[STIX]{x1D719}}(r,\unicode[STIX]{x1D703},t)\boldsymbol{e}_{\unicode[STIX]{x1D719}})+B_{\unicode[STIX]{x1D719}}(r,\unicode[STIX]{x1D703},t)\boldsymbol{e}_{\unicode[STIX]{x1D719}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{U}(r,\unicode[STIX]{x1D703})=\boldsymbol{u}_{p}(r,\unicode[STIX]{x1D703})+r\sin \unicode[STIX]{x1D703}\unicode[STIX]{x1D6FA}(r,\unicode[STIX]{x1D703})\boldsymbol{e}_{\unicode[STIX]{x1D719}}. & \displaystyle\end{eqnarray}$$
$\boldsymbol{B}$
is decomposed with the poloidal streamfunction
$A_{\unicode[STIX]{x1D719}}$
and toroidal field
$B_{\unicode[STIX]{x1D719}}$
. The velocity field is time independent and prescribed by profiles of the meridional circulation
$\boldsymbol{u}_{p}$
and differential rotation
$\unicode[STIX]{x1D6FA}$
.
We rewrite the induction equation in the framework of mean-field theory in terms of
$A_{\unicode[STIX]{x1D719}}$
and
$B_{\unicode[STIX]{x1D719}}$
and we introduce a poloidal term
$S$
based on the Babcock–Leighton mechanism. We finally normalize the equations by using the solar radius
$R_{\odot }$
as the length scale and the diffusion time
$R_{\odot }^{2}/\unicode[STIX]{x1D702}_{t}$
as the time scale, and obtain the following two coupled partial differential equations:
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}A_{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}_{t}}\left(\unicode[STIX]{x1D735}^{2}-\frac{1}{\unicode[STIX]{x1D71B}^{2}}\right)A_{\unicode[STIX]{x1D719}}-R_{e}\frac{\boldsymbol{u}_{p}}{\unicode[STIX]{x1D71B}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D71B}A_{\unicode[STIX]{x1D719}})+C_{s}S, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}B_{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}t} & = & \displaystyle \frac{\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D702}_{t}}\left(\unicode[STIX]{x1D735}^{2}-\frac{1}{\unicode[STIX]{x1D71B}^{2}}\right)B_{\unicode[STIX]{x1D719}}+\frac{1}{\unicode[STIX]{x1D71B}}\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D71B}B_{\unicode[STIX]{x1D719}})}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D702}/\unicode[STIX]{x1D702}_{t})}{\unicode[STIX]{x2202}r}-R_{e}\unicode[STIX]{x1D71B}\boldsymbol{u}_{p}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{B_{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x1D71B}}\right)\nonumber\\ \displaystyle & & \displaystyle -\,R_{e}B_{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{p}+C_{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D71B}[\unicode[STIX]{x1D735}\times (A_{\unicode[STIX]{x1D719}}\boldsymbol{e}_{\unicode[STIX]{x1D719}})]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FA},\end{eqnarray}$$
where
$\unicode[STIX]{x1D702}$
is the effective magnetic diffusivity,
$\unicode[STIX]{x1D702}_{t}$
is the turbulent diffusivity in the convection zone and
$\unicode[STIX]{x1D71B}=r\sin \unicode[STIX]{x1D703}$
. These equations are controlled by three dimensionless parameters:
$C_{\unicode[STIX]{x1D6FA}}=\unicode[STIX]{x1D6FA}_{0}R_{\odot }^{2}/\unicode[STIX]{x1D702}_{t}$
which characterizes the shear by differential rotation (i.e. the omega effect);
$C_{s}=s_{0}R_{\odot }/\unicode[STIX]{x1D702}_{t}$
which characterizes the Babcock–Leighton source term;
$R_{e}=u_{0}R_{\odot }/\unicode[STIX]{x1D702}_{t}$
(the Reynolds number) which characterizes the intensity of the meridional circulation.
$\unicode[STIX]{x1D6FA}_{0}$
,
$s_{0}$
and
$u_{0}$
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
$C_{\unicode[STIX]{x1D6FA}}=1.4\times 10^{5}$
,
$C_{s}=30$
and
$R_{e}=1.20\times 10^{3}$
, which corresponds to the parameters in DeRosa et al. (Reference DeRosa, Brun and Hoeksema2012). The rotation rate thus corresponds 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 depth, vanishing at the bottom radial boundary. The equatorward flow penetrates slightly beneath the tachocline (Dikpati & Charbonneau Reference Dikpati and Charbonneau1999).
The rotation profile is inspired by solar angular velocity profile deduced from helioseismic inversions (Thompson et al.
Reference Thompson, Christensen-Dalsgaard, Miesch and Toomre2003). We assume solid-body rotation below
$r=0.66R_{\odot }$
and a differential rotation above:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}(r,\unicode[STIX]{x1D703})=\unicode[STIX]{x1D6FA}_{c}+\frac{1}{2}\left[1+\text{erf}\left(\frac{2(r-r_{c})}{d_{1}}\right)\right]\left(\unicode[STIX]{x1D6FA}_{eq}+a_{2}\cos ^{2}\unicode[STIX]{x1D703}+a_{4}\cos ^{4}\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FA}_{c}\right), & & \displaystyle\end{eqnarray}$$
with the parameters:
$\unicode[STIX]{x1D6FA}_{eq}=1$
,
$\unicode[STIX]{x1D6FA}_{c}=0.93944$
,
$r_{c}=0.7R_{\odot }$
,
$d_{1}=0.05R_{\odot }$
,
$a_{2}=-0.136076$
and
$a_{4}=-0.145713$
.
We assume different diffusivities in the envelope and in the stable interior, smoothly matching the two:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}(r)=\unicode[STIX]{x1D702}_{c}+\frac{(\unicode[STIX]{x1D702}_{t}-\unicode[STIX]{x1D702}_{c})}{2}\left[1+\text{erf}\left(\frac{r-r_{c}}{d}\right)\right], & & \displaystyle\end{eqnarray}$$
with the parameters:
$\unicode[STIX]{x1D702}_{c}=10^{9}~\text{cm}^{2}~\text{s}^{-1}$
,
$\unicode[STIX]{x1D702}_{t}=10^{11}~\text{cm}^{2}~\text{s}^{-1}$
and
$d=0.03R_{\odot }$
.
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
$B_{\unicode[STIX]{x1D719}}$
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
$\unicode[STIX]{x1D716}$
:
$$\begin{eqnarray}\displaystyle S(r,\unicode[STIX]{x1D703},B_{\unicode[STIX]{x1D719}}) & = & \displaystyle \frac{1}{2}\left[1+\text{erf}\left(\frac{r-r_{2}}{d_{2}}\right)\right]\left[1-\text{erf}\left(\frac{r-R_{\odot }}{d_{2}}\right)\right]\left[1+\left(\frac{B_{\unicode[STIX]{x1D719}}(r_{c},\unicode[STIX]{x1D703},t)}{B_{0}}^{2}\right)\right]^{-1}\nonumber\\ \displaystyle & & \displaystyle \times \,(\cos \unicode[STIX]{x1D703}+\unicode[STIX]{x1D716}\sin \unicode[STIX]{x1D703})\sin ^{3}\unicode[STIX]{x1D703}B_{\unicode[STIX]{x1D719}}(r_{c},\unicode[STIX]{x1D703},t),\end{eqnarray}$$
with the parameters:
$r_{2}=0.95R_{\odot }$
,
$d_{2}=0.01R_{\odot }$
,
$B_{0}=10^{5}$
and
$\unicode[STIX]{x1D716}=10^{-3}$
.
The STELEM code is used to solve (2.3) and (2.4): it uses a finite-element method in space and a third-order scheme in time. The temporal scheme used is adapted from Spalart, Moser & Rogers (Reference Spalart, Moser and Rogers1991) 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. Reference Jouve, Brun, Arlt, Brandenburg, Dikpati, Bonanno, Käpylä, Moss, Rempel and Gilman2008).
The numerical domain is an annular meridional cut of the Sun with the colatitude
$\unicode[STIX]{x1D703}\in [0,\unicode[STIX]{x03C0}]$
and the radius
$r\in [0.6,1]R_{\odot }$
(i.e. from slightly below the tachocline located at
$r\approx 0.7R_{\odot }$
up to the solar surface). We use a uniform grid with 64 points in radius and cosine of the latitude. At the latitudinal boundaries (
$\unicode[STIX]{x1D703}=0$
and
$\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}$
) and at the bottom radial boundary (
$r=0.6R_{\odot }$
),
$A_{\unicode[STIX]{x1D719}}$
and
$B_{\unicode[STIX]{x1D719}}$
are set to 0. At the upper radial boundary (
$r=R_{\odot }$
), the solution is matched to an external potential field. Usual initial conditions involve setting a confined dipole field configuration, i.e.
$A_{\unicode[STIX]{x1D719}}$
is set to
$\sin \unicode[STIX]{x1D703}/r^{2}$
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 that described in Réville et al. (Reference Réville, Brun, Matt, Strugarek and Pinto2015a ), but adapted to spherical coordinates.
We solve the set of the conservative ideal MHD equations composed of the continuity equation for the density
$\unicode[STIX]{x1D70C}$
, the momentum equation for the velocity field
$\boldsymbol{u}$
with its momentum written
$\boldsymbol{m}=\unicode[STIX]{x1D70C}\boldsymbol{u}$
, the energy equation which is noted
$E$
and the induction equation for the magnetic field
$\boldsymbol{B}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D70C}\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\boldsymbol{m}+\unicode[STIX]{x1D6FB}\boldsymbol{\cdot }(\boldsymbol{m}\boldsymbol{v}-\boldsymbol{B}\boldsymbol{B}+\mathsf{I}p)=\unicode[STIX]{x1D70C}\boldsymbol{a}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}E+\unicode[STIX]{x1D735}\boldsymbol{\cdot }((E+p)\boldsymbol{u}-\boldsymbol{B}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{B}))=\boldsymbol{m}\boldsymbol{\cdot }\boldsymbol{a}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\boldsymbol{B}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{B}-\boldsymbol{B}\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$
where
$p$
is the total pressure (thermal and magnetic),
$\mathsf{I}$
is the identity matrix and
$\boldsymbol{a}$
is a source term (gravitational acceleration in our case). We use the ideal equation of state:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}\unicode[STIX]{x1D700}=p_{th}/(\unicode[STIX]{x1D6FE}-1), & & \displaystyle\end{eqnarray}$$
where
$p_{th}$
is the thermal pressure,
$\unicode[STIX]{x1D700}$
is the internal energy per mass and
$\unicode[STIX]{x1D6FE}$
is the adiabatic exponent. This gives for the energy:
$E=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D700}+\boldsymbol{m}^{2}/(2\unicode[STIX]{x1D70C})+\boldsymbol{B}^{2}/2$
.
PLUTO solves normalized equations, using three variables to set all the others: length, density and speed. If we note with
$\ast$
the parameters related to the star and with 0 the parameters related to the normalization, we have
$R_{\ast }/R_{0}=1$
,
$\unicode[STIX]{x1D70C}_{\ast }/\unicode[STIX]{x1D70C}_{0}=1$
and
$u_{kep}/U_{0}=\sqrt{GM_{\ast }/R_{\ast }}/U_{0}=1$
, where
$u_{kep}$
is the Keplerian speed at the stellar surface and
$G$
the gravitational constant. By choosing the physical values of
$R_{0}$
,
$\unicode[STIX]{x1D70C}_{0}$
and
$U_{0}$
, one can deduce all of the other values given by the code in physical units. In our set-up, we choose
$R_{0}=R_{\odot }=6.96\times 10^{10}$
cm,
$\unicode[STIX]{x1D70C}_{0}=\unicode[STIX]{x1D70C}_{\odot }=6.68\times 10^{-16}~\text{g}~\text{cm}^{-3}$
(which corresponds to the density in the solar corona above 2.5 km, cf. Vernazza, Avrett & Loeser (Reference Vernazza, Avrett and Loeser1981)) and
$U_{0}=u_{kep,\odot }=4.37\times 10^{2}~\text{km}~\text{s}^{-1}$
. Our wind simulations are then controlled by three parameters: the adiabatic exponent
$\unicode[STIX]{x1D6FE}$
for the polytropic wind, the rotation of the star normalized by the escape velocity
$u_{rot}/u_{esc}$
and the speed of sound normalized also by the escape velocity
$c_{s}/u_{esc}$
. Note that the escape velocity is defined as
$u_{esc}=\sqrt{2}u_{kep}=\sqrt{2GM_{\ast }/R_{\ast }}$
. For the rotation speed, we take the solar value, which gives
$u_{rot}/u_{esc}=2.93\times 10^{-3}$
. We choose to fix
$c_{s}/u_{esc}=0.243$
, which corresponds to a
$1.3\times 10^{6}$
K hot corona for solar parameters and
$\unicode[STIX]{x1D6FE}=1.05$
. This choice of
$\unicode[STIX]{x1D6FE}$
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
$\unicode[STIX]{x1D6FE}\neq 5/3$
is a simplified way of taking into account heating, which is not modelled here.
We assume axisymmetry and use the spherical coordinates
$(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$
. 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 (Reference Einfeldt1988)). PLUTO uses a reconstruct–solve–average approach using a set of primitive variables
$(\unicode[STIX]{x1D70C},\boldsymbol{u},p,\boldsymbol{B})$
to solve the Riemann problem corresponding to the previous set of equations.
The numerical domain is an annular meridional cut with the colatitude
$\unicode[STIX]{x1D703}\in [0,\unicode[STIX]{x03C0}]$
and the radius
$r\in [1,20]R_{\odot }$
. We use a uniform grid in latitude with 512 points, and a stretched grid in radius with 512 points; the grid spacing is geometrically increasing from
$\unicode[STIX]{x0394}r/R_{\ast }=0.001$
at the surface of the star to
$\unicode[STIX]{x0394}r/R_{\ast }=0.01$
at the outer boundary. At the latitudinal boundaries (
$\unicode[STIX]{x1D703}=0$
and
$\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}$
), we set axisymmetric boundary conditions. At the top radial boundary (
$r=20R_{\ast }$
), we set an outflow boundary condition which corresponds to
$\unicode[STIX]{x2202}/\unicode[STIX]{x2202}r=0$
for all variables, except for the radial magnetic field where we enforce
$\unicode[STIX]{x2202}(r^{2}B_{r})/\unicode[STIX]{x2202}r=0$
. 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 (
$r=R_{\ast }$
), we set a condition similar to the one described in Zanni & Ferreira (Reference Zanni and Ferreira2009) to be as close as possible to a perfect rotating conductor. We also implement the same differential rotation as described in (2.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.
Reference Dedner, Kemm, Kröner, Munz, Schnitzer and Wesenberg2002).
2.3 Coupling method
To couple the two codes described above, we use the projection of the surface magnetic field produced by STELEM onto spherical harmonics, which PLUTO can read. To proceed so, we select a given time
$t_{0}$
, then we project the surface radial magnetic field at that time
$B_{r}(t_{0},r=R_{\odot },\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$
onto spherical harmonics
$Y_{\ell }^{0}(\unicode[STIX]{x1D703})$
up to a degree
$\ell _{max}$
:
$$\begin{eqnarray}\displaystyle B_{r}(r=R_{\odot },\unicode[STIX]{x1D703},t_{0})=\mathop{\sum }_{\ell =1}^{\ell _{max}}\unicode[STIX]{x1D6FC}_{\ell ,0}(r=R_{\odot })Y_{\ell }^{0}(\unicode[STIX]{x1D703}). & & \displaystyle\end{eqnarray}$$
Note that there is no dependency in
$\unicode[STIX]{x1D719}$
due to axisymmetry (
$m=0$
). In this study, we choose
$\ell _{max}=40$
, 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
$R_{ss}$
equal to
$15R_{\odot }$
. For more information about the PFSS method, for the original computation see Altschuler & Newkirk (Reference Altschuler and Newkirk1969) and Schatten, Wilcox & Ness (Reference Schatten, Wilcox and Ness1969), for a more explicit calculation see Schrijver & De Rosa (Reference Schrijver and De Rosa2003). We obtain:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}_{\ell ,0}(r)=\unicode[STIX]{x1D6FC}_{\ell ,0}(r=R_{\odot })\frac{\ell (R_{\odot }/R_{ss})^{2\ell +1}(r/R_{\odot })^{\ell -1}+(\ell +1)(r/R_{\odot })^{-(\ell +2)}}{\ell (R_{\odot }/R_{ss})^{2\ell +1}+(\ell +1)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}_{\ell ,0}(r)=(\ell +1)\unicode[STIX]{x1D6FC}_{\ell ,0}(r=R_{\odot })\frac{(R_{\odot }/R_{ss})^{2\ell +1}(r/R_{\odot })^{\ell -1}-(r/R_{\odot })^{-(\ell +2)}}{\ell (R_{\odot }/R_{ss})^{2\ell +1}+(\ell +1)}. & \displaystyle\end{eqnarray}$$
In Réville et al. (Reference Réville, Brun, Strugarek, Matt, Bouvier, Folsom and Petit2015b
), it was explained that the fiducial value of 2.5
$R_{\odot }$
usually found in the literature for
$R_{ss}$
was an approximation, and that it is possible to define an optimal source surface radius that allows recovery of the best value for the open flux in the 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
$R_{ss}$
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:
$$\begin{eqnarray}\displaystyle \boldsymbol{B}^{PLUTO}=\left\{\begin{array}{@{}l@{}}\displaystyle B_{r}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\mathop{\sum }_{\ell =1}^{\ell _{max}}\unicode[STIX]{x1D6FC}_{\ell ,0}(r)Y_{\ell }^{0}(\unicode[STIX]{x1D703})\\ \displaystyle B_{\unicode[STIX]{x1D703}}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\mathop{\sum }_{\ell =1}^{\ell _{max}}\unicode[STIX]{x1D6FD}_{\ell ,0}(r)\sqrt{\frac{2\ell +1}{4\unicode[STIX]{x03C0}}}\frac{1}{\ell +1}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}P_{\ell }^{0}(\unicode[STIX]{x1D703})\\ \displaystyle B_{\unicode[STIX]{x1D719}}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=0.0.\end{array}\right. & & \displaystyle\end{eqnarray}$$
The magnetic field provided by STELEM is dimensionless (cf. (2.3) and (2.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.
Reference McComas, Ebert, Elliott, Goldstein, Gosling, Schwadron and Skoug2008); this yields values between
$2.3\times 10^{-14}$
and
$3.1\times 10^{-14}~M_{\odot }/\text{yr}$
(Réville & Brun Reference Réville and Brun2017). The magnetic field is then normalized by the PLUTO units described in § 2.2.
From our dynamo run, we thus obtain a times series of 54 snapshots
$B_{r}(t_{i})$
, sampled over an 11-year cycle with a time step of approximately 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
$t_{P}=R_{\ast }/u_{kep}=R_{\ast }^{3/2}/\sqrt{GM_{\ast }}$
), 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
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 with the latitude. (a) Shows the general aspect of the magnetic field over 3 cycles. (b) 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 1(a) 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 Reference Clette and Lefèvre2012). 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, as in cycle 18 of the Sun). This is in qualitative agreement with the current observed delay of 1 year (Temmer et al. Reference Temmer, Rybák, Bendík, Veronig, Vogler, Otruba, Pötzi and Hanslmeier2006). 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. Figure 1(b) 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 § 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 (SSN) based on the generation of toroidal magnetic field at the tachocline, given by:
$$\begin{eqnarray}\displaystyle SSN_{proxy}=SSN_{0}\int _{\unicode[STIX]{x1D703}=0}^{\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}}B_{\unicode[STIX]{x1D719}}(r_{c},\unicode[STIX]{x1D703},t)^{2}r_{c}^{2}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}, & & \displaystyle\end{eqnarray}$$
where
$SSN_{0}$
is a scaling constant to find values for our SSN proxy of the same order of magnitude as the Sun (Hung et al.
Reference Hung, Brun, Fournier, Jouve, Talagrand and Zakari2017).
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.
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.
Time evolution of the coefficients of the surface magnetic field projection onto spherical harmonics over 2.5 dynamo cycles, grouped as equatorially symmetric and antisymmetric families. Only the
$m=0$
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
$\ell$
for the projection onto the spherical harmonics when considering only
$m=0$
modes. The first thing we notice is that the amplitude of the antisymmetric family modes is similar to that of the Sun (between
$-$
4 and 4 G, as shown in DeRosa et al. (Reference DeRosa, Brun and Hoeksema2012)). The
$\ell =3$
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 favouring higher
$\ell$
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 approximately 10 % of the antisymmetric family mean amplitude. This is different from simple 2.5-D 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.5 G) is close to what has been observed since cycle 21 (see for instance DeRosa et al. (Reference DeRosa, Brun and Hoeksema2012)).
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. (2.16)), we can define the total radial surface magnetic energy as:
$$\begin{eqnarray}\displaystyle ME=\mathop{\sum }_{\ell =1}^{\ell _{max}}\unicode[STIX]{x1D6FC}_{\ell ,0}^{2}. & & \displaystyle\end{eqnarray}$$
Then we can define the energy of a specific harmonic component of the radial surface field as:
$$\begin{eqnarray}\displaystyle f_{\ell }=\unicode[STIX]{x1D6FC}_{\ell ,0}^{2}/ME. & & \displaystyle\end{eqnarray}$$
Using data from the Wilcox Observatory it can be shown that
$ME$
and
$f_{1}$
were anticorrelated from cycle 20 to 23 (Réville & Brun Reference Réville and Brun2017), 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 figure 4(a). In our case,
$f_{1}$
and
$ME$
have a phase delay of one quarter of a period, so not fully correlated nor anticorrelated 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 figure 4(b), we show the same comparison but for the
$\ell =5$
mode, and here we have a phase delay of one third of a period, which is closer to anticorrelation. 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.
Comparison of the evolution of the surface magnetic energy and the energy of 2 mode components over several cycles:
$\ell =1$
for (a),
$\ell =5$
for (b).
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
$t=0$
, the different panels (a,b,c,d,e) and (f) correspond respectively to times
$t=0.0$
, 2.9, 3.8, 4.5, 6.0 and 11.9 years. These times were indicated as black dashed lines in figure 1(a). The first 4
$R_{\odot }$
close to the star are visible for both hemispheres in order to see the effect of asymmetry. The colour scale represents the following quantity:
$\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{B}/(c_{s}\Vert \boldsymbol{B}\Vert )$
, which is the solar wind velocity projected onto 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 colour table, 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
$A_{\unicode[STIX]{x1D719}}$
.
Snapshots of the evolution of the coronal magnetic field at different times: if we set
$t=0$
years at the first minimum, then we have
$t=0$
, 2.9, 4.2, 4.9, 5.8 and 12.0 years. The colour scale represents the following quantity:
$\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{B}/(c_{s}\Vert \boldsymbol{B}\Vert )$
, which is the solar wind velocity projected onto 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 analyse 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, Sheeley & Rich (Reference Wang, Sheeley and Rich2007).
We start at a minimum of activity, go through a maximum and then return to a minimum. At the minimum of activity (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 (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 (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 (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 (e), the southern hemisphere has now reversed, since the overall magnetic field, especially at the poles, is now the opposite of the one in (a). In the northern hemisphere, the reversal is not complete, the pole is the last region to which the coronal hole of opposite polarity has not yet spread. In (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 (a), hence completing an activity cycle.
Radial cuts at minimum and maximum of activity, for
$u_{r}$
,
$u_{\unicode[STIX]{x1D719}}$
and
$B_{r}$
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
$u_{r}$
,
$u_{\unicode[STIX]{x1D719}}$
and
$B_{r}$
to see if we recover the expected behaviours. 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 the 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 (Reference Keppens and Goedbloed1999). 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 radius 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
$1/r$
dependency. Far from the star,
$B_{r}$
displays an expected dependency of
$1/r^{2}$
, typical of a purely radial field stretched by the wind to satisfy the divergence-free property (in spherical coordinates:
$\unicode[STIX]{x2202}(r^{2}B_{r})/\unicode[STIX]{x2202}r=0$
). This is true both at the 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 affects planets or satellites. Figure 7 shows the evolution of the radial wind velocity at
$r=20R_{\odot }$
in
$\text{km}~\text{s}^{-1}$
over an 11-year cycle in our simulation. We show only the radial velocity because this far from the star, the wind has become mostly radial, so it is
$u_{r}$
which contains most of the wind total velocity. We chose
$r=20R_{\odot }$
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.
Time–latitude diagram of the wind speed over an 11-year cycle in
$\text{km}~\text{s}^{-1}$
at
$r=20R_{\odot }$
.
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~\text{km}~\text{s}^{-1}$
, while here we reach between 200 and
$235~\text{km}~\text{s}^{-1}$
at 0.1 AU, which gives between 420 and
$470~\text{km}~\text{s}^{-1}$
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 Reference Leer and Holzer1980). Such complex developments will be undertaken in a later study.
We will focus then on the latitudinal distribution of the wind. At the 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 a 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
$20R_{\odot }$
, 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, Fujiki & Iju (Reference Tokumaru, Fujiki and Iju2015), 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 Reference Svalgaard and Kamide2013). So, like in the model, there seems to be a reduction of the north–south asymmetry between the surface and far away from the surface in the wind.
When we compare this result with the same figure from Pinto et al. (Reference Pinto, Brun, Jouve and Grappin2011), 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 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. (Reference Tokumaru, Kojima and Fujiki2010) and Sokół et al. (Reference Sokół, Swaczyna, Bzowski and Tokumaru2015) displays the same behaviour, with some kind of relaxation 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. (Reference Sokół, Swaczyna, Bzowski and Tokumaru2015).
Finally, we can notice that during the 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 due to the presence of high-latitude current sheets (Khabarova et al. Reference Khabarova, Malova, Kislov, Zelenyi, Obridko, Kharshiladze, Tokumaru, Sokół, Grzedzielski and Fujiki2017).
Figure 8 displays the Alfvén surface profile at the minimum and maximum (respectively a and b). The Alfvén surface corresponds to the surface at which the wind speed becomes faster that the Alfvén speed
$v_{A}=\sqrt{\Vert \boldsymbol{B}\Vert /(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C})}$
. 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 rather large, due to both the magnetic field still being stronger and the drop in the wind velocity observed in figure 7 between years 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.
Meridional cuts of the Alfvén surface (black line) at the minimum and maximum of activity (respectively a,b). The colour in the background shows the velocity projected onto the magnetic field in units 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.
The Alfvén radius
$r_{A}(\unicode[STIX]{x1D703})$
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. (Reference Pinto, Brun, Jouve and Grappin2011)) as the average cylindrical radius of the Alfvén surface weighted by the local mass flux
$r^{2}\sin \unicode[STIX]{x1D703}\unicode[STIX]{x1D70C}\boldsymbol{u}$
crossing the surface of a sphere:
$$\begin{eqnarray}\displaystyle \langle r_{A}\rangle =\frac{\displaystyle \int r^{2}\sin \unicode[STIX]{x1D703}\unicode[STIX]{x1D70C}u_{r}r_{A}(\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}}{\displaystyle \int r^{2}\sin \unicode[STIX]{x1D703}\Vert \unicode[STIX]{x1D70C}\boldsymbol{u}\Vert \,\text{d}\unicode[STIX]{x1D703}}. & & \displaystyle\end{eqnarray}$$
The mass loss rate is defined as:
$$\begin{eqnarray}\displaystyle {\dot{M}}=2\unicode[STIX]{x03C0}R_{0}^{2}\int _{0}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70C}u_{r}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}, & & \displaystyle\end{eqnarray}$$
where
$R_{0}$
is the radius at which the mass loss rate is measured (i.e. the radius of the spherical integration surface). In theory, the mass loss rate should be independent of
$R_{0}$
. In our case, the mass loss is evaluated at the outer boundary of the domain.
Time evolution of several global quantities of interest over an 11-year cycle: average Alfvén radius (a), mass loss (b) and angular momentum loss (c). 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.
In the same way, we define the angular momentum loss rate as:
$$\begin{eqnarray}\displaystyle \dot{J}=2\unicode[STIX]{x03C0}R_{0}^{3}\int _{0}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70C}\sin ^{2}\unicode[STIX]{x1D703}u_{r}\left(u_{\unicode[STIX]{x1D719}}-\frac{B_{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D70C}}\frac{\boldsymbol{B}_{p}\boldsymbol{\cdot }\boldsymbol{u}_{p}}{\Vert \boldsymbol{u}_{p}\Vert ^{2}}\right)\,\text{d}\unicode[STIX]{x1D703}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{B}_{p}$
and
$\boldsymbol{u}_{p}$
are respectively the poloidal magnetic field and velocity flow.
Figure 9 shows the evolution of these three quantities over a solar cycle, starting at the 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 (3.2), and the parity factor defined as follows:
$$\begin{eqnarray}\displaystyle P=\frac{\unicode[STIX]{x1D6FC}_{2,0}^{2}-\unicode[STIX]{x1D6FC}_{1,0}^{2}}{\unicode[STIX]{x1D6FC}_{2,0}^{2}+\unicode[STIX]{x1D6FC}_{1,0}^{2}}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FC}_{2,0}^{2}$
is the energy of the quadrupole and
$\unicode[STIX]{x1D6FC}_{1,0}^{2}$
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. (Reference Pizzo, Schwenn, Marsch, Rosenbauer, Muehlhaeuser and Neubauer1983) and Marsch & Richter (Reference Marsch and Richter1984), 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
$R_{\odot }$
: isothermal and polytropic models tend to give a lower estimate for the Alfvén radius (Pneuman & Kopp Reference Pneuman and Kopp1971). When we compare with Pinto et al. (Reference Pinto, Brun, Jouve and Grappin2011) 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
$R_{\odot }$
at the minimum to 2.2
$R_{\odot }$
at the 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 (Reference Réville and Brun2017) for the impact of two versus three dimensions and theoretical dynamo versus observational maps, we see that their Alfvén radius goes from 5
$R_{\odot }$
at minimum to 7
$R_{\odot }$
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. (Reference Pinto, Brun, Jouve and Grappin2011), we do have an anticorrelation, which is not perfect but does capture most of this effect, which could explain our intermediate results.
Mass flux in the meridional plane at different times in the cycle: (a) is at
$t=0$
years from the first minimum of activity, (b) at
$t=3.1$
years, (c) at
$t=4.9$
years and (d) at
$t=9.5$
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\times 10^{-14}~M_{\odot }/yr$
at maximum. From the data in McComas et al. (Reference McComas, Ebert, Elliott, Goldstein, Gosling, Schwadron and Skoug2008) and the calculations from Réville & Brun (Reference Réville and Brun2017), the mass loss inferred by Ulysses reached a minimum value of
$2.3\times 10^{-14}~M_{\odot }/yr$
in 1991 and a maximum value of
$3.1\times 10^{-14}~M_{\odot }/yr$
in 1992–1993. Our minimum mass loss is correct, but the amplitude of the variations is smaller. Wang (Reference Wang, Donahue and Bookbinder1998) 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 behaviour: it decreases until 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 behaviour 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 anticorrelated with the evolution of the magnetic surface energy, which is also the case in the model of Réville & Brun (Reference Réville and Brun2017). 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. (Reference Pinto, Brun, Jouve and Grappin2011), 20 % in Réville & Brun (Reference Réville and Brun2017)). 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 2-D 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. (Reference Pinto, Brun, Jouve and Grappin2011): the mass loss comes mainly from the coronal holes located at the poles near the equatorial streamer. However, since the reversals of the two hemispheres happen at different times but still close to one another, we can see that there are fewer coronal holes than expected than 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 with 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
$3.5\times 10^{29}~\text{g}~\text{cm}^{2}~\text{s}^{-1}$
at minimum of activity to
$1.0\times 10^{29}~\text{g}~\text{cm}^{2}~\text{s}^{-1}$
at maximum of activity, hence a variation of 70 %. This is the same trend as found in Pinto et al. (Reference Pinto, Brun, Jouve and Grappin2011) and Réville et al. (Reference Réville, Brun, Matt, Strugarek and Pinto2015a
): the star loses less angular momentum when it is more multipolar, as shown by the anticorrelation with the parity factor. It can be used to compute the magnetic spin-down time scale, defined as:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}t_{sd}=\frac{J_{\odot }}{\dot{J}}, & & \displaystyle\end{eqnarray}$$
where
$J_{\odot }$
is the Sun’s angular momentum, estimated at
$1.84\times 10^{48}~\text{g}~\text{cm}^{2}~\text{s}^{-1}$
in Pinto et al. (Reference Pinto, Brun, Jouve and Grappin2011) from a 1-D seismic solar model (Brun et al.
Reference Brun, Antia, Chitre and Zahn2002). This yields a magnetic spin-down time scale varying from
$1.7\times 10^{11}$
years to
$5.7\times 10^{11}$
years, which is of the same order of magnitude to that found in Pinto et al. (Reference Pinto, Brun, Jouve and Grappin2011).
Time evolution of the north–south asymmetry over an 11-year cycle for the average Alfvén radius (a), mass loss (b) and angular momentum loss (c). 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 (5.1), (5.2) and (5.3), except that the integral goes from 0 to
$\unicode[STIX]{x03C0}/2$
or from
$\unicode[STIX]{x03C0}/2$
to
$\unicode[STIX]{x03C0}$
. 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
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D70C}\boldsymbol{u}=0$
, 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 corresponds 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.
Reference Tokumaru, Fujiki and Iju2015).
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. (Reference Pinto, Brun, Jouve and Grappin2011), 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 the Babcock–Leighton flux-transport model described in Jouve & Brun (Reference Jouve and Brun2007) and including an asymmetry parameter
$\unicode[STIX]{x1D716}$
as in DeRosa et al. (Reference DeRosa, Brun and Hoeksema2012). 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, which 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. (Reference Pinto, Brun, Jouve and Grappin2011) the wind was considered isothermal) in spherical geometry (adapted from Réville et al. (Reference Réville, Brun, Matt, Strugarek and Pinto2015a
)). 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 of the wind.
We find that an anticorrelation 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 (Reference Réville and Brun2017). In Pinto et al. (Reference Pinto, Brun, Jouve and Grappin2011), 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 (Reference Réville and Brun2017) 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. Reference McComas, Ebert, Elliott, Goldstein, Gosling, Schwadron and Skoug2008). 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 changes 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 to that usually computed (see, e.g. Pinto et al. Reference Pinto, Brun, Jouve and Grappin2011). 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 anticorrelated with the parity factor, because they are highly sensitive to topology, 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. (Reference Sokół, Swaczyna, Bzowski and Tokumaru2015) using IPS data from Tokumaru et al. (Reference Tokumaru, Kojima and Fujiki2010), 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 had a 9 month delay for the dynamo generated field, and this 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 nonlinear coupling (Shukuya & Kusano Reference Shukuya and Kusano2017). 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 take it into account when inferring these 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 the direct consequence that all physical quantities may not necessarily exhibit a solar-like behaviour. For instance the polytropic approach is a simple approximation to 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 as 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 an interesting relation between dynamo and wind properties.
For the future, several aspects remain to be explored. A more realistic coronal heating, such as those prescribed in Cranmer et al. (Reference Cranmer, van Ballegooijen and Edgar2007) or Riley et al. (Reference Riley, Lionello, Linker, Cliver, Balogh, Beer, Charbonneau, Crooker, DeRosa and Lockwood2015), should be used to improve the realism of the wind solutions. Réville & Brun (Reference Réville and Brun2017) have shown that considering the full 3-D topology of the solar magnetic field also brings a lot more information and can help recover 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 3-D set-up of 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 time scales 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.
Acknowledgements
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.
