# Age Dependence of Wind Properties for Solar-type Stars: A 3D Study

Victor Réville, Colin P. Folsom, Antoine Strugarek, Allan Sacha Brun

► To cite this version:

Victor Réville, Colin P. Folsom, Antoine Strugarek, Allan Sacha Brun. Age Dependence of Wind Properties for Solar-type Stars: A 3D Study. The Astrophysical Journal, 2016, 832, 10.3847/0004-637X/832/2/145. insu-03669398

HAL Id: insu-03669398

https://insu.hal.science/insu-03669398v1

Submitted on 21 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

AGE DEPENDENCE OF WIND PROPERTIES FOR SOLAR-TYPE STARS: A 3D STUDY

VICTOR RÉVILLE<fnref n="1" />, COLIN P. FOLSOM<fnref n="2" />,<fnref n="3" />, ANTOINE STRUGAREK<fnref n="1" />,<fnref n="4" />, AND ALLAN SACHA BRUN<fnref n="1" />

ABSTRACT

Young and rapidly rotating stars are known for intense, dynamo-generated magnetic fields. Spectropolarimetric observations of those stars in precisely aged clusters are key input for gyrochronology and magnetochronology. We use Zeeman Doppler imaging maps of several young K-type stars of similar mass and radius but with various ages and rotational periods to perform three-dimensional (3D) numerical MHD simulations of their coronae and follow the evolution of their magnetic properties with age. Those simulations yield the coronal structure as well as the instant torque exerted by the magnetized, rotating wind on the star. As stars get older, we find that the angular momentum loss decreases with $\Omega_{*}^{3}$, which is the reason for the convergence on the Skumanich law. For the youngest stars of our sample, the angular momentum loss shows signs of saturation around $8\Omega_{*}$, which is a common value used in spin evolution models for K-type stars. We compare these results to semianalytical models and existing braking laws. We observe a complex wind-speed distribution for the youngest stars with slow, intermediate, and fast wind components, which are the result of interaction with intense and nonaxisymmetric magnetic fields. Consequently, in our simulations, the stellar wind structure in the equatorial plane of young stars varies significantly from a solar configuration, delivering insight about the past of the solar system interplanetary medium.

**Key words:** magnetohydrodynamics (MHD) – stars: winds, outflows – stars: magnetic field – stars: rotation – Sun: evolution – methods: numerical

1. INTRODUCTION

Among all the stellar properties, the characteristics of solar-like stars' winds are probably the most difficult to probe. Emissions throughout the electromagnetic spectrum unveil some of the properties of the photosphere and the coronae of stars, and internal structures can be probed with asteroseismology. Winds, on the contrary, produce very few detectable signals, although they are likely to exist in all stars possessing a hot corona, as supersonic outflows are the only stable way to balance coronal pressure with the near-zero interstellar medium pressure (Parker 1958; Velli 1994). Linsky & Wood (1996) have shown that absorption by neutral hydrogen at the astropause could be detected in Ly$\alpha$ spectra around astrophysical of nearby solar-like stars, unraveling properties of the stellar wind shocking against the interstellar medium. A growing sample of solar-type stars with positive detection for stellar winds led to a relationship between X-ray fluxes originating from coronal loops and mass loss rates (Wood et al. 2002). The “strength” of stellar winds, the mass loss rate $\dot{M}$, has consequently been related to the magnetic activity of the star. Wood et al. (2005) have obtained the relation $\dot{M} \propto F_{X}^{1.34\pm0.18}$, for $F_{X} \le 10^{6}$ erg cm$^{-2}$ s$^{-1}$, where $F_{X}$ is the X-ray flux. Beyond this value, weaker mass loss rates are observed, suggesting a saturation effect that is below the usual $F_{X}$ saturation value (Randich 2000; Pizzolato et al. 2003; Güdel 2004).

In parallel, the development of Zeeman Doppler imaging (ZDI) (Semel 1989; Donati & Brown 1997; Piskunov & Kochukhov 2002) has opened the study of surface magnetic fields for cool stars. Large-scale magnetic fields have been detected in the whole mass range that is thought to correspond to the existence of a convective envelope ($0.1 M_{\odot} - 1.4 M_{\odot}$).

The study of the geometrical and topological properties of the field in the context of stellar evolution is still in progress (Donati & Landstreet 2009; See et al. 2015) and raises theoretical questions about their generation through dynamo processes in convective envelopes (see Brun et al. 2015, and references therein). Nonetheless, the magnetic field amplitude of these stars has been shown to be a growing function of the rotation rate (Noyes et al. 1984; Brandenburg & Saar 2000; Vidotto et al. 2014a). This is necessary to explain the rotational braking of cool main-sequence dwarfs, as evolutionary models need the wind to carry angular momentum at a rate proportional to $\Omega_{*}^{3}$ (Kawaler 1988; Bouvier et al. 1997; Matt et al. 2015) all along the main sequence. However, recent studies suggest that the wind braking could stop or strongly decay for evolved stars, around a solar Rossby number $R_{\odot} \approx 2.5$ (van Saders et al. 2016), perhaps because of a change in magnetic topology (Garraffo et al. 2015a; Réville et al. 2015a). Hence wind, magnetism, and rotation are likely to evolve coherently through the whole life of solar-like stars.

After Schatzman (1962) understood that a magnetized outflow would carry away stellar angular momentum, Weber & Davis (1967) demonstrated that this loss is proportional to the Alfvén radius squared. Several studies followed to try to estimate the Alfvén radius from stellar parameters and thus give scaling laws for the angular momentum loss (Mestel 1968; Kawaler 1988). The latest braking laws have been successfully introduced in stellar evolution models (Matt et al. 2012; Gallet & Bouvier 2013), and we recently demonstrated that the magnetic topology could be included in those formulations through a simple scalar parameter, the open magnetic flux (Réville et al. 2015a, 2015b).

Most studies (see, e.g., Matt et al. 2012; Réville et al. 2015a) have been made in two dimensions with axisymmetric

<footnote n="1">¹Laboratoire AIM, DSM/IRFU/SAP, CEA Saclay, 91191 Gif-sur-Yvette Cedex, France; victor.reville@cea.fr, sacha.brun@cea.fr</footnote>

<footnote n="2">²IPAG, Université Joseph Fourier, B.P.53 F-38041 Grenoble Cedex 9, France; colin.folsom@obs.ujf-grenoble.fr</footnote>

<footnote n="3">³IRAP, CNRS—Université de Toulouse, 14 avenue Edouard Belin 31400 Toulouse, France</footnote>

<footnote n="4">⁴Département de physique, Université de Montréal, C.P. 6128 Succ. Centre-Ville, Montréal, QC H3C-3J7, Canada; strugarek@astro.umontreal.ca

*Received 2016 July 22; revised 2016 September 19; accepted 2016 September 19; published 2016 November 28*</footnote>
configurations (see Cohen & Drake 2014; Garraffo et al. 2015a, for a three-dimensional study of angular momentum loss with idealized magnetic field topologies) and were not able to capture the structure of complex magnetic fields observed by ZDI. Three-dimensional (3D) MHD simulations are now taking into account this complexity (Cohen et al. 2011; Vidotto et al. 2014b; do Nascimento et al. 2016; Alvarado-Gómez et al. 2016a, 2016b) to derive a self-consistent coronal structure. The complex structure of the corona is needed to study the interaction between stars and close-in planets, which has been shown to be very sensitive to 3D effects (Strugarek et al. 2015). Yet, to our knowledge, the influence of realistic magnetic fields on the long time variation of the wind properties has not been studied.

This work proposes to include observed, realistic magnetic
fields in the context of stellar evolution. We used spectro-
polarimetric observations of the surface fields of solar-like stars
to constrain 3D MHD simulations of stellar winds. The stars of
our sample share similar properties except their rotational
periods and their ages, which range from 25 Myr to 4.5 Gyr.
We developed a coherent framework to take into account the
evolution of the coronal properties with time, inspired by X-ray
flux observations, spin evolution models, and theoretical,
ab initio models (see Holzwarth & Jardine 2007; Cranmer &
Saar 2011; Suzuki et al. 2013). We confirm that the evolution
with age of global properties of the wind, such as the mass and
angular momentum loss, follows simple prescriptions in
agreement with the spin evolution models. These prescriptions
can be recovered by the semianalytical model we developed in
Réville et al. (2015b), except for the saturation of angular
momentum that appears only in our simulations. Also, the
three-dimensional structure of the young stars' winds shows
interesting features, such as a trimodal speed distribution, that
we explain through various interactions with the intense
magnetic field. We show that superradial expansion is a key
factor to explain the fastest wind components of a young stars'
wind. We also observe regions of fast wind encountering
slower streams in the equatorial plane, the so-called corotating
interaction regions (CIRs), that could be more common in the
wind of young stars.

This paper is organized as follows. The ingredients of our
model are described throughout Section 2. In Section 2.1, we
describe the observations used to constrain the surface
magnetic fields of our simulations. Section 2.2 describes our
prescriptions for the coronal properties and Section 2.3 our
numerical setup. The results are presented in two parts: Section
3, where we look at global properties such as the angular
momentum and mass loss over time, and Section 4, where the
tridimensional structure of the wind is detailed, with a special
focus on the velocity distribution. We summarize and reflect
upon our findings in Section 5.

2. MODEL INGREDIENTS AND DESCRIPTION

2.1. *Observational Data: ZDI Maps*

We consider in this paper six different stars whose ages are
precisely determined. Five of those stars belong to the study
made by Folsom et al. (2016) and share a mass and radius close
to (0.9 *M*⊙, 0.9 *R*⊙). Their rotational periods and surface
magnetic fields have been determined using observations from
the spectropolarimeters NARVAL (operating at Télescope
Bernard Lyot, France) and ESPaDOnS (operating at the

Canada-France-Hawaii Telescope, Hawaii). Their ages are
determined by studies of the open clusters and associations
they belong to. They were chosen to span a range of ages and
rotation rates but have similar masses, and bright stars were
selected to have sufficient signal-to-noise ratio.

The sixth star we include in our study is the Sun, for which we used a magnetogram obtained at the Wilcox observatory (DeRosa et al. 2012), and which will serve as a reference and oldest star in this paper. It is well known that the solar magnetic topology of the Sun and the solar wind structure vary over an 11-year cycle (see e.g., Pinto et al. 2011; Sokól et al. 2015). We considered the Sun in its minimum of activity in late 1996, during cycle 22. The ZDI maps of the five young stars are able to describe the surface magnetic fields of the stars up to a degree $l_{max} = 15$ in a spherical<fnref n="5" /> harmonics decomposition, which represents a sum of 135 different modes. The solar magnetograms made at the Wilcox observatory are able to reach much higher resolution ($l_{max} = 50$). We chose to cut the solar map to $l_{max} = 15$ to keep the same resolution for all stars in our sample. This partly justifies the choice of the solar minimum, whose energy spectrum is more concentrated in large-scale structure than during a maximum of activity. The first 15 spherical harmonics represent 95% of the magnetic energy during the minimum, and 80% during the maximum for cycle 22. Although the ages are sampled logarithmically, the sampling of the rotation rates is enough to follow closely the changes in coronal parameters and in magnetic field amplitude (see Section 2.2). The stellar parameters of all the targets are listed in Table 1.

Figure 1 shows the surface radial magnetic field reconstructed from the spherical harmonics decomposition of the Zeeman Doppler analysis (Folsom et al. 2016). The field is presented as orthographic projections on three different angles, with views facing the equator and the two poles. The color scales are chosen according to the amplitude of the magnetic field of the star. The average radial magnetic field of each case is given in Table 1 and mostly increases with the rotational frequency, as expected by dynamo theory (Durney & Latour 1978; Weiss 1994; Brun et al. 2015). We can see that, for each case, a dominant dipolar component is present, alongside smaller scale modes. This also motivated our choice of the solar minimum of activity, which exhibits a mostly dipolar field. Nonetheless, the dipolar components can show large inclinations with respect to the rotation axis, thus being far from an axisymmetric configuration and making a 3D approach necessary.

These magnetic maps are used as boundary conditions and specify the surface magnetic field of our simulations. However, they do not properly describe the stellar parameters by themselves. Like the magnetic field amplitude and topology, the thermodynamical properties of the base of the corona are likely to change with age and rotation (Güdel 2004; Giampapa 2005). The next section is dedicated to describing the model we used to take into account those variations.

2.2. Evolution of Coronal Properties with Age

Our numerical model needs, in addition to the surface
magnetic field, assumptions for the coronal base temperature
and density. Several studies have addressed the evolution of

<footnote n="5">⁵ The actual resolution of ZDI is somewhat lower than this ($l_{max}$ = 8–10), but the fields are reconstructed using 15 spherical harmonics.</footnote>
<table><caption>**Table 1**

Stellar Parameters of the Young Solar Analogs</caption><thead><tr><th>Name</th><th>Age (Myr)</th><th>Period (days)</th><th>Mass (M<sub>☉</sub>)</th><th>Radius (R<sub>☉</sub>)</th><th>T<sub>eff</sub> (K)</th><th>⟨B<sub>r</sub>⟩ (G)</th><th>% dipole</th><th>% axis.</th></tr></thead><tbody><tr><td>BD-16351</td><td>42 ± 6</td><td>3.3</td><td>0.9</td><td>0.9</td><td>5243</td><td>33.0</td><td>60</td><td>5</td></tr><tr><td>TYC 5164-567-1</td><td>120 ± 10</td><td>4.7</td><td>0.9</td><td>0.9</td><td>5130</td><td>48.8</td><td>78</td><td>78</td></tr><tr><td>HII 296</td><td>125 ± 8</td><td>2.6</td><td>0.9</td><td>0.9</td><td>5322</td><td>52.0</td><td>57</td><td>50</td></tr><tr><td>DX Leo</td><td>257 ± 47</td><td>5.4</td><td>0.9</td><td>0.9</td><td>5354</td><td>21.3</td><td>69</td><td>1</td></tr><tr><td>AV 2177</td><td>584 ± 10</td><td>8.4</td><td>0.9</td><td>0.9</td><td>5316</td><td>5.4</td><td>57</td><td>4</td></tr><tr><td>Sun 1996</td><td>4570</td><td>28</td><td>1.0</td><td>1.0</td><td>5778</td><td>1.1</td><td>35</td><td>75</td></tr></tbody></table>

Note. See Folsom et al. (2016) and references therein. The study of AV 2177 belongs to the second part of the aforementioned study.

those parameters with age and other stellar properties (Holzwarth & Jardine 2007; Cranmer & Saar 2011; Suzuki et al. 2013). Among others, their objective was to explain the mass loss rate signatures in the astropheres’ Lyα absorption spectra (see Wood 2004; Wood et al. 2005). Those studies are also nourished with a long history of X-ray coronal emissions (Pallavicini et al. 1981; Pallavicini 1989; Pizzolato et al. 2003; Güdel 2004; Wright et al. 2011), which show that coronal densities and temperatures tend to increase with the rotational period in solar-like stars. For instance, Holzwarth & Jardine (2007) gave scaling laws for the evolution of the temperature T and number density n as a function of the rotational frequency Ω* only assuming a power-law dependence:

$$T = T_{\odot} \left( \frac{\Omega_{*}}{\Omega_{\odot}} \right)^{0.1}, \quad n = n_{\odot} \left( \frac{\Omega_{*}}{\Omega_{\odot}} \right)^{0.6}. \qquad (1)$$

The values of the exponent for the power law correspond to their reference case, which aims to match the lowest branch of mass loss measurements. In their model, the mass loss is obtained computing one-dimensional (1D) polytropic and magnetocentrifugal wind from the coronal parameters (see Weber & Davis 1967; Sakurai 1985; Réville et al. 2015b). The study of Wood et al. (2005) shows, however, mass loss rates that can reach 100 times the solar value for rather slow rotators (see the case of 70 Oph, with a period of ≈20 days; Wood et al. 2005). Those extreme values need a stronger increase of the density with the rotation rate. Suzuki et al. (2013), using simulations of flux tubes heated by Alfvén wave dissipation, showed that such values could be reached and that a saturation could be obtained through the increase of the coronal density that increases the radiative losses. The dependence of the coronal density on Ω* is, however, also constrained by the observed X-ray fluxes and spin evolution models that suggest that 0.6 is a good estimate for the exponent (Ivanova & Taam 2003).

Moreover, the reference case fits the supposedly weak mass loss of more active stars (corresponding to most of the ages and rotation rates selected in our study) without invoking an additional transition below the $F_X$ saturation threshold, which still requires further theoretical understanding (Vidotto et al. 2016). We thus chose to keep the same prescriptions as Holzwarth & Jardine (2007), changing the solar reference temperature and density, as we use a different value for the polytropic index γ. In our model, $T_☉ = 1.5 \times 10^6$ K and $n_☉ = 10^8$ cm⁻³ are calibrated such that a wind with γ = 1.05 recovers a wind velocity of 444 km s⁻¹ at 1 AU and a mass loss rate of $3.2 \times 10^{-14} M_☉/yr^{-1}$. This value for γ has been widely used in the literature, including our works (Washimi &

Shibata 1993; Matt et al. 2012; Réville et al. 2015a), to describe the quasi-isothermal expansion of the wind with a polytropic model.

The temperature thus varies from $1.5 \times 10^6$ K to $1.9 \times 10^6$ K and the density from $1 \times 10^8$ g/cm⁻³ to $4.2 \times 10^8$ g/cm⁻³ throughout our sample (see Table 2).

## 2.3. Computational Methods and Boundary Conditions

In this study, we numerically solve the time-dependent ideal magnetohydrodynamics equations until a steady state is reached in our wind simulations. We use the PLUTO code (Mignone et al. 2007), using a finite-volume Godunov-type scheme and a Harten, Lax, van Leer, and Einfeldt (HLLE) solver (Einfeldt 1988) in three dimensions. Finite-volume methods aim to provide fully compressible and shock-capturing numerical methods that consider fluxes of conservative quantities through volumes. Hence, they formulate the MHD equations as a set of eight conservation equations defined as follows:

$$\frac{\partial}{\partial t} \rho + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (2)$$

$$\frac{\partial}{\partial t} \mathbf{m} + \nabla \cdot (\mathbf{m}\mathbf{v} - \mathbf{B}\mathbf{B} + \mathbf{I}\mathbf{p}) = -\rho \nabla \Phi + \rho \mathbf{a}, \qquad (3)$$

$$\frac{\partial}{\partial t} (E + \rho \Phi) + \nabla \cdot ((E + p + \rho \Phi)\mathbf{v} - \mathbf{B}(\mathbf{v} \cdot \mathbf{B})) = \mathbf{m} \cdot \mathbf{a}, \qquad (4)$$

$$\frac{\partial}{\partial t} \mathbf{B} + \nabla \cdot (\mathbf{v}\mathbf{B} - \mathbf{B}\mathbf{v}) = 0, \qquad (5)$$

where the energy $E \equiv \rho\epsilon + \rho v^2/2 + B^2/2$, the magnetic field **B**, the mass density ρ, and the momentum **m** ≡ ρ**v** are the conservative variables. Here, **v** is the velocity field, **p** = $p_{\text{th}} + B^2/2$ is the total (thermal plus magnetic) pressure, and **I** is the identity matrix. The potential Φ accounts for the gravitational attraction of the star, and **a** is a source term that contains the Coriolis and centrifugal forces, because we solve the equations in a rotating frame. The magnetic field is split into a background and a variable component for computational purposes (see Powell 1994). An ideal equation of state is used to close the set of MHD equations, and the internal energy is written

$$\epsilon = \frac{p}{\rho(\gamma - 1)}, \qquad (6)$$
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Reconstruction of the radial magnetic field for each star. The saturation value of the color scale changes with the amplitude of the radial magnetic field of each case. For all stars, we observe a dominant dipole field with smaller-scale features. The overall amplitude of the field decreases with the rotation rate of the star, and maps are sorted by age.</figcaption></figure>

with $\gamma = 1.05$, the ratio of specific heats, which differs from the usual value of 5/3 for a hydrogen gas in order to mimic the extended coronal heating.

We solve the equation in a Cartesian geometry with a grid centered on the star that extends from -30 to 30 stellar radii in each direction. The grid is uniform in a cube of $[-1.5 R_*, 1.5 R_*]^3$ with 192 grid points in each direction and then stretched up to $30 R_*$ with an additional 256 points for each direction. The resolution at the surface of the star is 50% larger than the one used in Réville et al. (2015a), but several tests have demonstrated that this resolution is enough to have numerical convergence.

The initialization is done by setting a spherically symmetric profile of a $\gamma = 1.05$ polytropic wind for the density, pressure, and poloidal velocity. This initial solution is obtained by a Newton-Raphson method on the normalized velocity and the critical radius of the polytropic wind solution. The MHD equations are then solved in a frame rotating with the star. We only initialize a solid body rotation inside the star so that it is the magnetic field that gives its rotating motion to the surrounding plasma. The magnetic field is initialized with a potential field extrapolation (Schatten et al. 1969) using the radial component of the ZDI map and a source surface radius $r_{ss} = 15 R_*$. This particular initial choice of the source surface radius has no impact on the final state since the extrapolated potential field then dynamically evolves with the stellar winds toward a relaxed state.

Boundary conditions at the surface of the star—which model the base of the corona in our<fnref n="6" /> case—are a sensitive point of our study. As in Réville et al. (2015a), we set three layers with different properties. For all layers, density and pressure are

maintained as the initial transonic polytropic wind solution. In the top layer, the poloidal velocity is set to be parallel to the magnetic field, while the toroidal velocity and the magnetic field are free to evolve. In the middle layer, the magnetic field is still free, the poloidal velocity is zero, and the toroidal velocity is set for solid body rotation. Finally, in the deepest layer, we enforce the reconstructed magnetic field, considering a perturbation in the toroidal field that self-adapts to minimize the overall current. More explicitly, the magnetic field solution interacting with the rotating wind in open regions will be in general different from the potential extrapolation we set at the initialization. Hence, to ensure a current-free magnetic field inside the star, which is supposed to be a perfect conductor, we dynamically modify $B_\phi$ in the deepest layer to get as close as possible to a curl-free magnetic field in the rotating frame (see Matt & Balick 2004; Zanni & Ferreira 2009). This boundary condition has a strong effect on the conservation of MHD invariants. For instance, the quantity that corresponds to the derivative of the electric field potential in axisymmetric configurations (see Ustyugova et al. 1999; Keppens & Goedbloed 2000; Matt et al. 2012; Réville et al. 2015a) remains in our 3D simulations close to the stellar rotation rate when our boundary condition is applied. We will see in Section 3.2 that this condition is key for an improved treatment of the angular momentum loss computation.

The outer edges of our domain are treated with *outflow* boundary conditions that set the derivative of each field normal to the boundary to zero.

### 3. RESULTS: GLOBAL PROPERTIES

#### 3.1. *Mass and Angular Momentum Loss*

Figure 2 shows the resulting steady-state wind solutions achieved in the six simulations. The convergence time is

<footnote n="6">⁶ See Matsumoto & Suzuki (2012) for a 2D model of the chromosphere and the transition region.</footnote>
<table><caption>**Table 2**
Input Parameters and Results of the Simulations</caption><thead><tr><th>Name</th><th>Ω<sub>*</sub>/Ω<sub>☉</sub></th><th>n (10<sup>8</sup> g cm<sup>-3</sup>)</th><th>T (10<sup>6</sup> K)</th><th>⟨R<sub>A</sub>/R<sub>*</sub>⟩</th><th>j̇/j̇<sub>☉</sub></th><th>Ṁ/Ṁ<sub>☉</sub></th><th>r<sub>ss, opt</sub></th><th>r<sub>ss, est</sub></th></tr></thead><tbody><tr><td>BD-16351</td><td>8.5</td><td>3.6</td><td>1.85</td><td>13.9</td><td>380</td><td>5.5</td><td>8.1</td><td>7.7</td></tr><tr><td>TYC 5164-567-1</td><td>5.9</td><td>2.9</td><td>1.8</td><td>16.7</td><td>190</td><td>2.75</td><td>10.7</td><td>10.5</td></tr><tr><td>HII 296</td><td>10.7</td><td>4.15</td><td>1.9</td><td>13.8</td><td>470</td><td>5.65</td><td>9.3</td><td>8.7</td></tr><tr><td>DX Leo</td><td>5.2</td><td>2.7</td><td>1.76</td><td>13.3</td><td>120</td><td>3.1</td><td>7.6</td><td>7.4</td></tr><tr><td>AV 2177</td><td>3.3</td><td>2.06</td><td>1.7</td><td>7.0</td><td>20</td><td>3.0</td><td>4.6</td><td>4.3</td></tr><tr><td>Sun</td><td>1.0</td><td>1.0</td><td>1.5</td><td>4.4</td><td>1.0</td><td>1.0</td><td>2.7</td><td>2.7</td></tr></tbody></table>

**Note.** Evolution of the input parameters and results of the simulation with age and rotational frequency. The effective Alfvén radii are coherent with the Alfvén surfaces shown in Figure 2, with comparable sizes for the four youngest stars and much smaller and decreasing values for AV 2177 and the Sun. The optimal and estimated $r_{ss}$ are very close to each other, with a slightly higher value for the optimal, although the change in resulting open flux is very small.

typically of the order of a few times the Alfvén timescale, that is, the time for Alfvén waves to cross the simulation’s domain. Given the high resolution, and the time step that can vary by one order of magnitude depending on the magnetic field amplitude, each of these simulations reaches a steady state after $10^5$ to $2 \times 10^6$ time steps (between $10^5$ to $5 \times 10^5$ core hours on supercomputers).

As usual, the topology of the coronal magnetic field in a steady state can be divided into open-field regions, or coronal holes, and closed-field regions, or dead zones, where the plasma corotates with the star. The surface magnetic field of all the stars in our sample includes a significant dipole, which can be strongly inclined. This mode gives its large-scale structure to the coronal magnetic field. However, the magnetic field inclination, amplitude, and smaller scale components lead to irregular shapes of the Alfvén surface, which is shown in Figure 2. For some cases, the Alfvén surface extends beyond the simulation domain. Indeed, for fast rotators, field collimation induces an increase of the poloidal magnetic field amplitude near the rotation axis, and the Alfvén surface is pushed farther away. A precise description of this phenomenon can be found in Washimi & Shibata (1993), Ferreira (2013), and Réville et al. (2015a), and we will address some of its consequences in Section 4. However, the global properties we are interested in, such as the mass and angular momentum loss rates, are constants within a small numerical variation once integrated over a surface that encloses the largest closed coronal loops. In the case of TYC 5164-561-1, which seems to have a significant part of its Alfvén surface out of the computation domain (more than any other case), this numerical variation of the mass and angular momentum loss is below 3%. It drops below 1% for cases where the Alfvén surface is fully inside the computation domain.

The mass loss rate $\dot{M}$ and the angular momentum loss rate $\dot{j}$ associated with the wind are computed as

$$ \dot{j} = \int_S \rho \Lambda v \cdot dS, \qquad (7) $$

where

$$ \Lambda = R \left( v_{\varphi} - B_{\varphi} \frac{B_p}{\rho v_p} \right), \qquad (8) $$

and

$$ \dot{M} = \int_S \rho v \cdot dS. \qquad (9) $$

The subscript $p$ and $\varphi$ stand for the poloidal and azimuthal components of each vectorial field, and $R$ is the cylindrical

radius. Those integrals can be computed from any surface $S$ that contains all of the closed coronal loops. From those outputs we define an effective Alfvén radius, which conveniently matches the relation given in the simplified model of Weber & Davis (1967):

$$ \langle R_A \rangle = \sqrt{\frac{j}{\Omega_* \dot{M}}}. \qquad (10) $$

This effective value matches quantitatively the average cylindrical radius on the irregularly shaped Alfvén surface of our simulations. All of the global quantities computed from our simulations are given in Table 2. We see that the angular momentum loss (AML) varies by a factor of 470 from the Sun to the fastest rotator HII 296. The variation of the mass loss is lower, with values that reach six times the solar value, here defined as $3.0 \times 10^{-14} M_{\odot} \text{yr}^{-1}$. The position of the Alfvén radius is globally increasing with rotation rate but is also very sensitive to temperature. For instance, the largest value we get is $16.6 R_{\star}$ for TYC 5164-567-1, which has a magnetic field similar to HII 296 but a slightly cooler coronal temperature, which makes the wind slower in our model.

One can recover the convergence on the Skumanich law for solar-like stars assuming a loss of angular momentum proportional to $\Omega_*^3$ in the nonsaturated regime. In the saturated regime, although there is no consensus,<fnref n="7" /> the dependence of the angular momentum loss is usually assumed to be linear with $\Omega_*$. Hence, a canonical expression for the stellar wind torque can be written as (Kawaler 1988; Bouvier et al. 1997)<fnref n="8" />

$$ \dot{j} = \dot{j}_{\odot} \left( \frac{\min(\Omega_{*}, \Omega_{\text{sat}})}{\Omega_{\odot}} \right)^2 \left( \frac{\Omega_{*}}{\Omega_{\odot}} \right) \left( \frac{M_{\odot}}{M_{*}} \frac{R_{*}}{R_{\odot}} \right)^{0.5} \qquad (11) $$

where $\Omega_{\text{sat}}$ is a saturation value of the angular momentum loss, which occurs at $\Omega_{\text{sat}} = 8\Omega_{\odot}$ for K-type stars. This saturation value corresponds to a Rossby number $Ro \approx 0.1$ (Wright et al. 2011), and higher mass stars have higher $\Omega_{\text{sat}}$ (for G stars the value is around $15\Omega_{\odot}$; see Gallet & Bouvier 2015). For clarity, we define $\tilde{\Omega}_{*} = \min(\Omega_{*}, \Omega_{\text{sat}})$, with $\Omega_{\text{sat}} = 8\Omega_{\text{sat}}$. Moreover, we have for every star in our sample $f(M_{*}, R_{*}) \equiv (M_{*}/M_{\odot})^{-0.5}(R_{*}/R_{\odot})^{0.5} \approx 1$, so formulation (11) reduces to $\dot{j}/\dot{j}_{\odot} = \tilde{\Omega}_{*}^2 \Omega_{*}/\Omega_{\odot}^3$.

In Figure 3, we compare the resulting torque from our simulations (in blue) to this empirical formulation (11) (in red).

<footnote n="7">⁷ For instance, a purely spherically symmetric radial Weber and Davis model yields a self-consistent saturation where $\dot{j} \propto \Omega_*^2$ (see Keppens et al. 1995).</footnote>

<footnote n="8">⁸ Recent work by Matt et al. (2015) proposed a more accurate description of the mass dependence.</footnote>
<figure><img src="image_5.png" /></figure>

BD-16351

<figure><img src="image_6.png" /></figure>

TYC 5164-567-1

<figure><img src="image_7.png" /></figure>

HII 296

<figure><img src="image_8.png" /></figure>

DX Leo

<figure><img src="image_9.png" /></figure>

AV 2177

<figure><img src="image_10.png" /></figure>

Sun 1996

**Figure 2.** Steady state of the 3D simulations. The Alfvén surfaces are shown in beige. The surface radial field is shown at the surface of the star, and field lines of positive polarity are shown in red, and negative polarity in blue. The Alfvén surfaces are largely shaped by the dominant dipole of each case. However, they are irregular, due to the dipole inclination and higher order components of the magnetic field. The size of the Alfvén surface grows with the amplitude of the field. We show a zoom-in on the stellar surface with smaller scale features at the surface in the top left corner of each panel.
All quantities have been normalized by the solar value. We look only at the $\Omega_*$-dependence of the angular momentum loss from the simulation, and the agreement is good. This agreement is the result of the combination of the observed magnetic fields, the hypothesis on the coronal temperature and density, and the simulation methods. Taking into account complex, three-dimensional magnetic fields is necessary, as some of our targets have a significant part of their magnetic energy in nonaxisymmetric modes. Moreover, we can see an inflection in the slope for the two fastest rotators BD-16351 and HII 296, which corresponds to the saturation value (which is not used whatsoever in our simulation). The saturation seems to appear self-consistently in our simulations, thanks to the plateau reached by the mass loss rates for our fast rotators (we will come back to this point in the next section).

The mass loss varies in our simulation (magenta line) approximately linearly with $\Omega_*$ (green line) and thus covers one order of magnitude over the sample. The variation of the angular momentum loss, which covers three orders of magnitude, is the result of three ingredients (see Equation (10)), the mass loss rate, the rotation rate, and the average Alfvén radius squared, and we can say that each of these ingredients accounts for one order of magnitude. Actually, the Alfvén radius squared seems to account for a bit more than the mass loss rate, but those are not independent parameters, and this statement must be handled with caution.

## 3.2. Semianalytical Model

Let us now compare these results with the semianalytical model we developed in Réville et al. (2015b) from the parameter study of Réville et al. (2015a). The Alfvén radii computed from the 3D simulations can be compared with the braking law we derived using a 2.5D axisymmetric setup in Réville et al. (2015a):

$$ \left\langle \frac{R_A}{R_*} \right\rangle = K_3 \left( \frac{\Upsilon_{\text{open}}}{(1 + (f/K_4)^2)^{0.5}} \right)^m, \quad (12) $$

where $K_3 = 0.65$, $K_4 = 0.06$, and $m = 0.3$ are fitted constants; $f = \Omega_* R_*/\sqrt{GM_*/R_*}$ is the breakup ratio; and $\Upsilon_{\text{open}} = \Phi_{\text{open}}^2/(R_*^2 \dot{M}_{\text{esc}})$ is a modified magnetization parameter (see ud-Doula & Owocki 2002). The open magnetic flux $\Phi_{\text{open}}$ is computed as the unsigned magnetic flux over a spherical surface beyond the largest closed magnetic loop. The value of the open flux should be constant whatever integration surface one chooses as long as it respects this latter condition (see Réville et al. 2015a, 2015b). Figure 4 shows how our simulations (green stars) fit this braking law. The blue line represents the braking law we derived in Réville et al. (2015a). Our set of simulations needs a small recalibration to be modeled by our braking law. Reducing the constant $K_3$ by 15% to a value of 0.55 gives an excellent agreement for all our cases.

The importance of the boundary conditions is illustrated in Figure 4. The two red stars are simulations for the Sun and TYC 5164-567-1 that were made by keeping fixed the initial extrapolated magnetic toroidal field inside the deepest layer of our boundary conditions. For the others (green stars), we set our self-adapting conditions ensuring a curl-free magnetic field inside the star. The deviation is around 25% on the Alfvén radius of TYC 5164-567-1 and the Sun, which induces a large

<figure><img src="image_11.png" /><figcaption>Figure 3. Evolution of the mass and angular momentum loss with age, compared with the rotational frequency and Equation (11). The angular momentum loss in blue follows the empirical law $\dot{J} \propto \dot{\Omega}_*^2 \Omega_*$. The variation of the mass loss is shown in magenta. The variation of the mass loss is close to linear with respect to $\Omega_*$ and is thus not enough to explain the total AML variation, which needs an increase of the Alfvén radius. The results of the semianalytical model are added with dashed lines.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 4. Comparison of the 3D cases and the braking law of Réville et al. (2015a). A new fit reducing by 15% the constant $K_3$ is needed for a good agreement with the data, due to higher coronal temperatures in our sample (green stars). Red stars in the figure stand for simulations of the Sun and TYC 5164-567-1 with simpler boundary conditions. Those Alfvén radii differ from the braking law by more than 25%, which is equivalent to an angular momentum loss divided by two.</figcaption></figure>

underestimation of the angular momentum loss that scales as $R_*^2$.

In Réville et al. (2015b), we proposed a method to compute the open flux based on an appropriate value for the source surface radius in a potential field extrapolation. We defined the optimal source surface as the source surface radius for which the potential field source surface model (Altschuler & Newkirk 1969; Schatten et al. 1969) gives the same open flux as the simulations. To estimate this optimal source surface radius, we considered a pressure balance between the thermal and ram pressure of a spherically symmetric magnetocentrifugal wind model with the magnetic pressure of the multipolar expansion of the surface magnetic field. We demonstrated the accuracy of this estimation, using 2.5D simulations, and we show in
Table 2 that it still holds for the 3D simulations performed here. The optimal ($r_{ss, opt}$) and estimated ($r_{ss, est}$) source surface radii are close even though $r_{ss, est}$ is systematically slightly smaller.

In terms of open flux, the results from the simulation and the potential extrapolation differ by less than 10% for all the cases. This good agreement is due to the right choice for the source surface radius. A potential extrapolation made with a constant $r_{ss}$ would have created a large discrepancy for part of the sample, given the large variation of optimal source surface radius. The large values of $r_{ss, opt}$ are consistent with the size of large coronal loops observed for the youngest stars of the sample with the largest magnetic fields that extend up to $10 R_*$ (see Figure 2). Moreover, $r_{ss, opt}$ values match the Alfvén points at the largest streamers' extremities and are always smaller than the effective Alfvén radius. This indicates there is less angular momentum at the cusp of the streamers than for the coronal holes (see Keppens & Goedbloed 2000; Garraffo et al. 2015b).

Hence, the semianalytical method described in Réville et al. (2015b) is likely to be successful for estimating the angular momentum loss from the mass loss and the open flux, if we adapt our formulation with the updated value of the $K_3$ constant. Going back to Figure 3, we superimposed the mass and angular momentum loss rates given by our semianalytical model with the dashed blue and magenta lines, respectively. The semianalytical model overestimates the torque compared to the simulations. This can be understood simply by looking at the mass loss of the spherically symmetric wind solution used in the semianalytical model and the mass loss of the simulation (dashed and solid magenta lines). In the simulations, the mass loss is always smaller than the spherically symmetric solution, since part of the plasma is contained in closed magnetic loops that cover a large part of the stellar surface. For strong magnetic fields, it can be two to three times less than the spherically symmetric ideal case. Hence the AML is consequently smaller. It is interesting to note that the saturation does not appear in the semianalytical model for the fastest rotators BD-16351 and HII 296. The saturation regime we observe in the simulations complies with a linear dependency of the AML with $\Omega_*$, and is thus different from the one that arises with a purely radial field in the Weber and Davis model (Keppens et al. 1995), which is used in the semianalytical model. Hence, in addition to the intrinsic saturation of the dynamo-generated magnetic field observed in X-rays, which should be contained in our magnetic maps, mass loss saturation due to confinement in large coronal loops could be involved in the saturation of angular momentum loss.

# 4. RESULTS: 3D STRUCTURE OF THE WINDS

## 4.1. Speed Distribution

The previous section showed that global properties of young stars' winds were coherent with simpler models parameterized on 3D MHD simulations. Integrated quantities such as the mass and angular momentum loss average out the complex and 3D structure of the wind. This structure is, however, relevant when it comes to studying the interaction of those winds with other objects such as companion stars or planets (see Cohen & Drake 2014; Vidotto et al. 2014b; Strugarek et al. 2015; do Nascimento et al. 2016).

The speed distribution of all the simulations is structured with slow and fast wind components, although no fast-wind-

related heating mechanism is included. Thermal heating is only provided by the coronal temperature and the equation of state through the choice of $\gamma = 1.05$, which is scaled with the slow solar wind component. Nevertheless, if the speed distribution is narrow in our solar case, the interaction with the strong magnetic field of our fast rotators yields broader distributions and can lead to interesting dynamical properties.

Figure 5 shows the histograms of the distribution of the wind speed projected on a sphere of radius $25 R_*$ for four stars of our sample. First, we observe a broadening of the speed distribution when the rotation rate of the star increases. While the speed distribution of the solar case is bracketed between 280 and $330 \text{ km s}^{-1}$, for the fast rotators TYC 5164 and HII 296, we observe a flattened distribution between 200 and $650 \text{ km s}^{-1}$. Since the temperature increases with rotation in our model, it is expected to reach higher speeds because of the thermal driving. The theoretical speed of the polytropic solution is indicated by the black dashed line and increases with rotation. We observe, however, both higher and lower components in the stellar winds of those fast rotators.

The distribution of fast rotators is organized with three peaks. In Figure 6, we show a volume rendering of the Mach number that corresponds to the four simulations of Figure 5. We can see that, for the young stars TYC 5164-567-1 and HII 296, the trimodal distribution appears with the blue (slow component), green (intermediate component), and red (fast component) colors. The first is due to the slow wind component emanating from the streamers. The strong dipolar magnetic fields of those stars create rings of slow winds located at the edge of the dead zones. Interestingly, those slow winds are positively correlated with the creation of currents in the simulation. The strong velocity gradient perpendicular to the dead zone boundary could be responsible for the creation of a current density even before the discontinuity in polarity that occurs beyond. This correlation needs, however, a more detailed analysis that is beyond the scope of this work.

The second peak corresponds to the theoretical speed at this distance from the star, given by the polytropic solution. Typically, this component will emanate from quiet polar regions, where the interaction with the magnetic field is the weakest. The third and fastest component appears between the streamers and slower winds blown at the poles. It seems to be confined in a flux tube at midlatitudes.

The additional acceleration given to the wind is likely to come from the magnetocentrifugal effect (Sakurai 1985; Washimi & Shibata 1993; Réville et al. 2015a, 2015b). Magnetocentrifugal acceleration is the consequence of the centrifugal force acting on the field lines embedded in the stellar rotation by the magnetic stress. Hence, the higher the rotation rate and the magnetic field, the higher is the magnetocentrifugal effect, which can be the dominant process in the wind acceleration (see Michel 1969; Washimi & Shibata 1993). For young stars, of periods of a few days, thermal driving and the magnetocentrifugal effect are comparable and must both be taken into account (Réville et al. 2015b). The magenta dashed line gives the wind speed obtained with the magnetocentrifugal wind solution (see Sakurai 1985; Réville et al. 2015b, for a detailed description of the solution calculation). This value seems to be coherent for the third peak of the speed distribution of the star HII 296. Open-field regions at midlatitudes are efficiently accelerated by this process (see Figure 6).
<figure><img src="image_13.png" /><figcaption>**Figure 5.** Speed distribution at $R = 25 R_*$ for four stars of our sample. A flattening and broadening of the distribution (all curves have a unitary integral) can be observed as the rotation and magnetic field increase. For fast rotators, three components can be identified: an average one around the polytropic solution at this radius (dashed black line), a slow component due to the slow wind emanating from streamers, and a fast component accelerated by the magnetocentrifugal effect and flux tube expansion. The speed of the magnetocentrifugal wind solution is shown in dashed magenta lines when different from the polytropic solution. The panels are sorted from left to right and top to bottom by rotation rate and not by age.</figcaption></figure>

In the case of TYC 5164-567-1, the theoretical Sakurai speed is, however, lower than the observed fast peak. This could be partly explained by the fact that the magnetocentrifugal wind solution is computed with the average surface field of the star. Hot magnetic spots conveniently located at the surface have a strong enough magnetic flux to explain the lowest part of the peak. However, another effect that is due to fast rotation and high magnetic field is also able to accelerate a stellar wind.

For a supersonic outflow, a superradial expansion of a given magnetic flux tube will necessarily accelerate the plasma (see Appendix). Here, two effects, both due to strong magnetic fields and fast rotation, can be accounted for by explaining the expansion of flux tubes. First, a flux tube located near the pole and yet close to the streamer boundary—a typical situation with an inclined dipole topology, with a maximum configuration precisely located in longitude—will have on one side its field lines driven by the streamers and thus bent downward, and on the other side, collimation of the field lines toward the rotation axis will bend them upward. This latitudinal expansion process is illustrated in Figure 7 in the left-hand panel.

A longitudinal expansion process can also occur when the fast flux tube originates at the boundary of the strong concentration of flux near the pole. The magnetic stress on each side of the flux tube is different and leads to a differential efficiency of the magnetocentrifugal effect on each side. Typically, the frozen-in magnetic field line that originates inside the “hot spot” will rotate faster than a magnetic field line originating outside. This effect is consequently responsible for a longitudinal expansion of the flux tube and is illustrated in the right-hand panel of Figure 7.

For the fast rotators, the acceleration of fast streams clearly starts just beyond the sonic surface. We have estimated the superradial expansion factor defined by

$$f_{\text{exp}} = \left(\frac{r_c}{r}\right)^2 \frac{A}{A_c}, \qquad (13)$$

where $r$ is the spherical radius, and $A$ is the surface of the section of the magnetic flux tube, computed at $r = 18 R_*$ and at the sonic critical point (subscript $c$). In the case of TYC 5164-567-1, we find $f_{\text{exp}}$ to increase between 6.5 and 10 from the
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 6.** Volume rendering of the Mach number for the stars of Figure 5. The three peaks in the velocity distribution of Figure 5 are visible in blue, green, and red, for the slow, intermediate, and fast wind components, respectively. The color scale is linear and the Mach number varies between 0 and 4.</figcaption></figure>

outer boundary to the core of the fastest flux tube. The maximum speed observed in the core of the flux tube is around $1300 \text{ km s}^{-1}$, meaning around three times the polytropic wind speed at this distance from the star.

The broadening of the speed distribution, by the interaction of a thermally heated corona with the strong magnetic field and the fast rotation, seems to be a reliable feature in our simulations. The trimodal distribution of speeds, using polytropic and magnetocentrifugal models, could be a simple input in other models, for instance the one used to compute mass loss rates (see Wood 2004, where they used constant

stellar wind velocities in their multispecies simulations of the termination shock of astrospheres).

It is hard, though, to predict the effect of additional accelerating mechanisms in coronal holes—which are necessary to explain the structure of the solar wind—in the case of young stars. We can imagine that the trimodal distribution would remain, but that the separation between the slowest mode (which comes from the streamers) and the two others (which originate from coronal holes) would be widened. Also, the two fastest peaks would show higher speeds than the one obtained in our simulations.
<figure><img src="image_18.png" /><figcaption>**Figure 7.** Sketch of the latitudinal and longitudinal superradial expansion processes. In the left-hand panel, the collimation of field lines toward the axis due to high rotation is responsible for the superradial expansion. On the right, the magnetic flux gradient at the stellar surface will generate a differential magnetic stress on both ends of the superradially expanded flux tube. See also Figure 6.</figcaption></figure>

## 4.2. Slow and Fast Winds in the Equatorial Plane

The magnetic topology of the Sun goes from a strong equatorial dipole at minimum to a more quadrupolar configuration during maximum (DeRosa et al. 2012). The heliospheric current sheet (HCS) is corrugated due to the rotationally modulated direction of the dipolar moment, and north and south sector polarities are observed at 1 AU depending on when the Earth is beneath or above it. Hence, streams of fast wind encounter slower and higher density wind that wraps the HCS. Indeed, the angle of the stream $\psi = \arctan(r\Omega_*/v)$ (Richardson 2004) is a decreasing function of the velocity of the flow. Streamlines of the fast component are thus less bent and compress against slower and denser wind streams. These so-called corotating interaction regions (Belcher & Davis 1971) were detected in the solar system thanks to Mariner 5. In the solar system, the spacecraft Pioneer 10 and 11 have shown that 75% of the CIRs have formed shocks between 3 and 5 AU (Smith & Wolfe 1976). Those shocks dissipate the energy of the solar wind and are able, for instance, to accelerate ions (Richardson 2004).

In the case of rapidly rotating stars, whose dipoles are significantly inclined (BD-16351 and DX Leo in our study), we should find similar features enhanced by rapid rotation and strong magnetic fields. Figure 8 shows the Mach number in the equatorial plane of BD-16351 and the Sun in 1996. The flow is structured as a Parker spiral in both cases (Parker 1958). Streamers can be seen in the equatorial plane of BD-16351 because of inclination of its dipolar magnetic field. The slow wind can be clearly identified following the streamers in that case. The flow is drastically more uniform in the case of the (quiet) Sun, and the Parker spiral shows weaker inclination of the field lines.

In the case of BD-16351, we observe adjacent fast and slow streams that encounter each other at the edge of the domain. The direction of the field (and thus the flow) seems different for the fast stream than for the slow stream, and a compression region occurs. Figure 9 shows the profiles of the density, the Mach number, the magnetic field, and the current density amplitude along a 25 $R_*$ radius circle in the equatorial plane of BD-16351. Two density peaks mark the slow wind components. These peaks correspond to troughs in the profiles of the Mach number and the magnetic field amplitude. The thin throat in the magnetic field profile matches the crossing of the

current sheet that also corresponds to the slow wind components. This structure is similar to the HCS, except for its inclination, which follows the strongly inclined dipole of BD-16351.

It is interesting to look at what happens to the flow just ahead of a density peak. The magnetic field is increasing, which means that a compression in the flow occurs. However, if the density increases at first, it decreases before the peak, compensated for by an acceleration of the flow. In these regions, the wind reaches its maximum Mach number, just below 5, as shown in red in Figure 8. The velocity then drops rapidly inside the streamer, where higher (by a factor of two to five) densities are achieved. This structure is repeated with the second streamer, beyond $\phi = 5.5$.

Nonetheless, the shock conditions are not yet fulfilled (see, for example, Gurnett & Bhattacharjee 2000), and the solution would likely require us to extend it in a larger domain to produce shocks. In our simulations, the discontinuities are separated by a contact layer with no matter exchanged between them. From a numerical point of view, the HLLE solver that we use is one of the most diffusive approximate Riemann solvers and may not be able to correctly capture such an oblique shock. It is likely, though, that shocks will form much earlier in the case of fast rotators, simply because of the enhanced helicity of the Parker spiral, thus changing the energetic budget in the astrospheres of young suns.

CIRs in the solar system are thought to be the most frequent cause for geomagnetic storms (Yermolaev et al. 2012). Our solar case does not show such features because we lack an additional heating mechanism for the acceleration of the fast wind. Yet, we can extrapolate our results stating that younger stars are likely to create more CIRs that shock in the interplanetary medium, adding up to other dynamical events that are thought to be enhanced in the environment of active stars, such as flares and coronal mass ejections (Schrijver et al. 2012).

## 5. DISCUSSIONS AND PERSPECTIVES

Our study addresses 3D simulations constrained by spectro-polarimetric observations of magnetic fields in the context of stellar evolution, along the main sequence. Our findings can be divided into two parts. First, considering the global and integrated properties, the mass and angular momentum loss
<figure><img src="image_19.png" /><figcaption>**Figure 8.** Color maps of the Mach number and the density logarithm in the equatorial plane of BD-16351. Magnetic field lines are plotted in black and create a Parker spiral. The fast wind component coming from coronal holes encounters the slow wind component created by streamers, and compression regions are visible at the edge of the domain.</figcaption></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 9.** Profiles of the density, Mach number, magnetic field, and current amplitude interpolated along a 25 $R_*$ equatorial circle in the BD-16351 case. Peaks of density associated with troughs of velocity mark the slow wind component. Compression and acceleration occur when the fast wind encounters the slow wind component.</figcaption></figure>

follow simple prescriptions thanks to an appropriate modeling of the evolution of the temperature and the coronal density with the rotation rate. An angular momentum proportional to the $\Omega_*^3$ can be obtained using our prescriptions with observed magnetic field amplitudes, which is required to observe the convergence of spin rates on the Skumanich law.

Our 3D simulations follow the braking law we derived in Réville et al. (2015a) if the $K_3$ constant is reduced by 15%. This can be understood because our simulations are made with a larger (and different for each case) coronal temperature compared to this previous work. Hence, the wind described is here faster, and for a given magnetic field strength, the Alfvén radius is closer to the star. However, the fit shows little deviation, and one constant $K_3$ in the scaling law is able to describe the whole range of temperature of our sample. A more general braking law should quantify the influence of the

temperature on $K_3$, since the variation remains limited in this study. Also, the use of a fully 3D geometry could be involved in the variation of $K_3$. Nonetheless, because this formulation expresses the dependence of the Alfvén radius on a magnetization parameter that includes the thermodynamics of the wind, it is likely to be valid for a wide range of magnetic fields, rotation rates, coronal temperatures, and densities, if one allows a small dependence of $K_3$ on the temperature.

With this adapted formulation, the semianalytical model we developed in Réville et al. (2015b) can be applied. The estimation method of the open flux of the simulation with a potential extrapolation that was tested on 2.5D configurations is perfectly operational with 3D nonaxisymmetric fields. Our semianalytical model is consequently able to estimate closely the evolution of angular momentum with the rotational period, as long as the mass loss rate of the simulations does not deviate too much from the spherically symmetric value used in the semianalytical model. This deviation grows as the stars rotate faster and possess more intense magnetic fields. Large coronal loops are able to confine more plasma, and the mass loss seems to plateau for $\Omega_* > 8\Omega_⊙$.

This behavior has consequences on the angular momentum loss that shows signs of saturation beyond this rotation rate, whose value is coherent with the saturation value used in rotation evolution models for K-type stars (Gallet & Bouvier 2015). Although the saturation of angular momentum loss is often associated with the saturation observed in the X-ray fluxes, the precise process behind this saturation remains unknown. Some works have suggested a stochastic change of the dynamo process generating the magnetic field could be involved in a topology switch from small scales to large scales that turns on the $\Omega_*^3$ braking law (Barnes 2003; Brown 2014; Garraffo et al. 2015b). However, no such transition is observed in our sample, as all our stars possess a strong dipolar field. In our simulations, the AML saturation seems to be due to the confinement of the outflow in large coronal loops that reduces the mass loss, which can be associated with the dependence of the wind braking on the filling factor (see Cranmer & Saar 2011; Gallet & Bouvier 2013, 2015). These results need,
however, to be confirmed by more simulations of fast rotators and are likely to be highly dependent on the prescriptions we used for the coronal base densities and temperature.

The mass loss rate of young stars in our study, although up to six to nine times the solar one,<fnref n="9" /> does not reach the highest values derived in Wood et al. (2005). A much higher dependence on rotation, for either the temperature or the coronal density, would have been necessary to observe 100 $\dot{M}_{\odot}$ values in our simulations. Our semianalytical model could be used to study the influence of different prescriptions on the variation of $\dot{J}$. A change in the exponents of the evolution laws, or the solar initial values, could be tested. A more physical description of the stellar wind acceleration, driven, for example, by Alfvén wave turbulence (see Suzuki 2006, 2013; Matsumoto & Suzuki 2012; Oran et al. 2013; Sokolov et al. 2013; Lionello et al. 2014), can help us to understand how coronal parameters evolve with age. Future works will be dedicated to including such processes in our simulations.

To solve this issue, more observations are also critical. In the work of Wood et al. (2005) (see also Wood 2004), the analysis of the Ly$\alpha$ absorption spectra is coupled with numerical simulations of the terminal shock, where the wind speed is kept constant at the slow solar wind value around 400 km s$^{-1}$. Our study, and this is the second part of our findings, brings more accurate constraints on the wind velocity amplitudes and distribution for solar-like stars that could be used to improve those calculations.

Indeed, several studies have now shown that the wind speeds are likely to increase for young stars, mostly because of higher coronal temperature and magnetocentrifugal acceleration (Washimi & Shibata 1993; Holzwarth & Jardine 2007; Matt et al. 2012; Suzuki et al. 2013; Réville et al. 2015a, 2015b). We show that the speed distribution of young and active stars follows a trimodal structure due to the interaction with the magnetic field. Moreover, the one-dimensional magnetocentrifugal wind solution (Weber & Davis 1967; Sakurai 1985) is at best a low estimate of the fastest component of the wind, as other magnetic processes are able to accelerate the wind. For instance, we observe in our simulations fast streams in the vicinity of the star caused by latitudinal or longitudinal superradial expansion of flux tubes due to the fast rotation and the nonaxisymmetry of the magnetic field.

As slow and fast winds exist around solar-like stars, they can interact in the equatorial plane, creating CIRs, which could be more common for younger stars and could be forming shocks closer to the star (the usual distance observed in the solar system is between 1 and 5 AU). This could have consequences on the energetics of expanding stellar winds of young stars and on exoplanetary space weather. The current sheet of young stars is also strongly corrugated when the dipole is inclined and polarity variations occur closer to one another, which has important consequences for the interaction with close-in planets, especially when they are within the Alfvén surface (Strugarek et al. 2015). This is, however, a very preliminary study, and more focused work needs to be done in that sense, improving, for instance, the numerical scheme and shock-absorbing method.

We thank Sean Matt and Jérôme Bouvier for continuous and fruitful discussions. We acknowledge funding by the ANR Blanc TOUPIES SIMI5-6 020 01, the ERC STARS2 207430, and CNES via Solar Orbiter. High-performance computations

were performed on the French machines Turing (IDRIS) and Curie (TGCC), within the GENCI 1623 program. We are grateful to Andrea Mignone and his team at University of Torino for the development and the maintenance of the PLUTO code.

AS is a National Postdoctoral Fellow at the Canadian Institute of Theoretical Astrophysics (CITA) and acknowledges support from the Natural Sciences and Engineering Research Council of Canada.

## APPENDIX

### FLUX TUBE EXPANSION IN THE SUPERSONIC REGIME

The mass conservation within a flux tube reduces to the equation

$$ \frac{1}{v} \frac{dv}{dr} + \frac{1}{\rho} \frac{d\rho}{dr} + \frac{1}{A} \frac{dA}{dr} = 0, \qquad (14) $$

where A is the section of the flux tube. The momentum equation, if we consider an isothermal flow for simplicity, can then be written as

$$ \left( v - \frac{c_s^2}{v} \right) \frac{dv}{dr} = \frac{c_s^2}{A} \frac{dA}{dr} - \frac{v_{\text{kep}}^2(r_c)r_c}{r^2}, \qquad (15) $$

where $c_s$ is the constant sound speed, $v_{\text{kep}}(r) = \sqrt{GM_*/r}$ is the Keplerian velocity, and $r_c$ is the critical radius where the wind becomes supersonic.

In terms of Mach number $M = u/c_s$, we can write

$$ \left( M - \frac{1}{M} \right) \frac{dM}{dr} = \frac{1}{A} \frac{dA}{dr} - \frac{v_{\text{kep}}^2(r_c)r_c}{c_s^2 r^2}. \qquad (16) $$

When the wind is subsonic ($M < 1$), the term $(M - 1/M)$ is negative, and thus what matters for the acceleration of the wind is the sign of the term on the right-hand side. The expansion of the flux tube can be locally described by $A \propto r^\alpha$ and yields $(1/A)dA/dr = \alpha/r$. The spherically symmetric solution, which corresponds to the case $A \propto r^2$, and $\alpha = 2$ gives

$$ \frac{2}{r} \le \frac{v_{\text{kep}}^2(r_c)r_c}{c_s^2 r^2}, \qquad (17) $$

as the Parker solution is always accelerating. It is shown that superradial expansion is globally anticorrelated with the wind speed (Wang & Sheeley 1990), because of the local inversion of the latter inequality below the sonic surface. Nonetheless, superradial expansion can accelerate the outflow as long as $(1/A)dA/dr$ remains smaller than $v_{\text{kep}}^2(r_c)r_c/(c_s^2 r^2)$, for instance if the superradial expansion is located near the surface or well below the sonic point (see Velli 2010).

However, in the supersonic regime, a superradial expansion will necessarily accelerate the outflow. Indeed, as $\alpha \ge 2$, and since the right-hand side is positive in the radial case, any superradial expansion will grow this term, and $dM/dr$ will be larger.

The treatment of this problem without the isothermal approximation is more complex (see Kopp & Holzer 1976). In our study, the sound speed variation may not be negligible, but the qualitative behavior remains, and a quantitative analysis is left for future works.

<footnote n="9">⁹ In Table 2, we have considered an upper value of the solar mass loss rate $\dot{M} = 3 \times 10^{-14} M_{\odot} \text{yr}^{-1}$.</footnote>
REFERENCES

Altschuler, M. D., & Newkirk, G. 1969, SoPh, 9, 131

Alvarado-Gómez, J. D., Hussain, G. A. J., Cohen, O., et al. 2016a, A&A, 588, A28

Alvarado-Gómez, J. D., Hussain, G. A. J., Cohen, O., et al. 2016b, A&A, 594, A95

Barnes, S. A. 2003, ApJ, 586, 464

Belcher, J. W., & Davis, L., Jr. 1971, JGR, 76, 3534

Bouvier, J., Forestini, M., & Allain, S. 1997, A&A, 326, 1023

Brandenburg, A., & Saar, S. H. 2000, in ASP Conf. Ser. 198, Stellar Clusters and Associations: Convection, Rotation, and Dynamics, ed. R. Pallavicini, G. Micela, & S. Sciortino (San Francisco, CA: ASP), 381

Brown, T. M. 2014, ApJ, 789, 101

Brun, A. S., García, R. A., Houdek, G., Nandy, D., & Pinsonneault, M. 2015, SSRv, 196, 303

Cohen, O., & Drake, J. J. 2014, ApJ, 783, 55

Cohen, O., Kashyap, V. L., Drake, J. J., et al. 2011, ApJ, 733, 67

Cranmer, S. R., & Saar, S. H. 2011, ApJ, 741, 54

DeRosa, M. L., Brun, A. S., & Hoeksema, J. T. 2012, ApJ, 757, 96

do Nascimento, J.-D., Jr., Vidotto, A. A., Petit, P., et al. 2016, ApJL, 820, L15

Donati, J.-F., & Brown, S. F. 1997, A&A, 326, 1135

Donati, J.-F., & Landstreet, J. D. 2009, ARA&A, 47, 333

Durney, B. R., & Latour, J. 1978, GApFD, 9, 241

Einfeldt, B. 1988, JCoPh, 25, 294, http://arrow.utias.utoronto.ca/~groth/ aer1319/Handouts/Additional_Reading_Material/SIAMNumAnal-1988-einfeldt.pdf

Ferreira, J. 2013, in EAS Publications Ser. 62, Role and Mechanisms of Angular Momentum Transport During the Formation and Early Evolution of Stars, ed. P. Hennebelle & C. Charbonnel (Les Ulis: EDP), 169

Folsom, C. P., Petit, P., Bouvier, J., et al. 2016, MNRAS, 457, 580

Gallet, F., & Bouvier, J. 2013, A&A, 556, A36

Gallet, F., & Bouvier, J. 2015, A&A, 577, A98

Garraffo, C., Drake, J. J., & Cohen, O. 2015a, ApJL, 807, L6

Garraffo, C., Drake, J. J., & Cohen, O. 2015b, ApJ, 813, 40

Giampapa, M. S. 2005, in Saas-Fee Advanced Course 34: The Sun, Solar Analogs and the Climate, ed. J. D. Haigh et al. (Berlin: Springer), 307

Güdel, M. 2004, A&Ar, 12, 71

Gurnett, D., & Bhattacharjee, A. 2000, Introduction do Plasma Physics: With Space and Laboratory Applications (Cambridge: Cambridge Univ. Press)

Holzwarth, V., & Jardine, M. 2007, A&A, 463, 11

Ivanova, N., & Taam, R. E. 2003, ApJ, 599, 516

Kawaler, S. D. 1988, ApJ, 333, 236

Keppens, R., & Goedbloed, J. P. 2000, ApJ, 530, 1036

Keppens, R., MacGregor, K. B., & Charbonneau, P. 1995, A&A, 294, 469

Kopp, R. A., & Holzer, T. E. 1976, SoPh, 49, 43

Linsky, J. L., & Wood, B. E. 1996, ApJ, 463, 254

Lionello, R., Velli, M., Downs, C., Linker, J. A., & Mikić, Z. 2014, ApJ, 796, 111

Matsumoto, T., & Suzuki, T. K. 2012, ApJ, 749, 8

Matt, S., & Balick, B. 2004, ApJ, 615, 921

Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJL, 799, L23

Matt, S. P., Pinzón, G., Greene, T. P., & Pudritz, R. E. 2012, ApJ, 745, 101

Mestel, L. 1968, MNRAS, 138, 359

Michel, F. C. 1969, ApJ, 158, 727

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

Noyes, R. W., Weiss, N. O., & Vaughan, A. H. 1984, ApJ, 287, 769

Oran, R., van der Holst, B., Landi, E., et al. 2013, ApJ, 778, 176

Pallavicini, R. 1989, A&Ar, 1, 177

Pallavicini, R., Golub, L., Rosner, R., et al. 1981, ApJ, 248, 279

Parker, E. N. 1958, ApJ, 128, 664

Pinto, R. F., Brun, A. S., Jouve, L., & Grappin, R. 2011, ApJ, 737, 72

Piskunov, N., & Kochukhov, O. 2002, A&A, 381, 736

Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147

Powell, K. G. 1994, NASA Technical Report, ICASE Report 94-24

Randich, S. 2000, in ASP Conf. Ser. 198, Stellar Clusters and Associations: Convection, Rotation, and Dynamics, ed. R. Pallavicini, G. Micela, & S. Sciortino (San Francisco, CA: ASP), 401

Réville, V., Brun, A. S., Matt, S. P., Strugarek, A., & Pinto, R. F. 2015a, ApJ, 798, 116

Réville, V., Brun, A. S., Strugarek, A., et al. 2015b, ApJ, 814, 99

Richardson, I. G. 2004, SSRv, 111, 267

Sakurai, T. 1985, A&A, 152, 121

Schatten, K. H., Wilcox, J. M., & Ness, N. F. 1969, SoPh, 6, 442

Schatzman, E. 1962, AnAp, 25, 18

Schrijver, C. J., Beer, J., Baltensperger, U., et al. 2012, JGRA, 117, A08103

See, V., Jardine, M., Vidotto, A. A., et al. 2015, MNRAS, 453, 4301

Semel, M. 1989, A&A, 225, 456

Smith, E. J., & Wolfe, J. H. 1976, GeoRL, 3, 137

Sokól, J. M., Swaczyna, P., Bzowski, M., & Tokumaru, M. 2015, SoPh, 290, 2589

Sokolov, I. V., van der Holst, B., Oran, R., et al. 2013, ApJ, 764, 23

Strugarek, A., Brun, A. S., Matt, S. P., & Réville, V. 2015, ApJ, 815, 111

Suzuki, T. K. 2006, ApJL, 640, L75

Suzuki, T. K. 2013, AN, 334, 81

Suzuki, T. K., Imada, S., Kataoka, R., et al. 2013, PASJ, 65, 98

ud-Doula, A., & Owocki, S. P. 2002, ApJ, 576, 413

Ustyugova, G. V., Koldoba, A. V., Romanova, M. M., Chechetkin, V. M., & Lovelace, R. V. E. 1999, ApJ, 516, 221

van Saders, J. L., Ceillier, T., Metcalfe, T. S., et al. 2016, Natur, 529, 181

Velli, M. 1994, ApJL, 432, L55

Velli, M. 2010, Twelfth International Solar Wind Conference, 1216, 14

Vidotto, A. A., Donati, J.-F., Jardine, M., et al. 2016, MNRAS, 455, L52

Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014a, MNRAS, 441, 2361

Vidotto, A. A., Jardine, M., Morin, J., et al. 2014b, MNRAS, 438, 1162

Wang, Y.-M., & Sheeley, N. R., Jr. 1990, ApJ, 355, 726

Washimi, H., & Shibata, S. 1993, MNRAS, 262, 936

Weber, E. J., & Davis, L., Jr. 1967, ApJ, 148, 217

Weiss, N. O. 1994, in Lectures on Solar and Planetary Dynamics, ed. M. R. E. Proctor & A. D. Gilbert (Cambridge: Cambridge Univ. Press), 59

Wood, B. E. 2004, LRSP, 1, 2

Wood, B. E., Müller, H.-R., Zank, G. P., & Linsky, J. L. 2002, ApJ, 574, 412

Wood, B. E., Müller, H.-R., Zank, G. P., Linsky, J. L., & Redfield, S. 2005, ApJL, 628, L143

Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48

Yermolaev, Y. I., Lodkina, I. G., Nikolaeva, N. S., & Yermolaev, M. Y. 2012, JGRA, 117, A08207

Zanni, C., & Ferreira, J. 2009, A&A, 508, 1117