Dynamical coupling of a mean-field dynamo and its wind: Feedback loop over a stellar activity cycle
Barbara Perri, Allan Sacha Brun, Antoine Strugarek, Victor Réville
DDraft version February 3, 2021
Typeset using L A TEX default style in AASTeX62
Dynamical coupling of a mean-field dynamo and its wind: Feedback loop over a stellar activity cycle
Barbara Perri,
1, 2
Allan Sacha Brun, Antoine Strugarek, and Victor R´eville Universit´e Paris Saclay and Universit´e de Paris, CEA, CNRS,
AIM , F-91190 Gif-sur-Yvette, France Centre for mathematical Plasma Astrophysics, KU Leuven, Celestijnenlaan 200b-box 2400, 3001 Leuven, Belgium IRAP, Universit´e Toulouse III - Paul Sabatier, CNRS, CNES, Toulouse, France (Received December 31, 2020; Revised January 26, 2021; Accepted February 1, 2021)
Submitted to ApJABSTRACTWe focus on the connection between the internal dynamo magnetic field and the stellar wind. If thestar has a cyclic dynamo, the modulations of the magnetic field can affect the wind which in turn canback-react on the boundary conditions of the star, creating a feedback loop. We have developed a 2.5-dimensional numerical set-up to model this essential coupling. We have implemented an alpha-Omegamean-field dynamo in the PLUTO code and then coupled it to a spherical polytropic wind modelvia an interface composed of four grid layers with dedicated boundary conditions. We present here adynamo model close to a young Sun with cyclic magnetic activity. First we show how this model allowsto track the influence of the dynamo activity on the corona by displaying the correlation between theactivity cycle, the coronal structure and the time evolution of integrated quantities. Then we addthe feedback of the wind on the dynamo and discuss the changes observed in the dynamo symmetryand the wind variations. We explain these changes in terms of dynamo modes: in this parameterregime, the feedback loop leads to a coupling between the dynamo families via a preferred growth ofthe quadrupolar mode. We also study our interface in terms of magnetic helicity and show that itleads to a small injection in the dynamo. This model confirms the importance of coupling physicallyinternal and external stellar layers, as it has a direct impact on both the dynamo and the wind.
Keywords: dynamo, magnetohydrodynamics (MHD), Sun:corona, solar wind, solar-terrestrial relations INTRODUCTIONStars are complex objects that require to combine many fields of physics to be understood as a whole. In particularthe stellar interiors have traditionally been associated with fluid dynamics and later on magnetohydrodynamics (MHD),while stellar winds and atmospheres are usually described using plasma physics. To go beyond this first approach,these various communities had to gather their knowledge to produce the most consistent models possible. The bestexample for this is of course the closest star from us and hence the easiest one to observe in thorough details, the Sun.For more than four centuries the Sun has been observed displaying a cyclic activity through sunspot observations.Later on, Hale (1908) linked the sunspots to magnetic activity, which led to realize it has an 11-year period foramplitude and 22-year period for polarity. The general framework adopted nowadays to explain such a generation oflarge-scale magnetic field is the interface dynamo (Moffatt 1978; Parker 1993; Browning et al. 2006): the differentialrotation profile in the convection zone of the star (Schou & Bogart 1998; Thompson et al. 2003) leads to the generationof strong toroidal fields at the tachochline (Spiegel & Zahn 1992), which in turn is used to regenerate poloidal fieldsthanks to the combination of turbulence, buoyancy and Coriolis force at the surface. This dynamo loop allows for theamplification of the initial magnetic field until saturation, thus sustaining it against ohmic dissipation (Miesch 2005;
Corresponding author: Barbara [email protected] a r X i v : . [ a s t r o - ph . S R ] F e b Perri et al.
Brun & Browning 2017). To analyse the dynamo field, it can be decomposed into modes characterized by degrees (cid:96) and m by projecting it on the spherical harmonics base. Under the assumption of axisymmetry (e.g m = 0), odd (cid:96) are calleddipolar modes and constitutes the primary dynamo family, while even (cid:96) are called quadrupolar modes and constitutesthe secondary family (McFadden et al. 1991). A simplified yet efficient approach to model solar magnetism is themean-field dynamo (Roberts 1972; Krause & Raedler 1980), focusing on large-scale fields and assuming axisymmetry.The generation of toroidal field through differential rotation is then deemed the Ω effect and the regeneration of thepoloidal or toroidal field via turbulence is deemed the α effect. The combination of these two effects can give α -Ω, α or α -Ω dynamo loops. This description has the advantages of being easily implemented in MHD simulations withlow computational costs and yielding realistic results (Charbonneau 2020; Tobias 2019).On the other hand, it was not until the observations of the Mariner 2 spacecraft that the existence of a dynamic solaratmosphere being accelerated into a transsonic and transsalfv´enic wind was accepted by the community (Neugebauer& Snyder 1962). The first hydrodynamical description was given by Parker (1958) and magnetism was added bySchatzman (1962), Weber & Davis (1967), Mestel (1968) and Sakurai (1985) to yield a better description of thecorresponding torque applied to the star. The solar wind can be described using empirical models such as the WSAmodel (Wang & Sheeley 1990; Arge & Pizzo 2000) linking the terminal speed of a flux tube with its expansion factor.Another way is to use MHD simulations in 1D along a flux tube (Lionello et al. 2001; Suzuki et al. 2013; Pinto &Rouillard 2017), in 2D with axisymmetry (Keppens & Goedbloed 1999; Matt & Pudritz 2008; R´eville et al. 2015) orin 3D with a full description of the corona (T´oth et al. 2012; Usmanov et al. 2014; Riley et al. 2015; R´eville & Brun2017). The wind can even be better described using a multi-fluid approach (Usmanov et al. 2000; Hollweg & Isenberg2002) to take into account the influence of the different populations of particles. There are still open problems likethe precise mechanism behind the heating of the corona which can be approximated by a polytropic state with a lowadiabatic index, but is described better with detailed heating mechanisms such as magnetic reconnection (Parker 1988)or the turbulent dissipation of waves (Cranmer et al. 2007; R´eville et al. 2018; Shoda et al. 2020) in a non-trivial way.Recent observations, by the satellite
Ulysses for cycle 22 and the beginning of cycle 23 (McComas et al. 2003, 2008)or by the satellite
OMNI for cycles 23 and 24 (Owens et al. 2017), have shown that there is a correlation betweenthe 11-year dynamo cycle and the evolution of the corona. During a minimum of activity, the magnetic field is low inamplitude and its topology is mostly dipolar. The corona is also very structured with fast wind at the poles (around800 km/s) associated with coronal holes, and slow wind at the equator (around 400 km/s) associated with streamers(Wang et al. 2007). During a maximum of activity, the magnetic field is high in amplitude and its topology is a mixtureof high modes, dominated by the quadrupolar modes with equatorial symmetry (DeRosa et al. 2012), while fast andslow solar streams can be found at all co-latitudes. Complementary observations from Earth using IPS (Inter-PlanetaryScintillations, see Hewish et al. (1964) and Asai et al. (1998)) can help reconstruct the latitudinal evolution of thesolar wind speed over time. They display an 11-year pattern in the form of a cross due to the correlation between theactivity cycle and the latitudinal distribution of slow and fast streams (Tokumaru et al. 2010; Sok´o(cid:32)l et al. 2015). Thissuggests that there is a coupling operating between the interior and the exterior of the Sun, but the exact physicalmechanisms at play as well as the relevant timescales are still unclear (Wedemeyer-B¨ohm et al. 2009).From a theoretical and numerical point of view, it is however very difficult to study all of these layers simultaneously.The magnetic field and the wind evolve over a broad range of scales: the magnetic field evolves over 11 year withvarious scales ranging from several Megameters to a solar radius (Brun & Browning 2017), while the solar wind adaptsin a few hours, accelerating over a couple of tens of solar radii but extending up to several Astronomical Units (AUs)(Meyer-Vernet 2007). The physical properties of their respective environment are also very different: the interior ofthe star is characterized with convection (Rieutord & Rincon 2010) and large scale flows (Thompson et al. 2003; Basu& Antia 2010), while the exterior witnesses eruptive events in a dynamical atmosphere (Schrijver 2005). The rapidlychanging β plasma parameter (ratio of the thermal pressure over the magnetic pressure) at the surface of the star, orthe Mach number going from subsonic to supersonic in the corona are also evidence that the relative importance ofphysical phenomena can vary between the interior of the star and its extended atmosphere (Gary 2001). All of thesedisparities make the modeling of this coupling a numerical challenge.There have been many attempts to resolve this problem with various approaches. A first approach is to use aquasi-static approach, meaning that the coupling is realized through a series of wind relaxed states corresponding to asequence of magnetic field configurations evolving in time. These models can be data-driven, using series of magneticfield observations (Luhmann et al. 2002; Merkin et al. 2016; R´eville & Brun 2017), or rely on the numerical couplingbetween two codes dedicated respectively to the inner and outer layers of the Sun (Pinto et al. 2011; Perri et al. 2018). ynamical coupling of mean-field dynamo and wind NUMERICAL DYNAMO-WIND SET-UPWe begin by presenting the equations solved and our numerical set-up dedicated to the study of dynamo-windcoupling. We use the PLUTO code (Mignone et al. 2007), which is a compressible multi-physics and multi-solvercode with an extended MHD module. We first present the wind model in section 2.1, already described in Perri et al.(2018). Then we present the implementation of the mean-field dynamo and its validation in section 2.2. Finally wedescribe the numerical interface between the two computational zones in section 2.3. We will limit ourselves here toan axisymmetric (2.5D) description of the dynamo-wind coupling. The main consequence is that any magnetic fluxemerging through the stellar boundary is axisymmetric, and we cannot cope with magnetic structures localised inlongitude such as starspots and active regions. This limitation does not prevent us to assess how the coupling operateson a large-scale basis, and we leave the extension to fully 3D mean-field models (Karak & Miesch 2017; Kumar et al.2019) for future work. 2.1.
Wind model
Our wind model is adapted from R´eville et al. (2015). We solve the set of the conservative ideal MHD equationscomposed of the continuity equation for the density ρ , the momentum equation for the velocity field u with itsmomentum written m = ρ u , the energy equation with the energy noted E and the induction equation for the magneticfield B : ∂∂t ρ + ∇ · ρ u = 0 , (1) ∂∂t m + ∇ · ( mu − BB π + I p ) = ρ a , (2) ∂∂t E + ∇ · (( E + p ) u − B π ( u · B )) = m · a − π ∇ · [( η w ( ∇ × B ) × B )] , (3) ∂∂t B + ∇ · ( uB − Bu ) = −∇ × ( η w ∇ × B ) , (4)where p is the total pressure (thermal and magnetic), I is the identity matrix, a is a source term (gravitationalacceleration in our case), and η w is the resistivity in the wind domain. We use the ideal equation of state: ρε = p th / ( γ − , (5) Perri et al. 𝜌, 𝑃 →
Polytropic 𝑉 ’ = 𝑟𝑠𝑖𝑛𝜃Ω ∗ 𝑉 ∥ 𝐵 𝐵 → Background Field 𝜕 𝐵 ’ 𝜕𝑟 = 0𝑉 ∥ 𝐵 Ω 𝑅 ∗ DYNAMO WIND
Figure 1.
Presentation of the wind model. On the left is an example of the relaxed state of the corona for a dipolar magneticfield: the colorscale represents the poloidal speed in code units ( U = 437 km/ where p th is the thermal pressure, ε is the internal energy per mass and γ is the adiabatic exponent. This gives forthe energy: E = ρε + m / (2 ρ ) + B / (8 π ).PLUTO solves normalized equations, using three quantities to set all the others: length R , density ρ and speed U . If we note with ∗ the parameters related to the star and with 0 the parameters related to the normalization,we have R ∗ /R = 1, ρ ∗ /ρ = 1 and u kep /U = (cid:112) GM ∗ /R ∗ /U = 1, where u kep is the Keplerian speed at the stellarsurface, G the gravitational constant and M ∗ the stellar mass. By choosing the physical values of R , ρ and U , onecan deduce all of the other values given by the code in physical units. In our set-up, we choose R = R (cid:12) = 6 .
96 10 cm, ρ = ρ (cid:12) = 1 .
68 10 − g / cm (which corresponds to the density in the solar corona above 2.5 Mm, cf. Vernazzaet al. (1981)) and U = u kep, (cid:12) = 4 .
37 10 km/s. Our wind simulations are then controlled by three parameters: theadiabatic exponent γ for the polytropic wind, the rotation of the star normalized by the escape velocity u rot /u esc and the sound speed normalized also by the escape velocity c s /u esc . Note that the escape velocity is defined as u esc = √ u kep = (cid:112) GM ∗ /R ∗ . For the rotation speed, we take the solar value, which gives u rot /u esc = 2 .
93 10 − . Wechoose to set c s /u esc = 0 . . K hot corona for solar parameters and γ = 1 .
05. Thischoice of γ is dictated by the need to maintain an almost constant temperature as the wind expands, which is what isobserved in the solar wind. Hence choosing γ (cid:54) = 5 / u A /u esc which correspondsto the Alfv´en speed u A normalized by the escape velocity to control the amplitude of the initial magnetic field. Thevalue of this parameter depends on the dynamo parameters, as explained in Appendix C, and thus will be presentedin section 3.We assume axisymmetry but solve for the 3 components of flow and magnetic field. We use the spherical coordinates( r, θ ). 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 ( ρ, u , p, B ) to solve theRiemann problem corresponding to the previous set of equations. The time evolution is then implemented via a secondorder Runge-Kutta method. 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 thedeviations from a divergence-free field (Dedner et al. 2002).The numerical domain dedicated to the wind computation is an annular meridional cut with the co-latitude θ ∈ [0 , π ]and the radius r ∈ [1 . , R ∗ , as shown in the left panel of figure 1. We use an uniform grid in latitude with 256 pointsand a stretched grid in radius with 400 points, with the grid spacing geometrically increasing from ∆ r/R ∗ = 0 . ynamical coupling of mean-field dynamo and wind 𝜌, 𝑃 → Constant 𝑉 ’ = 𝑟𝑠𝑖𝑛𝜃Ω /011 𝑉 = 0𝐵 → Potential Field 𝐵 ’ = 0𝛼 effect and diffusivity 𝜂 Ω 𝑅 ∗ DYNAMO WIND
Figure 2.
Presentation of the dynamo model. On the left is an example of the dynamo evolution: the colorscale representsthe normalized toroidal magnetic field B φ in code units and the white lines are the contours of the poloidal magnetic field. Onthe right is a scheme of the top boundary conditions used for the dynamo computational domain alone. The wind domain is ingrey because it is not considered yet. at the surface of the star to ∆ r/R ∗ = 0 .
02 at the outer boundary. At the latitudinal boundaries ( θ = 0 and θ = π ) weset axisymmetric boundary conditions. At the top radial boundary ( r = 20 R ∗ ) we set an outflow boundary conditionwhich corresponds to ∂/∂r = 0 for all variables, except for the radial magnetic field where we enforce ∂ ( r B r ) /∂r = 0.Because the wind has opened the field lines and under the assumption of axisymmetry, this ensures the divergence-freeproperty of the magnetic field. The bottom boundary conditions that we would use with the wind computationaldomain only are shown in the right panel of figure 1. In the ghost cells (in green), the density ρ and pressure p areset to a polytropic profile, the rotation is uniform and the poloidal speed U pol is aligned with the poloidal magneticfield B pol . The latter is held fixed while the toroidal magnetic field B φ is linearly extrapolated from the numericaldomain. In the first point of the computational domain (in blue) all physical quantities are free to evolve, except forthe poloidal speed U pol which is forced to be aligned with the poloidal magnetic field B pol to minimize the generationof currents at the surface of the star and keep it as close as possible to a perfect conductor. We initialize the velocityfield with a spherically symmetric polytropic wind solution. In the grey area we do not solve yet the dynamo, whichwe will now describe. 2.2. Dynamo model
As we have seen in equation 4, PLUTO solves the full non-linear induction equation. For the dynamo inside the star,we will consider that the magnetic field B is the large-scale mean field and that the velocity u d is now only composedof a toroidal component u φ (we neglect the meridional flow for now) and thus implement an alternative form of theinduction equation: ∂∂t B + ∇ · ( u d B − Bu d ) = ∇ × ( α B ) − ∇ × ( η × ∇ × B ) , (6)where η is the effective magnetic diffusivity and α is a parametrized coefficient (the α -effect, that can be obtainedthrough the First Order Smoothing Approximation (FOSA) of the electro-motive force (Pouquet et al. 1976). Hencethe Ω effect is taken into account with the second term and the α effect with the third one. We can choose to have an α or α effect by disabling or enabling the third term for the equation on B φ . This form of the induction equation isonly active inside the star ( r < R ∗ ). The other equations 1, 2 and 3 are not solved inside the star.This new induction equation follows the same normalization as described before. However the mean field dynamocommunity usually refers to the control parameters C α and C Ω obtained with a different normalization. To make itmore convenient we will use in this article the traditional control parameters of the dynamo models when discussingour solution, but we have described in Appendix A the difference between the two and hence which conversion factorwe apply in our code to respect the PLUTO normalization.The flow in the dynamo domain will consist here only of a differential rotation such that u φ = Ω r sin θ Ω(r , θ ). Therotation in this zone is solar-like (Thompson et al. 2003) with a solid body rotation below 0 . R ∗ and differential Perri et al.
Case C critα ωA bench . ± .
002 158 . ± . A PLUTO A (cid:48) bench . ± .
002 157 . ± . A (cid:48) PLUTO B bench . ± .
003 172 . ± . B PLUTO B (cid:48) bench . ± .
002 168 . ± . B (cid:48) PLUTO
Table 1.
Comparison of the critical dynamo number C critα and the cycle period ω between the benchmark described in Jouveet al. (2008) and our dynamo set-up for various cases. The letters A and B refer to cases without and with a jump in diffusivity.Cases with a prime refer to radial boundary conditions, while cases without refer to potential boundary conditions. rotation above, with the equator rotating faster than the poles (25 versus 35 days). We use the following normalizedprofile: Ω( r, θ ) = Ω c + 12 (cid:18) (cid:18) r − r c d (cid:19)(cid:19) (cid:0) − Ω c − c cos θ (cid:1) , (7)where Ω c = 0 . r c = 0 . R ∗ , d = 0 . R ∗ and c = 0 .
2. The corresponding control number C Ω = Ω R ∗ /η t is alwaysfixed at 1 . × which corresponds to an equatorial rotation rate of Ω / π = 456 nHz for η t = 10 cm . s − , whichis in good agreement with observation of the actual solar rotation rate.As a first simple approximation the α effect is constant in the convection zone and zero in the radiative zone belowwith a smooth transition between the two zones. We use the following profile studied by the benchmark of Jouve et al.(2008): α ( r, θ ) = α √
34 sin θ cos θ (cid:18) (cid:18) r − r c d (cid:19)(cid:19) quench . (8)The physical amplitude of the α effect is given by the C α = α R (cid:63) /η t parameter, which value will vary. The quenchingterm T quench allows for saturation at the reference magnetic field value B quench : T quench = 1 + (cid:18) B φ ( r, θ, t ) B quench (cid:19) . (9)We have a swift transition in diffusivity between the radiative and the convection zone of two orders of magnitudes,following this profile: η ( r ) = η c + 12 ( η t − η c ) (cid:18) (cid:18) r − r c d (cid:19)(cid:19) , (10)with η c = 10 cm . s − and η t = 10 cm . s − .The magnetic field is initialized with a dipole confined in the convection zone: A φ = A R (cid:63) sin θr (cid:18) − (cid:18) . (cid:18) r c − rd (cid:19)(cid:19)(cid:19) , (11)where A = u A √ πρ ∗ R (cid:63) . Hence B φ is initially equal to 0.The numerical domain dedicated to the dynamo computation is an annular meridional cut with the co-latitude θ ∈ [0 , π ] and the radius r ∈ [0 . , . R ∗ , as shown in the left panel of figure 2. We use a uniform grid in latitudewith 256 points and a uniform grid in radius with 200 points, which yields a grid spacing of ∆ r/R ∗ = 0 . θ = 0 and θ = π ) we set axisymmetric boundary conditions. For the bottom boundarycondition ( r = 0 . R ∗ ) we use a perfect conductor condition: A φ = 0 , ∂ ( rB φ ) ∂r = 0 . (12)The top boundary conditions ( r = R ∗ ) that we use when only the dynamo model is running are shown in the rightpanel of figure 2. The poloidal magnetic field B pol is matched to a potential field following the method described inappendix B, while the toroidal magnetic field B φ is set to 0. ynamical coupling of mean-field dynamo and wind 𝜌, 𝑃 → Polytropic 𝑉 ’() = 0𝑉 , = 𝑟𝑠𝑖𝑛𝜃Ω 𝐵 ’() → Potential field 𝐵 , = 0 𝜌, 𝑃 → Polytropic 𝑉 ’() ∥ 𝐵 ’() 𝑉 , = 𝑟𝑠𝑖𝑛𝜃Ω 𝐵 ’() → Potential field 𝜕 𝐵 , 𝜕𝑟 = 0𝜌, 𝑃 → Constant 𝑉 ’() = 0𝑉 , = 𝑟𝑠𝑖𝑛𝜃Ω 𝑉 ’() ∥ 𝐵 ’() 𝑉 , = 𝑟𝑠𝑖𝑛𝜃Ω :;< 𝑅 ∗ + 4Δ𝑟 Ω 𝑅 ∗ WINDDYNAMOINTERFACE
Figure 3.
Description of the interface region. On the left is a scheme of the four layers between the dynamo and windcomputational domains. On the right are the descriptions of the various physical quantities in each corresponding region on thescheme. There both dynamo and wind domains are jointly considered, contrary to figures 2 and 1.
Since this kind of kinematic dynamo had never been implemented before in the PLUTO code, we validated ourmodel by comparing it with the benchmark described in Jouve et al. (2008). Table 1 compares the critical thresholdfor the dynamo number C critα to have a dynamo solution growing exponentially and the period ω of the correspondingcycle obtained in our dynamo set-up to the published benchmark cases. The letters A and B refer to cases withoutand with a jump in diffusivity. Cases with a prime refer to radial boundary conditions while cases without refer topotential boundary conditions. Our cases show a good agreement within the error bars found in the benchmark andwe are thus confident that our dynamo set-up is numerically robust.2.3. Wind-dynamo interface
The interface buffer between the dynamo and the wind computational domains is shown in figure 3. This interfaceis crucial to control and understand the interactions between the two domains. We will describe the chosen interfacein this section, combining the traditional dynamo and wind boundary conditions shown in figures 1 and 2. Now bothgrey zones are becoming active.The interface is divided into four layers of one grid point each because of the 2-point stencil used for our linearreconstruction. That way, the first two layers constitute the boundary condition for the dynamo and the last twolayers constitute the boundary condition for the wind. We also describe the first two points of the wind computationaldomain (blue in figure 3) where we add extra conditions for more stability.The first layer is very similar to the dynamo boundary conditions shown in the right panel of figure 2, except thatnow we prescribe the density and pressure to follow a polytropic law (but continuous with the constant value inside thestar). The magnetic field is extrapolated as a potential field using the value of the last point of the dynamo zone. Themathematical principle and numerical implementation of the potential field extrapolation are detailed in appendix B.In the last two layers of the interface, similar to the right panel of figure 1, the magnetic field is still extrapolatedbeyond R ∗ , the poloidal speed is aligned with the poloidal magnetic field and ∂ r B φ is set to 0 to limit the growth ofcurrents at the surface of the star.The second layer is special because we can enforce two types of boundary condition. It can either be similar to thefirst layer, so closer to the dynamo boundary condition: if we do so, the wind computational region has absolutely noimpact on the dynamo boundary conditions. This leads to what we call the one-way coupling where only the dynamocan influence the wind, without feedback from the latter. On the contrary we can choose to enforce a condition similarto the last two layers, so closer to the wind boundary condition: if we do so, the wind can back-react on the dynamoboundary conditions via the B φ condition which is derived from the wind computational domain. This is what we callthe two-way coupling where both the dynamo and the wind influence each other. Perri et al.
In the first two layers of the wind computational domain, we still impose the poloidal speed to limit again thegeneration of currents. We also impose some diffusivity a bit further than the star surface using the following expression: η w = 12 η (cid:18) (cid:18) r η − rd η (cid:19)(cid:19) , (13)with r η = 1 . R ∗ and d η = 0 . R ∗ . This allows the diffusivity to drop only after around ten grid points above theinterface. As we explain in section C, this additional diffusivity helps for the stability of our coupled dynamo-windmodel and allows the information from the dynamo domain to be smoothly transmitted to the wind domain.The final simulation box obtained is thus a combination of a dynamo computational domain, an interface bufferregion and a wind computational domain, as can be seen in figure 4. Thus the poloidal magnetic field generated bythe dynamo effect inside the star can impact directly the co-rotating idealized corona, with propagating timescalescharacterized in Appendix C. CONSISTENT EVOLUTION OF DYNAMO AND WIND QUANTITIES ALONG AN ACTIVITY CYCLETo validate the coupling from a theoretical point of view we focus on a model that develops a fast magnetic cycle,thus probably more related to a young Sun with the same rotation rate but enhanced activity. This allows us to reducethe discrepancy between the dynamo and the wind timescales and thus compute the evolution along several activitycycles using available computation resources. The main characteristics of this model, which we will now refer to asDW (for Dynamo-Wind), are given in table 2. We use the magnetic diffusivity to set the cycle period by adjusting thediffusive time t η = R ∗ /η t . Here it corresponds to a 5-day period for the magnetic cycle. Then we set C Ω to have thesolar rotation rate and finally we set C α to have the proper dynamo number D = C Ω × C α to trigger the dynamo. Thisyields a case where C Ω is set to a smaller value than C α (which is less frequent than the reverse), which correspondsto a strong generation of poloidal field due to turbulence. Because of this parameter range we initialize the magneticfield with a small value to prevent growth up to unrealistic values, with a value of u A /u esc of 0.05. Name η t (cm . s − ) C α C Ω B (G) u a /u esc γ c s /u esc ρ ∗ DW 3 . × . × . × . × − .
05 1 .
05 0 .
243 0 . Table 2.
Table of the main parameters for the test case DW.
Another important parameter to set is the quenching value of the toroidal field B quench . It is indeed crucial toensure that the large-scale magnetic field does not reach unrealistic values of u A /u esc (see R´eville et al. (2015)). Wehave thus derived an empirical criterion to stabilize the poloidal maximum surface energy close to the input value ofthe simulation. To do so, we have run 4 simulations with the same dynamo number D but with different ratios of C α and C Ω , and have found the following relationship: B quench = 0 . B surfpol (cid:114) C α C Ω . (14)We can then use our input parameters to compute an estimate of B quench , given a desired surface poloidal field.We will first present a one-way coupling version of the DW model, which we will call DW1. This allows us to showthe impact of the dynamo-generated magnetic field on the corona. For the first time we obtain a self-consistent solutionusing one simulation domain. This model is indeed already a great improvement in regards of the quasi-static methods(Perri et al. 2018) because it allows for a dynamical influence of the dynamo magnetic field on the wind without anyrisk of breaking causality (given that we respect the criteria derived in appendix C).We start by presenting the dynamo solution. The butterfly diagram obtained through this one-way coupling ispresented in the left panel of figure 4. We have chosen to display the evolution between approximately 1 .
07 and 1 . t η to avoid initial transients and focus on typical structures. The top panel shows the toroidal magnetic field B φ atthe base of the tachocline (0 . R ∗ ), the middle panel shows the radial magnetic field B r at the surface of the star andthe bottom panel shows the latitudinal magnetic field B θ also at the surface of the star. They are all in Gauss units.At the base of the tachocline (top panel), we can clearly see the reversal of polarity of the toroidal magnetic field,along the cyclic evolution of the surface poloidal field into as well (bottom panels). The period of the dynamo cycleis approximately 0.04 t η . The cycle displays both equatorial and polar branches at the surface due to our choice of ynamical coupling of mean-field dynamo and wind (a) (b) Figure 4.
Left panel: Butterfly diagram of the dynamo solution for case DW1. We show the toroidal magnetic field B φ at thebase of the convective zone ( r = 0 . R ∗ ) in the top panel, the radial magnetic field B r at the star surface in the middle panel,and the latitudinal magnetic field B θ at the surface in the bottom panel, all in Gauss units.Right panel: Meridional cuts of the dynamo-wind solution at different times for case DW1. We show the norm of the poloidalmagnetic field over the first 1.5 solar radii, in code units. Magnetic field lines are in white (full line when clockwise, dashedwhen counter clockwise). The surface of the star is indicated as a black dashed line. The corresponding times are shown asblack dashed lines on the butterfly diagram in the left panel. α effect. It corresponds to a so-called Parker-Yoshimura dynamo wave (Yoshimura 1975). Over the temporal rangeconsidered the cyclic magnetic field is always anti-symmetric with respect to the equator. The right panel of figure4 shows meridional cuts taken at moments of interest during the evolution of the solution. We can see the norm ofthe poloidal magnetic field over the first 1.5 solar radii, in code units ( B = 2 G). The surface of the star is markedby a black dashed line while the magnetic field lines are displayed in white (full line when clockwise, dashed whencounter clockwise). Snapshot (a) allows us to see the regeneration of the poloidal magnetic field due to the dynamoeffect inside the star, more precisely at the tachocline (0 . R ∗ ) due to our distribute α -effect. The magnetic structureappears mostly dipolar at large-scale, but more octupolar closer to the surface (both of anti-symmetric parity). We canalready see the field lines crossing the interface continuously into the corona. Snapshot (b) shows then the advection ofthe magnetic structures and their transmission to the corona thanks to our potential field boundary condition. Oncein the corona the structures interact with the stellar wind which tends to open the field lines and quickly decreasein amplitude. We begin to see the regeneration of the next cycle structures inside the tachocline, with the oppositepolarity. The corresponding times are shown as black dashed lines on the butterfly diagram in the left panel.Now we will present the response of the corona to this cyclic dynamo for the one-way coupling. The left panel offigure 5 shows meridional cuts taken at moments of interest during the evolution of the DW1 solution. We can seethe Mach number projected on the normalized magnetic field up to 20 R ∗ . This allows us to see at the same timethe various streams of the stellar wind and the corresponding polarity of the magnetic field associated. The blackline shows the Alfv´en surface which corresponds to the distance at which the wind speed becomes equal to the Alfv´enspeed. The surface of the star is marked by a black dashed line. We can see that the corona is highly dynamical withtransients generated at the surface of the star and carried by the wind to the external boundary of the computationaldomain. Due to the cyclic nature of the dynamo, the hemispheres reverse in polarity. In snapshot (c) the northernhemisphere is divided between two zones, one positive near the pole and one negative near the equator. In snapshot(d) however, we have a disposition similar but with reverse polarity: the coronal hole near the pole is now negativeand the equator is positive. This behavior is the same with opposite polarity in the southern hemisphere. Finally wecan notice that the Alfv´en surface is displaying various shapes along the cycle, which is a marker of the topologicalchanges of the dynamo field passing into the corona. Occasionally this surface can come very close to the star surface.This is due to both the low value of the magnetic field that we have to maintain with the quenching to ensure thecorrect coupling of the stellar interior and exterior, and also the dynamical evolution of the dynamo which can generatelocalized region of very low-amplitude surface field. We now turn to studying the long-term evolution of the corona0 Perri et al. (c) (d)
Figure 5.
Left panel: Meridional cuts of the dynamo-wind solution at two different times for case DW1. We show the quantity u · B / ( c S | B | ) which corresponds to the Mach number projected on the normalized magnetic field. The black line is the Alfv´ensurface which is the distance at which the wind speed becomes equal to the Alfv´en speed. The surface of the star is indicatedas a black dashed line. We clearly notice various polarity sectors.Right panel: Time-latitude diagram of the wind solution for case DW1. We show the wind speed u r at the far end of ourcomputational box ( r = 20 R ∗ ) in km/s. The corresponding times for the meridional cuts in the left panel are shown as blackdashed lines on the wind diagram. by considering the time-latitude diagram of the wind speed in the right panel of figure 5. The radial velocity u r istaken at the outer boundary of the domain (20 R ∗ ) and is shown in km/s. The corresponding times for the meridionalcuts in the left panel are shown as black dashed lines on the wind diagram. We observe latitudinal variations as aresponse to the coupling with the dynamo-generated magnetic field. They also present themselves as cyclic with aperiod which is half of the dynamo cycle (around 0 . t η vs 0 . t η ), because it is dictated by the topological changesrather than the polarity changes. We observe an alternation of slow (in blue) and fast (in red) streams at the polesand at the equator. Because of the rapidly changing magnetic field, we also observe regimes where fast streams are atall latitudes while the corona is adapting to the stellar interior. This is reminiscent of the solar wind observed duringthe minimum and maximum of the solar cycle (McComas et al. 2003): we don’t have the same structures because wedo not model the Sun but we still achieve a coupling between the dynamo field and the wind structures. The speedsreached are between 200 and 270 km/s at 20 R ∗ , which corresponds to between 400 and 520 km/s if we extrapolatethem at 1 AU. Thus we recover mostly the slow component of the wind, which is typical of polytropic winds and isexpected given our choice of c s /u esc .We now want to understand how the one-way coupling operates. Figure 6 shows time-radius diagrams of B r and B φ at two latitudes: the top panel is at co-latitude θ = 30 ◦ (closer to the northern pole), the bottom panel at co-latitude θ = 60 ◦ (mid-latitude in the northern hemisphere). This allows us to focus on the transmission of the magnetic fieldfrom the stellar interior to the corona, with the stellar surface marked by a black dashed line at 1 . R ∗ . The timescorresponding to the meridional cuts shown in figures 4 and 5 are shown as black dashed lines on the time-radiusdiagrams with the associated letter. We can clearly see the magnetic structures being generated at the base of theconvection zone at 0 . R ∗ , then being advected to the stellar surface and finally being transmitted to the corona withthe potential field boundary condition for B r on the left panel. The structures then are influenced by the solar windto the end of the computational domain but also quickly decrease in intensity. We can also see the reversals in polaritydue to the cyclic dynamo which translates into reversals of polarity in the corona as well, but with a delay due tothe transmission time. The various latitudes show that the continuity of the stellar corona is different depending onthe strength and geometry of the field. For example, θ = 60 ◦ corresponds to less intense dynamo structures whichdecrease rapidly at the surface of the star, so the transmission leads to less intense structures in the wind. For B φ we can see clearly the discontinuity created by our interface: the dynamo converges to a zero intensity toroidal fieldat the surface while the wind needs to have a stellar surface with small electric currents. This leads however to theperiodic generation of toroidal magnetic structures, which happen at the same time that the reversal in polarity for B r . It is thus very likely a result of the adjustment of the corona because of our condition on B φ and the necessity to ynamical coupling of mean-field dynamo and wind (a) (b) (c) (d) (a) (b) (c) (d) Figure 6.
Time-radius diagrams of B r and B φ (in code units, B = 2 G) for case DW1 at different latitudes: the top panelis at co-latitude θ = 30 ◦ (closer to the northern pole), the bottom panel at co-latitude θ = 60 ◦ (mid-latitude in the northernhemisphere). The stellar surface is marked by a dashed black line at 1 R ∗ . We only show the radial evolution from 0 . R ∗ to1 . R ∗ to focus on the most intense structures. The times corresponding to the meridional cuts shown in figures 4 and 5 areshown as black dashed lines on the time-radius diagrams with the associated letter. keep a divergence-free field in the domain. In this case this leads to the fact that the polarity of the structures of B φ in the corona is synchronized with the one from the interior. Here in the one-way coupling there is no feedback fromthe wind, so the dynamo itself is not different from a model without atmosphere.To complement our analysis of the DW1 case, we can further compute some key quantities of our stellar dynamomodel.We start with the magnetic characteristic length scale l b defined as: l b = (cid:80) (cid:96) (cid:96)β (cid:96), (cid:80) (cid:96) β (cid:96), , (15)where B r ( r = R ∗ , θ ) = (cid:80) (cid:96) max (cid:96) =1 β (cid:96), Y (cid:96) ( θ ), so that the coefficients β (cid:96), correspond to the decomposition of the surfaceradial magnetic field on the spherical harmonics. We only use the m = 0 modes because in our 2.5D simulationswe consider only axisymmetric fields. This length scale allows us to characterize the dominant scale in the magneticenergy power spectrum.We can also check the energy associated to a certain mode of the surface radial field, defined as follows: f (cid:96) = β (cid:96), (cid:80) (cid:96) max (cid:96) =1 β (cid:96), . (16)We will use these two quantities to characterize the minima and maxima of activity in our dynamo solution.The time evolutions of the magnetic characteristic length l b (black) and the fraction of energy in the dipolar com-ponent f (in red) are shown in figure 7. The characteristic length l b is evolving between 2.5 and 5, showing that thesurface magnetic energy of the star is best described by the first 5 modes. The maxima of l b will correspond to themaxima of the activity cycle because these are the moments where the surface magnetic energy is at its peak value.The dipole coefficient f evolves between 10 − and 50%, which means that the dynamo cycle oscillates between stateswhere the energy is mostly contained in the dipolar mode and states where the configuration is mostly multipolar justlike the Sun (DeRosa et al. 2012). The maxima of f will thus correspond to the minima of the activity cycle becausethese are the moments when the magnetic geometry is the most dipolar. We can observe an anti-correlation betweenthe two quantities just like what is observed for the Sun (DeRosa et al. 2012). However the minima of activity derivedfrom f are in phase-quadrature with the maxima of activity derived from l b , which is a characteristic of this dynamomodel.We can compare the evolution of integrated quantities to what has been found in quasi-static studies (Pinto et al.2011; R´eville & Brun 2017; Perri et al. 2018). We will focus on the average Alfv´en radius (cid:104) r A (cid:105) and the mass loss ˙ M .2 Perri et al.
Figure 7.
Time evolution of the magnetic characteristic length scale l b (in black) and the dipole coefficient f (in red) for caseDW1. Minima of activity are marked using red downward triangles, derived from f , and maxima with black upward triangles,derived from l b . Figure 8.
Time evolution of the integrated quantities of the wind for the case with one-way coupling. From left to right: theaverage Alfv´en radius (cid:104) r A (cid:105) in R ∗ and the mass loss ˙ M in M (cid:12) / yr. Minima of activity are marked using red downward triangles,derived from f , and maxima with black upward triangles, derived from l b . The Alfv´en radius is defined as the distance at which the wind speed equals the Alfv´en speed v A = B/ √ πρ . Theaverage Alfv´en radius, as defined in Pinto et al. (2011), corresponds to the Alfv´en radius averaged by the mass fluxthrough the surface of a sphere: (cid:104) r A (cid:105) = (cid:82) π r sinθρu r r A ( θ ) dθ (cid:82) π r sinθ || ρ u || dθ . (17)The mass loss is defined as follows: ˙ M = 2 πR s (cid:90) π ρu r sinθdθ, (18)where R s is the radius at which the mass loss is computed by integration through a spherical surface. In theory, whena stationnary state is reached the mass loss should be independent of R s , as long as it is big enough to include all theclosed magnetic structures. In our simulations, we perform an average of the mass losses computed with R s ∈ [9; 20] R ∗ to have a more precise value.The temporal evolution of these quantities is shown in figure 8. We have added the minima and maxima of activityas defined previously by l b and f with the same symbols as before: minima are marked using red downward triangles ynamical coupling of mean-field dynamo and wind (e)(f) Figure 9.
Left panel: Butterfly diagram of the dynamo solution for case DW2. We show the toroidal magnetic field B φ at thebase of the convective zone ( r = 0 . R ∗ ) in the top panel, the radial magnetic field B r at the star surface in the middle panel,and the latitudinal magnetic field B θ in the bottom panel, all in Gauss units.Right panel: Meridional cuts of the dynamo-wind solution at different times for case DW2. We show the norm of the poloidalmagnetic field over the first 1.5 solar radii, in code units. Magnetic field lines are in white (full line when positive, dashed whennegative). The surface of the star is indicated as a black dashed line. The corresponding times are shown as black dashed lineson the butterfly diagram in the left panel. and maxima with black upward triangles. The average Alfv´en radius varies between 1.5 and 3.0 R ∗ . This is smallerthan what was found in all the previous studies cited above, but this is mainly because we did not try to use solarparameters and because of the criterion C18 (see appendix C for more details) which limits us to weaker magneticfields. The mass loss varies between 4.5 and 6.4 10 − M (cid:12) / yr, which corresponds to a 42% variation. For comparison,the Sun has a variation of about 35% (McComas et al. 2008; R´eville & Brun 2017; Finley et al. 2019).In this study, we are more interested in the relative variations of these quantities than their values. We see mod-ulations in time with the evolution of the activity cycles, obtained for the first time in a completely self-consistentway. Based on similar studies, at minimum of activity we expect a bigger Alfv´en radius but a smaller mass loss. Atmaximum of activity, it is the opposite trend. The large-scale modulations are mostly in agreement with the cycles.At minimum of activity, when the field is mostly dipolar, we indeed see a higher Alfv´en radius and a lower mass loss.The maximum of activity is slightly less correlated with the integrated quantities, but most of the time it still doescorrespond to a smaller Alfv´en radius and a bigger mass loss. We also see many small-scale variations correspondingto fluctuations and transients, which explains why the correlation is not always perfect. This may be due to the factthat we have a rapidly changing magnetic field at the surface of the star which means the wind is constantly adjustingto the magnetic configuration. FEEDBACK-LOOP BETWEEN DYNAMO AND WINDWe will now present what happens when we enable the feedback of the wind, which means that the wind is now ableto back-react on the dynamo via the second layer of the interface and its condition on B φ (see figure 3). We will runthe same model as presented in table 2 but with full two-way coupling, so now we call this two-way coupling modelDW2. We will begin by showing the evolution of the global quantities as was done for DW1 in the last section. Thenwe will analyse the differences with the one-way coupling and quantify them in terms of dynamo modes and helicity.4.1. Impact of two-way coupling on global quantities
The corresponding butterfly diagram for the two-way coupling DW2 case is displayed in the left panel of figure 9with the same disposition as in figure 4. This new butterfly diagram is slightly different from the DW1 case. Thetoroidal magnetic field at the tachocline tends now to be more symmetrical with respect to the equator, starting at1 . t η , indicating that the quadrupolar modes (symmetric with respect to the equator) are dominating the dynamicsof the dynamo. At the surface the radial magnetic field is showing reversals of polarity but in an irregular way. It4 Perri et al.
Figure 10.
Time-latitude diagram of the wind solution for case DW2. We show the wind speed u r at the far end of ourcomputational box ( r = 20 R ∗ ) in km/s. We clearly see the alternance of slower and faster wind streams, and a huge degree ofasymmetry compared to figure 5. now presents north-south asymmetry most of the time and is alternating between symmetrical and anti-symmetricalconfigurations, which suggests an effective coupling of the dynamo families (DeRosa et al. 2012). The same effect ispresent as well for B θ . This asymmetry is highlighted by the snapshots in the right panel of figure 9. With the samepresentation as in figure 4, we can see in snapshot (e) a more symmetrical configuration when the poloidal field isgenerated at the tachocline. But when the field rises in snapshot (f), we see that the magnetic field lines crossingthe surface and influencing the corona have an asymmetric structure. The field exhibits more deformed closed loopsin the northern hemisphere because of the influence of the wind than in the southern hemisphere. The magneticstructures inside the star are also not symmetrical, with the northern hemisphere rising faster than the southern one.The corresponding times are shown as black dashed lines on the butterfly diagram in the left panel.In figure 10, we present the corresponding time-latitude diagram of the radial wind velocity u r for the two-waycoupling model. We do not display the corresponding meridional snapshots of the corona like in figure 5 because theevolution of the corona is very similar, apart from the north-south asymmetry we have already reported. The evolutionof the corona is still highly dynamical with transients. When the butterfly diagram is more symmetrical at the surface(for example between 1.15 and 1.25 t η ), we can see that the coronal structure becomes more irregular: there is a bigpatch of slow wind at the northern pole at 1.17 t η which is not reproduced elsewhere; fast streams seem also to bemore dominant in the southern hemisphere over this period of time considered. This shows already that the feedbackis introducing more variability in both the dynamo and the associated wind.Figure 11 shows time-radius diagrams of B r and B φ for case DW2 with two-way coupling. The times correspondingto the meridional cuts shown in figure 9 are shown as black dashed lines on the time-radius diagrams with the associatedletter. Once again, we clearly see the alternance of polarity and the efficient transmission of B r at θ = 30 ◦ . We canalso see that the feedback of the wind has modified the dynamo boundary conditions, so that at θ = 60 ◦ the magneticstructures do not decrease as much in intensity. For B φ , we are now left with only one point where B φ = 0, theother points depend on the dynamics of the wind. This results in more continuous and extended radial magnetic fieldstructures at the star surface, especially at θ = 60 ◦ . These structures are once again correlated with the changesof polarity of B r in the corona. In case DW2, this leads to an opposition of polarity between the structures in thestar and those in the corona. This is due to the specific period we have selected which focuses on a time with visiblenorth-south asymmetry. We have verified that at different times when the quadrupole modes are weaker, the polarityof B φ over the interface matches as in figure 6.Figure 12 shows the time evolution of l b and f for case DW2 with the two-way coupling. In this case, l b reachesslightly higher values, going up to 5.5. This is an indication that the cycle is slightly more multipolar than in caseDW1. We see that the two curves are still in anti-correlation over the considered time period, however the relativeposition of the minima and maxima is affected by the two-way coupling. We can still see a phase quadrature betweenthem at the beginning when the cycle is more anti-symmetrical like in case DW1. However when we start to be more ynamical coupling of mean-field dynamo and wind (e)(f) (e)(f) Figure 11.
Time-radius diagrams of B r and B φ (in code units, B = 2 G) for case DW2 at different latitudes: the top panelis at co-latitude θ = 30 ◦ (closer to the northern pole), the bottom panel at co-latitude θ = 60 ◦ (mid-latitude in the northernhemisphere). The stellar surface is marked by a dashed black line at 1 R ∗ . We only show the radial evolution from 0 . R ∗ and . . R ∗ to focus on the most intense structures. The times corresponding to the meridional cuts shown in figure 9 are shown asblack dashed lines on the time-radius diagrams with the associated letter. Figure 12.
Time evolution of the magnetic characteristic length scale l b (in black) and the dipole coefficient f (in red) for caseDW2. Minima of activity are marked using red downward triangles, derived from f , and maxima with black upward triangles,derived from l b . symmetrical (at 1.15 t η , when f reaches only 40% at maximum instead of 50%), the minima and maxima of activityderived begin to be in anti-correlation because l b is also affected by the two-way coupling and the feedback of the wind.We can then see in figure 13 the evolution of the integrated quantities defined above, namely the average Alfv´enradius and the mass loss as the various cycles proceed. The Alfv´en radius evolves between 1.5 and 2.75 R ∗ , which isa bit smaller as an interval with a smaller maximum value compared to case DW1. The mass loss evolves between4.7 and 6.4 10 − M (cid:12) / yr, which is again a bit smaller with a higher minimum value compared to case DW1. Thecorrelation with minima and maxima of activity is still working well when the cycle is anti-symmetrical, but whenit tends to be more symmetrical (between 1.15 and 1.25 t η ) the trends are not so clear anymore. This could showthat the trends derived in previous studies (described in section 3) are not generically true when the symmetry of thedynamo cycle changes, as f stops to be the most relevant marker of the dynamics of the cycle.6 Perri et al.
Figure 13.
Time evolution of the integrated quantities of the wind for case DW2. From left to right: the average Alfv´en radius (cid:104) r A (cid:105) in R ∗ and the mass loss ˙ M in M (cid:12) / yr. Minima of activity are marked using red downward triangles, derived from f , andmaxima with black upward triangles, derived from l b . Impact of the interface region on dynamo modes
As explained in the previous section, we see a clear difference between the two butterfly diagrams presented infigures 4 and 9. Without the wind influence, the tachocline toroidal and surface radial magnetic fields are equatoriallyanti-symmetric and the cycle is regular with a period of 0 . t η . With the influence of the wind, the cycle has the sameperiod but evolves to a more equatorially symmetric shape after 1.15 t η before returning to a more anti-symmetricpattern after 1.25 t η . It however never reaches full anti-symmetry, which suggests a coupling between the dynamofamilies (DeRosa et al. 2012).To understand this difference, we show in Figure 14 the evolution of the dipolar and quadrupolar modes ( (cid:96) = 1and (cid:96) = 2) for these two cases. The magnetic field is taken at the surface of the dynamo computational domain. Wealso show it for the entirety of the simulation, so between 0.0 and 1.75 t η , while previous figures only focused on thetime period between 1.0 and 1.3 t η . Without the influence of the wind (left panel), the dipolar mode has a ratherconstant maximal amplitude around 0.2 code units with cyclic variations with a period of approximately 0.03 t η (sameas the dynamo cycle), while the quadrupolar mode has an amplitude 4 orders of magnitude smaller. Hence, althoughit is continuously growing, the symmetric family remains negligible. But with the influence of the wind in case DW2,the quadrupolar mode grows to an amplitude equivalent of the dipolar mode (0.1 versus 0.2). This is however notconstant because we can clearly see a long-term modulation on the quadrupolar mode with periods of high intensity(between 1.15 and 1.25 t η as pointed before for example) and periods of very low values compared to the dipolar mode(between 1.0 and 1.15 t η for example). The feedback of the wind in this case has a visible influence by favoring thegrowth of the symmetric family by influencing the boundary conditions of the dynamo. This is a very interestingresult, demonstrating the need to couple both ways the dynamo and the wind.We nevertheless recall that the considered case has unusual values for a dynamo with C α >> C Ω . Tavakol et al.(1995) have shown that in such parameter regimes, the dynamo is highly non-linear and can switch easily fromsymmetric to anti-symmetric regimes because of the quenching or an asymmetry which results in the coupling of thedynamo families. We have also performed a threshold study similar to Jouve & Brun (2007) and have determined thatfor this set of parameters the quadrupolar mode has a growth threshold lower than the dipolar mode ( C Qα = 80 versus C Dα = 100). This case is therefore mostly a proof of concept to show that the feedback-loop between the dynamo andthe wind plays an important role in our simulations. More investigations are needed to assess whether this feedbackwould be as efficient in other parameter regimes or for other types of dynamos. Another possible explanation is thetreatment of magnetic helicity in our interface which we will investigate in the next section.4.3. Analysis of helicity in coupled dynamo-wind simulations
To quantify the exchange of information between the stellar interior and atmosphere, we now turn to magnetichelicity. Magnetic helicity is the measure of the warping and coiling of the magnetic field and is defined as: H = (cid:90) A · B dV, (19) ynamical coupling of mean-field dynamo and wind Figure 14.
Evolution of the dipolar and quadrupolar modes for model DW1 (without wind feedback) and DW2 (with windfeedback) (respectively on the left and on the right). This is displayed in code units. where B is the magnetic field and A is the vector potential associated to B such as B = ∇ × A . There is no unicityof the vector potential A , which means that magnetic helicity has to be computed with respect to a certain gauge. Toavoid this problem we can define the relative magnetic helicity which is gauge-invariant (Berger 1999): H rel = (cid:90) ( A + A p ) · ( B − B p ) dV, (20)where B p is a reference component, chosen to be potential ( ∇ × B p = 0) and which components match the componentsof B at the boundaries of the domain ( B p · n | S = B · n | S ). A p is the vector potential associated to B p such as B p = ∇ × A p .To compute B p and A p we use the fact that B p is chosen to be potential, which implies that it is deriving from ascalar magnetic potential φ ( B p = ∆ φ ). φ can be projected onto spherical harmonics and the boundary conditions thengive us enough constraints to determine it. To avoid discontinuity problems at the interface we compute a referencefield for the stellar interior and another for the corona, and we have checked after computation that the two match.To compute A , we combine direct integration of the magnetic field and a projection onto spherical harmonics to solveequations on the current j . For more details, please refer to appendix D.Figure 15 shows the time evolution of the relative magnetic helicity for case DW1 (on the left) and case DW2 (onthe right). We have selected a different time interval as what was previously shown, this time from 1 . . t η . Thisis to avoid the region where the dynamo modes are too coupled to really study only the influence of the interface onmagnetic helicity in the most basic situation. We have divided it into stellar interior (in blue/green colors) and exterior(in red/orange colors). We also decompose the helicity between its northern and southern components, respectively inblue and green for the interior and in orange and red for the exterior. This allows us to check the most basic propertyof magnetic helicity, which is its conservation in the numerical domain when magnetic diffusion is neglected. In thefigure, we can see that indeed the northern and southern components of helicity apparently compensate each otherrelatively well at each time for the star interior, and this for both coupling cases. This shows that the diffusivity insidethe star is not dominant enough to suppress this property on the time-scales considered here. We also can see that thehelicity is slightly more dominant in the exterior of the star with one order of magnitude of difference with the interior.However in the stellar corona, the conservation of relative helicity is visibly not as good as inside the star. We can seein the third panel, which shows the amplitude of the difference between northern and southern hemisphere relativehelicity for the stellar interior (in blue) and the exterior (in red) that the conservation of helicity in the exterior is notachieved in DW1 nor DW2. Ohmic diffusion does affect the relative helicity inside the star, but is unlikely to affect itin the corona where it weakly acts only in the very first points of the wind domain. This points to a possible loss orinjection of helicity at the interface and/or outer boundary. We will then analyze the fluxes of helicity to check thisassumption.Our set-up was designed to let only the poloidal field go through the dynamo-wind interface. Now we want to seehow our boundary condition affects magnetic helicity transfer into the idealized corona, as it has been advocated forinstance by Low (2013). To do so, we compute the evolution equation for H rel (Barnes 1988): ∂H rel ∂t = − πηc (cid:90) j · B dV + 2 (cid:90) ∂V [ A p × ( − u × B + η ∇ × B )] · n dS = − πηc (cid:90) j · B dV + 2 F H , (21)8 Perri et al.
Figure 15.
Time evolution of relative magnetic helicity for case DW1 (on the left) and DW2 (on the right), in code units. Thetop panel shows the helicity inside the star, and for the northern and southern hemisphere (in blue and green, respectively).The middle panel shows the same separation but for the helicity outside the star (northern hemisphere in orange, southernhemisphere in red). The last panel shows the difference between the northern and southern internal relative helicity in blue,and the difference between the northern and southern external relative helicity in red. where j is the current associated to the magnetic field B .Based on this equation, we can define a helicity flux through a surface associated to a temporal variation of helicity.Under the assumptions of axisymmetry, 2.5D, and spherical coordinates of this study, this yields: F H = 2 π (cid:90) π (( u φ B r − u r B φ ) − η ( ∇ × B ) θ ) A p,φ R s sin θ d θ, (22)where R s is the radius of the spherical surface through which the flux is computed. For more details about helicitybudget, please refer to appendix D.Figure 16 shows the time evolution of F H at the three different radii for case DW1 (on the left) and DW2 (on theright): in blue we have the top boundary condition for the wind at r = 20 R ∗ ; in green we have the bottom of the winddomain (just above the interface at r = R ∗ ); finally, in red, we have the top of the dynamo domain (just below theinterface at r = R ∗ ). The two upper panel show fluxes for the helicity budget in the wind domain: a positive (negative)value means that the wind helicity decreases. The bottom panel show the flux at the top of the dynamo domain: apositive (negative) value means that the dynamo helicity decreases. In all cases we show the decomposition of thehelicity fluxes on the northern (in dashed) and southern (in dots) hemisphere. As a sanity check, we have evaluatedthe flux loss at the bottom boundary condition of the dynamo, and found it to be completely negligible. At the topwind boundary condition, helicity is lost due to the modified outflow boundary condition for the wind. We find a lossof the order of 5 10 − (in code units) at maximum in both cases. With these numbers in mind, we can then see thatfor case DW1, the injection at the interface is negligible, being practically zero at the top of the dynamo and reachingonly 10 − at maximum at the bottom of the wind domain. However for case DW2, helicity transfer is enhanced at theinterface with typical amplitudes of 10 − at the bottom of the wind domain. We also notice that the northern andsouthern fluxes are different from one another in case DW2 due to the asymmetry discussed before (see section 4.2),while they compensate perfectly in case DW1. Note that in both cases we have verified that ohmic dissipation playsa negligible role, as shown in Appendix D.We therefore see that our interface plays a non-trivial role in helicity conservation and injection. The net helicityflux at the top of the dynamo domain is small in case DW1, but a substantial helicity is injected at the bottom of thewind domain due to our boundary condition on B φ . In case DW2, this effect is even more pronounced as our boundarydictates the helicity budget in the wind domain. It furthermore leads to helicity injection in the dynamo domain,which adapts strongly in case DW2 with the two-way coupling. This is probably another factor which explains whythe coupling of the modes is so strong in this case. This means that we need to find a more appropriate way to limitthe generation of currents at the bottom of our wind domain without jeopardizing the conservation of helicity. We ynamical coupling of mean-field dynamo and wind Figure 16.
Relative magnetic helicity fluxes for case DW1 (on the left) and DW2 (on the right) in code units. They are shownat three different radii: top of the wind computational domain ( r = 20 R ∗ ) in blue, bottom of the wind domain (just abovethe interface at r = R ∗ ) in green and top of the dynamo domain (just below the interface at r = R ∗ ) in red. The top andmiddle panels show the fluxes out of the wind domain (a positive/negative value means helicity is decreasing/increasing in thewind domain). The bottom panel shows the flux out of the dynamo top boundary (a positive/negative value means helicityis decreasing/increasing in the dynamo domain). Note how the magnetic helicity fluxes at the top dynamo boundary (bottompanel) differs between DW1 and DW2 models. In DW1 there is no net transfer, while there is a clear inward flux in DW2. plan to thus investigate further this condition in the interface for a more solar parameter space where the physicalquantities may be less sensitive to this effect. CONCLUSION AND DISCUSSIONWe have presented here one of the first example of self-consistent numerical coupling between an α -Ω stellar dynamoand a thermally driven stellar wind. We used the PLUTO code (Mignone et al. 2007) to couple the wind set-up ofR´eville et al. (2015) adapted for spherical coordinates in Perri et al. (2018), and the dynamo model described in thebenchmark of Jouve et al. (2008) and implemented for the first time in the PLUTO code. The key-point of this modelis the interface between the two physical domains: instead of being a software coupling between two independent codes,we have used a single numerical domain where a four-mesh thick layer structure at the surface of the star has beenimplemented. This allows us to solve self-consistently for the dynamo-wind coupling inside a computational domainencompassing the star’s interior and atmosphere from 0 . R ∗ up to 20 R ∗ . We can thus study dynamically the physicalexchanges between the dynamo and the wind along a dynamo cycle. The magnetic field generated inside the star isextrapolated using a potential field extrapolation across the interface layers to be transmitted to the corona and affectthe structure of the wind by its changes of amplitude and topology. We present here one of the first models able toquantify the feedback loop between the wind and the internal magnetic field.Here we lay out only the first developments of this model which must be seen as a proof of concept to demonstratethe efficiency and interest of this method. To do so, we presented a model with an activity cycle having a muchshorter period than the Sun. This model demonstrates cyclic variations of the corona along the dynamo cycle. Thesevariations are associated with the opening and closing of coronal holes with the corresponding generation of streamersand variations of the wind integrated quantities such as the average Alfv´en radius or the mass loss. These correlatewith the minima and maxima of activity of the cycles, reproducing results obtained so far only by data-driven models(Merkin et al. 2016; R´eville & Brun 2017) or quasi-static approaches (Pinto et al. 2011; Perri et al. 2018). Thisapproach however may generate transients from which it is difficult to evaluate if they carry physical information.By adjusting the interface layer, we are able to enable a one-way or two-way coupling between the dynamo and thewind, and thus we are able to illustrate a possible feedback loop between the dynamo and the wind. By enablingthe feedback of the wind, we observe a different dynamo behaviour where the wind feedback loop seems to enhancethe secondary family modes (Ossendrijver 2003; Svalgaard & Kamide 2013). We also demonstrate that our interfaceregion may generate additional helicity due to the various conditions on B φ , and are exploring new methods to modify0 Perri et al. this property. By focusing next on slower dynamo cycles we may have less transients and reconnection, which couldalso lead to more stable models.There are many different directions to broaden the accuracy and the applications of this model. To have a betterdescription of the physics, we have begun to implement a Babcock-Leighton mechanism with flux transport for thedynamo: this means that instead of the deep-seated α effect, we take into account the generation of strong toroidalmagnetic tubes at the base of the convective zone and their rise in this model by buoyancy, then twist by the Coriolisforce, and then allow the regeneration of poloidal field at the surface. The meridional circulation redistributes thepoloidal field at the base of the convective zone to be sheared by the differential rotation (Babcock 1961; Leighton 1969;Dikpati & Charbonneau 1999; Jouve & Brun 2007). This model needs the implementation of a non-local term whichis more complex given the architecture of the PLUTO code; we will report in a future work this new dynamo set-up.The description of the wind can also be improved by using a more realistic heating: instead of a polytropic equation,there are models based on turbulence and Alfv´en-wave heating that give results in good agreement with observations(Riley et al. 2015; R´eville et al. 2020; Hazra et al. 2021). For our next study, we will add these new elements andtackle the description of a coupled model between a solar-like dynamo with an 11-year cycle and a solar-like turbulentcorona, with an interface better designed to avoid injection of helicity. This is of course a numerical challenge giventhe discrepancy of the timescales with a timestep of the order of the hour to simulate at least 11 years of solar activity.On the other hand it will allow easier comparison with observational data to better calibrate the model (McComaset al. 2008). The parameter space description given in Appendix C will help us ensure that the coupling is operatingin an efficient and physical way. One of the more ambitious objectives would be also to describe more precisely thetransition region between the photosphere and the corona. In our model, the resolution is too coarse to allow aseparate description of the photosphere, chromosphere and corona, and the transitions in density and temperaturebetween these regions (Gary 2001; Meyer-Vernet 2007). Finally we remind the reader that this model could actuallybe very easily used for other stars of F, G or K spectral type with an internal structure similar of that of the Sun(with a convective zone at the surface and a radiative zone deep inside) and a thermal-driven wind. With elementssuch as the rotation profile, the internal poloidal flows or the dynamo mechanism prescribed with input parameters,we can easily change them to reproduce other stars. For example the two cases described in this article were chosenpurely theoretically but can be assimilated to young stars of T-Tauri type, with results similar to the work of vonRekowski & Brandenburg (2006). Acknowledgments
We thank Patrick Hennebelle and Thierry Foglizzo for useful suggestions. This work was supported by a CEA’Th`ese Phare’ grant, by CNRS and INSU/PNST program, by CNES Solar Orbiter and M´et´eo de l’espace funds andby the ERC Synergy grant WholeSun A. NORMALIZATION OF THE EQUATIONSDynamo and wind numerical codes usually have different normalisations for the MHD equations, leading to differentcontrol parameters for simulations. For the coupling, as dynamo and wind happen in the same numerical domain, wehad to choose one normalisation and convert the corresponding parameters. We chose the wind normalisation to focusonly on the induction equation in the PLUTO code. In this section, the two different normalisations are presentedalong with the conversion factors between the two for the induction equation.We start with the induction equation from the mean field theory (Moffatt 1978): ∂ B ∂t = ∇ × ( U pol × B + U φ × B − η ∇ × B + α B ) . (A1)Each quantity can then be written as a product between a constant physical amplitude with units and a normalizedprofile describing the dependency on ( r, θ, φ, t ): t = t ˜ t, r = L ˜ r, B = B (cid:101) B , U pol = U (cid:101) U , U φ = Ω L (cid:101) U φ , η = η ˜ η, α = α ˜ α. (A2) ynamical coupling of mean-field dynamo and wind ∂ (cid:101) B ∂ ˜ t = ∇ × (cid:18) t U L (cid:101) U pol × (cid:101) B + t Ω (cid:101) U φ × (cid:101) B − t η L ˜ η ∇ × (cid:101) B + α t L ˜ α (cid:101) B (cid:19) . (A3)To obtain the usual dynamo normalisation, we set the time and length constants to respectively the diffusive timeand the solar radius: t = R ∗ η t , L = R ∗ , (A4)with R ∗ = 6 .
96 10 cm and η t = 10 cm s − . By injecting A4 into equation A3, this leads to the following standardform of the induction equation: ∂ (cid:101) B ∂ ˜ t = ∇ × (cid:18) U R ∗ η t (cid:101) U pol × (cid:101) B + Ω R ∗ η t (cid:101) U φ × (cid:101) B − η η t ˜ η ∇ × (cid:101) B + α R ∗ η t ˜ α (cid:101) B (cid:19) = ∇ × (cid:18) R e (cid:101) U pol × (cid:101) B + C Ω (cid:101) U φ × (cid:101) B − η η t ˜ η ∇ × (cid:101) B + C α ˜ α (cid:101) B (cid:19) , (A5)with the standard control parameters being the magnetic Reynolds number R e , the normalized rotation rate C Ω andthe normalized alpha effect amplitude C α .For the wind normalisation used in the PLUTO code, we usually set the length, speed and time constants torespectively the solar radius, the Kepler speed at the surface of the star and the time ratio obtained by these twovalues: L = R ∗ , U = U K = (cid:114) GM ∗ R ∗ , t = L U , (A6)with G = 6 .
67 10 − dyne . cm . g − and M ∗ = M (cid:12) = 1 .
99 10 g. By injecting A6 into equation A3, we then obtain anew normalisation for the induction equation: ∂ (cid:101) B ∂ ˜ t = ∇ × (cid:18) U U K (cid:101) U pol × (cid:101) B + Ω R ∗ U K (cid:101) U φ × (cid:101) B − η U K R ∗ ˜ η ∇ × (cid:101) B + α U K ˜ α (cid:101) B (cid:19) = ∇ × (cid:16) R Pe (cid:101) U pol × (cid:101) B + C P Ω (cid:101) U φ × (cid:101) B − η P ˜ η ∇ × (cid:101) B + C Pα ˜ α (cid:101) B (cid:17) . (A7)To switch from dynamo to wind normalisation ( ie from equation A5 to equation A7), the following conversion factorcan be applied: ∂ (cid:101) B ∂ ˜ t = η t R ∗ U K ∇ × (cid:18) R e (cid:101) U pol × (cid:101) B + C Ω (cid:101) U φ × (cid:101) B − η η t ˜ η ∇ × (cid:101) B + C α ˜ α (cid:101) B (cid:19) . (A8) B. EXTRAPOLATION OF THE POLOIDAL MAGNETIC FIELDIn our interface region between the dynamo and the wind computational domains, we need to extrapolate themagnetic field generated by the dynamo inside the star to quantify the impact of its amplitude and topology on thewind. To do so, we have implemented a potential field extrapolation using the following mathematical method.We choose to perform a potential field extrapolation, meaning that the field at the surface of the star derives from ascalar potential Φ such as B = −∇ Φ which corresponds to a magnetic field in vacuum. Such a magnetic field has theproperty that both its divergence and curl are equal to 0, thus Φ is a solution to the Laplace equation which means itcan be decomposed as a sum of spherical harmonics Y m(cid:96) :Φ( r, θ, φ ) = (cid:96) max (cid:88) (cid:96) =0 (cid:96) (cid:88) m =0 (cid:16) Φ a(cid:96),m r (cid:96) + Φ b(cid:96),m r − ( (cid:96) +1) (cid:17) Y m(cid:96) ( θ, φ ) , (B9)where Y m(cid:96) ( θ, φ ) = c m(cid:96) P m(cid:96) ( θ ) e imφ , c m(cid:96) = (cid:113) (cid:96) +14 π ( (cid:96) − m )!( (cid:96) + m )! and P m(cid:96) ( θ, φ ) are the associated Legendre polynomials. As thefirst branch in r (cid:96) is divergent, we keep only the second one, choosing this form for the scalar potential:Φ( r, θ, φ ) = (cid:96) max (cid:88) (cid:96) =0 (cid:96) (cid:88) m =0 Φ b(cid:96),m r − ( (cid:96) +1) Y m(cid:96) ( θ, φ ) . (B10)2 Perri et al.
Applying the gradient operator to obtain the magnetic field and projecting on the vectorial spherical harmonics R m(cid:96) = Y m(cid:96) ( θ, φ ) e r and S m(cid:96) = ∇ ⊥ Y m(cid:96) ( θ, φ ) = ∂Y m(cid:96) ∂θ ( θ, φ ) e θ + θ ∂Y m(cid:96) ∂φ ( θ, φ ) e φ yields: B ( r, θ, φ ) = (cid:96) max (cid:88) (cid:96) =0 (cid:96) (cid:88) m =0 Φ b(cid:96),m ( (cid:96) + 1) r − ( (cid:96) +2) R m(cid:96) ( θ, φ ) − Φ b(cid:96),m r − ( (cid:96) +2) S m(cid:96) ( θ, φ ) , = + ∞ (cid:88) (cid:96) =0 (cid:96) (cid:88) m =0 α (cid:96),m ( r ) R m(cid:96) ( θ, φ ) + β (cid:96),m ( r ) (cid:96) + 1 S m(cid:96) ( θ, φ ) . (B11)We then assume that we know the value of the magnetic field at a certain radius named R b : B ( r = R b ) = + ∞ (cid:88) (cid:96) =0 (cid:96) (cid:88) m =0 α (cid:96),m ( R b ) R m(cid:96) + β (cid:96),m ( R b ) (cid:96) + 1 S m(cid:96) . (B12)This yields a system of two equations but with only one variable that is Φ b(cid:96),m . We have to choose to base ourselveson either B r or B θ to compute the solution. Based on previous benchmarks like Jouve et al. (2008) we realized thatmost dynamo models use B r as a reference, so that is what we chose as well. This yields then the solution:Φ b(cid:96),m = α (cid:96),m ( R b ) (cid:96) + 1 R (cid:96) +2 b . (B13)By injecting B13 into equation B11, this yields in the end: α (cid:96),m ( r > R b ) = α (cid:96),m ( R b ) (cid:18) rR b (cid:19) − ( (cid:96) +2) , (B14) β (cid:96),m ( r > R b ) = − α (cid:96),m ( R b ) (cid:18) rR b (cid:19) − ( (cid:96) +2) . (B15) C. PRELIMINARY STUDY WITH AN OSCILLATORY DIPOLEBefore diving into the full coupling with a dynamo cycle, we started with a preliminary study where we did not useyet the coupling model described above. The goal of this section is to have a more precise idea of the parameter spaceallowed for our model for the coupling to be operational. This will give a first indirect characterization of our couplingmodel by discussing the issue of time scales.In this simplified model, we have only the wind computational domain with the boundary conditions described infigure 1. The magnetic field is imposed at the surface of the star, but it varies in time: at each timestep, it oscillatesat a certain frequency ω which is a control parameter. The amplitude of the field is then given in code units by thefollowing relation: B osc = u A u esc (cid:112) πρ ∗ (1 + 0 . tω )) , (C16)where u A /u esc is the Alfv´en speed at the equator and at the surface of the star divided by the escape velocity, ρ ∗ is thedensity at the surface of the star. We recall that the Alfv´en speed u A is given by | B | / √ πρ . That way the magneticfield never reaches a null amplitude, it oscillates between 80% and 120% of its initial value. These oscillations areshown in the left panel of Figure 17. The wind then adapts to the magnetic field configuration: the right panel ofFigure 17 shows the radial speed u r of the wind at r = 14 . R ∗ (so at 75% of our numerical box) and at three differentlatitudes ( θ = 7 ◦ which is near the pole, θ = 45 ◦ which is at mid-latitude and θ = 84 ◦ which is near the equator). Wecan see oscillations with a frequency corresponding to the ω frequency of the magnetic field given in the left panel.Since we studied only a dipolar topology for the magnetic field, we have a latitudinal dependency, which means thatthe oscillations are more easily transmitted near the poles where the magnetic field influence is weaker so that thewind can adapt quicker, while the oscillations are less visible at the equator where the dipole is strong and thus createsa dead zone where the wind speed is reduced (between 0.45 and 0.48 vs 0.50 and 0.55 in PLUTO units).To quantify the propagation of information, we introduce three characteristic times: t A = ∆ ru A , t cyc = 2 πω , t η = (∆ r ) η , (C17) ynamical coupling of mean-field dynamo and wind Figure 17.
Time evolution in 1D of our oscillatory dipole. The left panel shows the evolution of B r which is imposed at thesurface of the star. The right panel shows the response in the wind speed u r at three different latitudes (near the poles, atmid-latitudes and near the equator). where t A is the local Alfv´en time, t cyc is the time associated with the oscillation of the surface magnetic field and t η is the magnetic diffusive time associated with η in the wind computational domain. Without diffusivity, our varioussimulations yield the following empirical criterion: t A < . t cyc . (C18)This is a first criterion concerning the resolution used for the numerical grid (cf. equation C17). For the testsperformed, we fixed the resolution with the following grid: we used an uniform grid of 256 points in θ between 0 and π , and a stretch grid of 256 points in r between 1 and 20 R ∗ with ∆ r = 0 .
01 at the surface of the star. We then usedvarious Alfv´en speeds to vary the corresponding Alfv´en time. Figure 18 shows how the wind is supposed to react ina case where the criterion C18 is respected and thus the coupling is operating properly. In the two meridional cuts,we can clearly see how the wind has adapted itself to the dipolar magnetic field: we have an equatorial streamer upto 4 R ∗ , elsewhere the magnetic field lines (in white) are open . The black line is the Alfv´en surface, which is thesurface at which u r = u A . This is the major visible change between the two configurations: on the left, the amplitudeof the magnetic field is minimal (80% of its initial value), thus the Alfv´en surface is further away from the star (onlatitudinal average, the corresponding Alfv´en radius r A is equal to 5.0 R (cid:63) ), while on the right the amplitude of themagnetic field is maximal (120% of its initial value) making the Alfv´en surface closer to the star ( r A = 3 . R (cid:63) ).If we choose to add magnetic diffusivity to the wind model, it can help propagate the information from inside the starto the outside. If the diffusivity is limited to a few grid points above the star surface (like with the profile described inequation 13), then it has little to no influence to the corresponding relaxed state of the wind. We still have to respectanother criterion to make sure that the wind has time to adapt to the magnetic cycle: t η < . t cyc . (C19)In the cases which did not respect criterion C18 and criterion C19, the coupling was not working properly. Whatwe observed was that the magnetic field was evolving too quickly at the star surface so that the wind could not adaptfast enough. This led to an uncoupling and thus a drop in the magnetic field, especially in the latitudinal componentof the magnetic field B θ . This leads to a strong back-reaction of the Lorentz force which is expressed as: F L = 14 π ∇ × B × B ⇒ πF L,θ = (cid:18) r ∂∂r ( rB θ ) − r ∂B r ∂θ (cid:19) B θ − r sin θ ∂∂θ (sin θ B φ ) B r . (C20) Please note that this is a misuse of language : the Maxwell laws prohibit the existence of a magnetic monopole, meaning that themagnetic field lines are not really ”open”, they simply reconnect far away from the star and way outside our numerical box, thus onlyappearing open in our simulations. Perri et al.
Figure 18.
Meridional cuts in 2D of the oscillatory dipole model. The left panel shows the case when the dipole amplitude isminimal (at 80% of its initial value), and the right panel the case when the dipole amplitude is maximal (at 120% of its initialvalue). The color bar represents the quantity u · B / ( c s || B || ), which is the wind speed projected on the magnetic field normalizedby the Mach number. The white lines represent the poloidal magnetic field lines. The black line represents the Alfv´en surfacewhich is the surface at which u r = u A . The star surface is represented by a black dashed line. This force triggers latitudinal flows, and in the end the wind speed is compressed to negative flows because of it, whichleads to accretion at the poles. D. COMPUTATIONS FOR HELICITYD.1.
Computation of A p with various boundary conditions D.1.1.
General principle
The relative magnetic helicity, which is gauge invariant, is defined as: H rel = (cid:90) ( A + A p ) · ( B − B p ) dV, (D21)where B p = ∇ × A p , but also where B p is potential ( ∇ × B p = 0) and is a reference component from B ( B p · n | S = B · n | S ).Thus that we want to find a potential magnetic field B p deriving from a scalar magnetic potential φ ( B p = −∇ φ )which we can project on spherical harmonics: φ ( r, θ, ϕ ) = (cid:88) (cid:96),m (cid:16) A m(cid:96) r (cid:96) + B m(cid:96) r − ( (cid:96) +1) (cid:17) Y m(cid:96) ( θ, ϕ ) . (D22)Because it is potential B p is also poloidal, thus we can introduce a quantity C p such as: B p = ∇ × ∇ × ( C p e r ) ⇒ A p = ∇ × ( C p e r ) . (D23) ynamical coupling of mean-field dynamo and wind B p and A p , we only need to compute C p , which can be linked to φ using the following relation: C p = − (cid:88) (cid:96),m (cid:18) A m(cid:96) (cid:96) + 1 r (cid:96) +1 − B m(cid:96) (cid:96) r − (cid:96) (cid:19) Y m(cid:96) . (D24)Under the assumption of axisymmetry, we have in the end: A p = − r ∂C p ∂θ e ϕ , B p = − ∂φ∂r e r − r ∂φ∂θ e θ , (D25) H rel = (cid:90) ( A · B − A · B p + A p B φ ) dV. (D26)The computation of A m(cid:96) and B m(cid:96) depends on the boundary conditions of the domain.D.1.2. Potential components with perfect conductor and potential field
For the dynamo, we have a perfect conductor boundary condition at the bottom, which leads to: ∂ r φ ( r = r b ) = 0 ⇔ (cid:96)r (cid:96) − b A m(cid:96) − ( (cid:96) + 1) r − ( (cid:96) +2) b B m(cid:96) = 0 . (D27)At the top we have a potential field boundary condition, which implies: ∂ r φ ( r = r t ) = − B r ( r = r t ) ⇔ (cid:96)r (cid:96) − t A m(cid:96) − ( (cid:96) + 1) r − ( (cid:96) +2) t B m(cid:96) = − G m(cid:96) , (D28)with B r ( r t ) = (cid:80) (cid:96),m G m(cid:96) Y m(cid:96) .We combine these two equations to obtain in the end: A m(cid:96) = G m(cid:96) (cid:96)r (cid:96) − t (cid:20)(cid:16) r b r t (cid:17) (cid:96) +1 − (cid:21) , B m(cid:96) = r (cid:96) +1 b G m(cid:96) ( (cid:96) + 1) r (cid:96) − t (cid:20)(cid:16) r b r t (cid:17) (cid:96) +1 − (cid:21) . (D29)D.1.3. Potential components with potential field and outflow boundaries
For the wind, we have a potential field boundary at the bottom, which leads to: ∂ r φ ( r = r b ) = − B r ( r = r b ) ⇔ (cid:96)r (cid:96) − b A m(cid:96) − ( (cid:96) + 1) r − ( (cid:96) +2) b B m(cid:96) = − G m(cid:96) , (D30)with B r ( r b ) = (cid:80) (cid:96),m G m(cid:96) Y m(cid:96) .At the top we have an outflow boundary condition, which is expressed using a modified outflow boundary conditionwith ∂ r ( B r r ) = 0 instead of ∂ r B r = 0, which leads to:2 r t B r ( r = r t ) + r t ∂ r B r ( r = r t ) = 0 ⇔ r t ∂ r φ ( r = r t ) + r t ∂ r φ ( r = r t ) = 0 . (D31)With expression D22, we obtain: r (cid:96)t A m(cid:96) + r − ( (cid:96) +1) t B m(cid:96) = 0 . (D32)We combine this condition with condition D30 to obtain in the end: B m(cid:96) = G m(cid:96) r (cid:96) +2 b (cid:16) r b r t (cid:17) (cid:96) +1 (cid:96) + ( (cid:96) + 1) , A m(cid:96) = − r − (2 (cid:96) +1) t B m(cid:96) . (D33)D.2. Computation of A We have a magnetic field B which derives from a vector potential A such that B = ∇ × A . Knowing B , we wantto recover A in 2.5D.To do so, we decompose the magnetic field B in a poloidal and a toroidal term in spherical coordinates: B = ∇ × ( C e ϕ ) + ∇ × ( A e r ) . (D34)6 Perri et al.
To obtain C , we can use direct integration, given that: Br = 1 r sin θ ∂∂θ (sin θ C) , B θ = − ∂∂ r (rC) . (D35)To obtain A , we use two different methods: direct integration for outside the star, projection on spherical harmonicsfor inside the star.With direct integration, we have: B ϕ = − r ∂ θ A. (D36)We can project A on spherical harmonics, which leads to: A = (cid:88) (cid:96),m A m(cid:96) Y m(cid:96) . (D37)Now if we take the r component of the projection of the current on spherical harmonics, we obtain: (cid:88) (cid:96),m (cid:96) ( (cid:96) + 1) r A m(cid:96) Y m(cid:96) = ∇ × B r . (D38)We can then link A and C to A using the following relation: A r = A, A θ = 0 , A ϕ = C. (D39)D.3. Balance for helicity fluxes
To make sure that all the helicity computations were carried out correctly, we have checked the balance describedby equation 21, which we rewrite here under the form: ∂H rel ∂t = Γ d + F H , (D40)where Γ d is the diffusion term due to the electric currents:Γ d = − πηc (cid:90) j · B dV, (D41)and F H is the helicity flux: F H = 2 (cid:90) ∂V [ A p × ( − u × B + η ∇ × B )] · n dS. (D42)The results are presented in figure 19. Case DW1 is on the left, case DW2 on the right. Top panel is the budget forthe dynamo domain, while bottom panel is the budget for the wind domain. For each domain, we plot in red ∂ t H rel , ingreen Γ d , in dashed blue the fluxes F H at the boundaries and in blue the sum Γ d + F H . We notice that in almost everycases, the diffusive term Γ d is negligible compared to the other terms. We include it however for more accuracy. Thebalance is qualitatively realized with the blue and red curves following each other in most instances. What we observeis that there are some minor errors, probably due to the numerical methods used to estimate the relative helicity andto derive the potential vector, and maybe to a lesser extend to small deviations from ∇ · B = 0 that are not correctedfast enough by our divergence cleaning method. The trends are nevertheless in good agreement, giving confidence inour helicity budget analysis. REFERENCES Arge, C. N., & Pizzo, V. J. 2000, Journal of GeophysicalResearch, 105, 10465. https://ui.adsabs.harvard.edu/abs/2000JGR...10510465A/abstract Asai, K., Kojima, M., Tokumaru, M., et al. 1998, Journalof Geophysical Research, 103, 1991. https://ui.adsabs.harvard.edu/abs/1998JGR...103.1991A/abstract ynamical coupling of mean-field dynamo and wind Figure 19.
Relative magnetic helicity budgets for case DW1 (on the left) and DW2 (on the right), in code units. Top panelis the budget for the dynamo domain while bottom panel is the budget for the wind domain. For each domain, we plot in red ∂ t H rel , in green Γ d , in dashed blue the fluxes F H at the boundaries and in blue the sum Γ d + F H . The balance is qualitativelyrealized with the blue and red curves following each other in most instances.Babcock, H. W. 1961, The Astrophysical Journal, 133, 572.https://ui.adsabs.harvard.edu/abs/1961ApJ...133..572B/abstractBarnes, D. C. 1988, The Physics of Fluids, 31, 2214,publisher: American Institute of Physics.http://aip.scitation.org/doi/10.1063/1.866622Basu, S., & Antia, H. M. 2010, The Astrophysical Journal,717, 488. https://ui.adsabs.harvard.edu/abs/2010ApJ...717..488B/abstractBerger, M. A. 1999, Plasma Phys. Control. Fusion, 41,B167. https://iopscience.iop.org/article/10.1088/0741-3335/41/12B/312Browning, M. K., Miesch, M. S., Brun, A. S., & Toomre, J.2006, ApJ, 648, L157.https://iopscience.iop.org/article/10.1086/507869Brun, A. S., & Browning, M. K. 2017, Living Reviews inSolar Physics, 14, 4. https://ui.adsabs.harvard.edu/abs/2017LRSP...14....4B/abstractCharbonneau, P. 2020, Living Rev Sol Phys, 17, 4.http://link.springer.com/10.1007/s41116-020-00025-6 Cranmer, S. R., van Ballegooijen, A. A., & Edgar, R. J.2007, The Astrophysical Journal Supplement Series, 171,520. https://ui.adsabs.harvard.edu/abs/2007ApJS..171..520C/abstractDedner, A., Kemm, F., Kr¨oner, D., et al. 2002, Journal ofComputational Physics, 175, 645. https://ui.adsabs.harvard.edu/abs/2002JCoPh.175..645D/abstractDeRosa, M. L., Brun, A. S., & Hoeksema, J. T. 2012, ApJ,757, 96, arXiv: 1208.1768.http://arxiv.org/abs/1208.1768Dikpati, M., & Charbonneau, P. 1999, The AstrophysicalJournal, 518, 508. https://ui.adsabs.harvard.edu/abs/1999ApJ...518..508D/abstractEinfeldt, B. 1988, SIAM Journal on Numerical Analysis, 25,294. https://ui.adsabs.harvard.edu/abs/1988SJNA...25..294E/abstractFinley, A. J., Deshmukh, S., Matt, S. P., Owens, M., & Wu,C.-J. 2019, The Astrophysical Journal, 883, 67.http://adsabs.harvard.edu/abs/2019ApJ...883...67FGary, G. A. 2001, Solar Physics, 203, 71. https://ui.adsabs.harvard.edu/abs/2001SoPh..203...71G/abstract Perri et al.