# Alfvén-wave-driven Magnetic Rotator Winds from Low-mass Stars. I. Rotation Dependences of Magnetic Braking and Mass-loss Rate

Munehito Shoda, Takeru K. Suzuki, Sean P. Matt, Steven R. Cranmer, Aline

A. Vidotto, Antoine Strugarek, Victor See, Victor Réville, Adam J. Finley, Allan Sacha Brun

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

## ► To cite this version:

Munehito Shoda, Takeru K. Suzuki, Sean P. Matt, Steven R. Cranmer, Aline A. Vidotto, et al., Alfvén-wave-driven Magnetic Rotator Winds from Low-mass Stars. I. Rotation Dependences of Magnetic Braking and Mass-loss Rate. The Astrophysical Journal, 2020, 896, 10.3847/1538-4357/ab94bf.insu-03636588

HAL Id: insu-03636588

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

Submitted on 6 Aug 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_2.png" /></figure>

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

CrossMark

# Alfvén-wave-driven Magnetic Rotator Winds from Low-mass Stars. I. Rotation Dependences of Magnetic Braking and Mass-loss Rate

Munehito Shoda<fnref n="1" />, Takeru K. Suzuki<fnref n="2" />, Sean P. Matt<fnref n="3" />, Steven R. Cranmer<fnref n="4" />, Aline A. Vidotto<fnref n="5" />, Antoine Strugarek<fnref n="6" />, Victor See<fnref n="3" />, Victor Réville<fnref n="7" />, Adam J. Finley<fnref n="3" />, and Allan Sacha Brun<fnref n="6" />

Received 2020 March 4; revised 2020 April 30; accepted 2020 May 18; published 2020 June 19

## Abstract

Observations of stellar rotation show that low-mass stars lose angular momentum during the main sequence. We simulate the winds of sunlike stars with a range of rotation rates, covering the fast and slow magneto-rotator regimes, including the transition between the two. We generalize an Alfvén-wave-driven solar wind model that builds on previous works by including the magneto-centrifugal force explicitly. In this model, the surface-averaged open magnetic flux is assumed to scale as $B_*f_*^{\text{open}} \propto \text{Ro}^{-1.2}$, where $f_*^{\text{open}}$ and $\text{Ro}$ are the surface open-flux filling factor and Rossby number, respectively. We find that, (1) the angular-momentum loss rate (torque) of the wind is described as $\tau_w \approx 2.59 \times 10^{30} \text{ erg} (\Omega_*/\Omega_\odot)^{2.82}$, yielding a spin-down law $\Omega_* \propto t^{-0.55}$. (2) The mass-loss rate saturates at $\dot{M}_w \sim 3.4 \times 10^{-14} M_\odot \text{ yr}^{-1}$, due to the strong reflection and dissipation of Alfvén waves in the chromosphere. This indicates that the chromosphere has a strong impact in connecting the stellar surface and stellar wind. Meanwhile, the wind ram pressure scales as $P_w \propto \Omega_*^{0.57}$, which is able to explain the lower envelope of the observed stellar winds by Wood et al. (3) The location of the Alfvén radius is shown to scale in a way that is consistent with one-dimensional analytic theory. Additionally, the precise scaling of the Alfvén radius matches previous works, which used thermally driven winds. Our results suggest that the Alfvén-wave-driven magnetic rotator wind plays a dominant role in the stellar spin-down during the main sequence.

*Unified Astronomy Thesaurus concepts: Stellar evolution (1599); Stellar rotation (1629); Stellar winds (1636); Solar wind (1534); Solar evolution (1492); Stellar mass loss (1613); Magnetohydrodynamical simulations (1966)*

## 1. Introduction

The dynamo process yields stellar magnetic fields (Leighton 1969; Brun et al. 2004; Hotta et al. 2016) that give rise to activity such as coronal heating (Alfvén 1947; Osterbrock 1961; Parker 1988; Rappazzo et al. 2008), stellar winds (Parker 1958; Velli 1994), flares, and coronal mass ejections (CMEs; Argiroffi et al. 2019; Notsu et al. 2019; Toriumi & Wang 2019). Stellar activity of this type is observed to decay over the lifetime of a star (Skumanich 1972; Güdel et al. 1997; Güdel 2007; Vidotto et al. 2014b). Understanding such long-term evolution is one of the most important challenges in astronomy, especially in the context of stellar influences on the habitability of exoplanets, as the erosion of planetary atmospheres is affected by stellar activity (Lammer et al. 2010; Johnstone et al. 2015b; Garraffo et al. 2016; Allan & Vidotto 2019; Johnstone et al. 2019; Airapetian et al. 2020; Vidotto & Cleary 2020). In order to understand activity evolution, we first need to understand the evolution of stellar rotation since (differential) rotation and convection are the ultimate origin of magnetic energy in low-mass stars.

It is widely known that low-mass stars spin down over their lifetimes (Schatzman 1962; Kraft 1967), approximately as $\Omega \propto t^{-1/2}$ over the age range of $\sim 10^8$ yr to the age of the Sun (Skumanich 1972). This stellar spin-down is due to the angular-momentum loss caused by magnetized stellar winds (magnetic braking; Weber & Davis 1967; Sakurai 1985; Kawaler 1988).

Magnetic braking governs the long-term variations in stellar rotation, and thus the stellar dynamo process, which in turn affects the intensity and structure of the stellar wind. In this way, the interplay between the stellar dynamo and stellar wind regulates the rotational evolution of stars (Brun & Browning 2017).

As well as being an indicator for dynamo efficiency, stellar rotation is an important fundamental quantity that can be used as a stellar age diagnostic. Since magnetic braking is stronger for faster rotators, the rotation periods of low-mass stars are often found to converge onto a sequence defined by mass and age, regardless of their initial rotation rates (Irwin & Bouvier 2009). For example, stars with $M_* \gtrsim 0.5 M_\odot$ are known to have rotationally converged by the age of the Hyades cluster (Radick et al. 1987; Delorme et al. 2011). The use of rotation as a proxy for age in this way is known as gyrochronology (Barnes 2003, 2007, 2010). An alternative diagnostic based on magnetic field strength instead of rotation rate has also been proposed (Vidotto et al. 2014b). Gyrochronology mainly appears to be applicable to middle-aged stars, i.e., $t_{age} \lesssim 2.5$ Gyr for sunlike stars (Meibom et al. 2015), while recent asteroseismic studies suggest that the stellar age–color–rotation relation may deviate from gyrochronology for stars older than the age of the Sun (Angus et al. 2015; Davies et al. 2015; van Saders et al. 2016). To understand what causes the breakdown of gyrochronology, we first need to correctly model the mechanism by which stars lose angular momentum.

<footnote n="1">¹National Astronomical Observatory of Japan, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588, Japan; munehito.shoda@nao.ac.jp</footnote>

<footnote n="2">²School of Arts & Science, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo, 153-8902, Japan</footnote>

<footnote n="3">³University of Exeter, Department of Physics & Astronomy, Physics Bldg., Stocker Road, Exeter EX4 4QL, UK</footnote>

<footnote n="4">⁴Department of Astrophysical and Planetary Sciences, Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, CO 80309, USA</footnote>

<footnote n="5">⁵School of Physics, Trinity College Dublin, The University of Dublin, Dublin-2, Ireland</footnote>

<footnote n="6">⁶Department of Astrophysics & UMR AIM, CEA Paris-Saclay, CNRS/INSU, University of Paris 7, Gif-sur-Yvette, France</footnote>

<footnote n="7">⁷IRAP, Université Toulouse III—Paul Sabatier, CNRS, CNES, Toulouse, France</footnote>
The evolution of stellar rotation periods is governed by several physical processes, such as disk-locking, core-envelope decoupling, internal-structure evolution, and magnetic braking (Gallet & Bouvier 2013, 2015). The magnetic braking plays a dominant role in the net angular-momentum loss during the main sequence. In contrast to observation-based approach to mass-loss (Johnstone et al. 2015a; Ahuir et al. 2020) and angular-momentum loss rate (Matt et al. 2015), we aim to model them in a physics-based way. One problem of physics-based stellar-wind models is that the scaling laws of the mass-loss rate (Schröder & Cuntz 2005; Holzwarth & Jardine 2007; Suzuki 2007, 2018; Cranmer & Saar 2011) and the Alfvén radius (Kawaler 1988; Matt & Pudritz 2008; Matt et al. 2012; Réville et al. 2015; Finley & Matt 2017, 2018) have been discussed independently (note that the torque is a function of the mass-loss rate, Alfvén radius, and rotation rate; see Weber & Davis 1967). However, both the mass-loss rate and the Alfvén radius vary with stellar-wind density, and therefore should be modeled simultaneously.

The mass-loss rate is determined by the energy balance in the chromosphere and the corona, while the Alfvén radius is related to the large-scale magnetism of the star and stellar-wind acceleration. Thus, in order to simultaneously model the mass-loss rate and Alfvén radius, we need (1) to resolve the chromosphere and waves therein (typical spatial scale ~a few 100 km) and (2) a simulation domain that is sufficiently large to cover the wind acceleration (typical spatial scale ~a few $10R_*$ or more), which requires typically $10^{-5}$ grid points in the radial direction. For this reason, we make use of a one-dimensional solar wind model that satisfies the aforementioned demand and generalize it to stellar wind by explicitly taking into account the rotation effect. This model allows us to investigate the dependence of stellar-wind parameters (mass-loss rate, Alfvén radius, torque) on the stellar rotation rate. A goal of this work is to derive the rotation dependence of stellar-wind characteristics and compare them with observations.

The remainder of this paper is organized as follows. In Section 2, we summarize the overview of the model in this work, including assumptions, basic equations, parameters, and numerical schemes. The numerical results are discussed in Section 3. The energetics of the stellar wind is discussed in Section 4, focusing on the mass-loss saturation and wave energetics. We discuss the overall results of our work in Section 5.

## 2. Model

### 2.1. Overview

We simulate equatorial stellar winds that extend from the stellar photosphere to beyond the fast-magnetosonic point. Our model is based on the magnetohydrodynamic equations including gravity, thermal conduction, and radiative cooling. For simplicity, and to reduce the numerical cost, we assume a one-dimensional geometry and axisymmetry. We therefore include the turbulent dissipation of Alfvén waves, which is a multidimensional effect, phenomenologically.

Several theoretical models explain the solar wind based on Alfvén-wave heating and acceleration (Suzuki & Inutsuka 2006; Cranmer et al. 2007; van der Holst et al. 2014; Shoda et al. 2019; Réville et al. 2020). We extend this Alfvén-wave modeling to winds from low-mass stars. In addition to Alfvén waves, for fast rotators the magneto-centrifugal force can

further accelerate the wind (Belcher & MacGregor 1976; Sakurai 1985; Réville et al. 2016; Johnstone 2017). In this work, we account for both effects and conventionally call our models “Alfvén-wave-driven magnetic rotator winds.” For simplicity, we fix the mass, luminosity, and metallicity of the star with solar values ($M_* = M_\odot$, $L_* = L_\odot$, $Z_* = Z_\odot$) and focus our interest on the rotation dependence.

A key factor in our model (and in the theory of magnetic braking in general) is the filling factor of open magnetic flux, defined as

$$f^{\text{open}}(r) = \frac{\Phi_{\text{open}}}{4\pi r^2 |B_r(r)|}, \quad (1)$$

where $|B_r(r)|$ is the unsigned radial magnetic field (averaged over solid angle) and $\Phi_{\text{open}}$ is the unsigned open magnetic flux. Note that $\Phi_{\text{open}}$ is constant in $r$ and is fixed for each simulation. Magnetic field lines eventually become open, due to their advection in the stellar wind; therefore, $f^{\text{open}}(r \to \infty) = 1$. At the stellar surface, $f^{\text{open}}$ is generally much lower than unity. For the solar surface, $f^{\text{open}}$ is typically $10^{-3}$ (Cranmer 2017). The radial increase of $f^{\text{open}}$ as a function of radius, i.e., the (super-radial) expansion of open magnetic field line, needs to be accounted for.

To follow the (super-radial) magnetic field expansion in our one-dimensional geometry, we make use of a field-aligned coordinate system (Hollweg et al. 1982; Kudoh & Shibata 1999; Suzuki & Inutsuka 2005). To be consistent with axisymmetry, all of the super-radial expansion is attributed to the poloidal ($r$ and $\theta$) components. The scale factors (that reflect the degree of expansion of magnetic flux in each direction) of the corresponding curvilinear coordinate system are given as

$$h_r = 1, \quad h_\theta = r f^{\text{open}}, \quad h_\phi = r. \quad (2)$$

For simplicity, the $\theta$ and $\phi$ components are attributed to (Alfvén) waves and rotation, respectively. By this simplification, the polarization of Alfvén waves is restricted to be linear. However, this restriction is unlikely to affect the conclusion because the wind structure and dynamics are weakly affected by the imposed polarization (Suzuki & Inutsuka 2006). Although stellar-wind outflows are far from symmetric (e.g., van der Holst et al. 2014), once the flux-tube expansion is appropriately prescribed by $f^{\text{open}}(r)$, our one-dimensional model is expected to give a reasonable estimation of wind parameters. In the wind acceleration region and below, due to the low-beta nature of the corona, the flux-tube expansion is essentially set by the global magnetic field and is not affected by the wind dynamics. Because the interactions between flux tubes are likely insignificant, in wave-driven winds, each flux tube behaves independently. Indeed, the three-dimensional structure of the solar wind is well reproduced from an ensemble of one-dimensional flux-tube models (Pinto & Rouillard 2017). Thus, by implementing a representative flux-tube expansion, we can reliably recover averaged properties of the stellar wind.

An overview of our model geometry is detailed in Figure 1. An equatorial magnetic flux tube is located on the stellar surface and expands super-radially into interplanetary space. MHD waves propagate along the background flux tube and partially dissipate in the atmosphere.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** A schematic picture of the stellar-wind geometry used in this study. Shown by the black symbols and characters are the numerical settings. The red symbols and characters refer to the physical processes considered in this work.</figcaption></figure>

## 2.2. Basic Equations

We assume that the system is one-dimensional ($\partial/\partial \theta = \partial/\partial \phi = 0$) and the scale factors are given by Equation (2). The MHD equations are then written as follows (see Appendix A for derivation):

$$ \frac{\partial}{\partial t} \rho + \frac{1}{r^2 f^{\text{open}}} (\rho v_r r^2 f^{\text{open}}) = 0, \quad (3) $$

$$ \begin{aligned} & \frac{\partial}{\partial t}(\rho v_r) + \frac{1}{r^2 f^{\text{open}}} \frac{\partial}{\partial r}[(\rho v_r^2 + p_T) r^2 f^{\text{open}}] \\ & = -\rho \frac{GM_*}{r^2} + p \frac{d}{dr} \ln(r^2 f^{\text{open}}) \\ & + \rho v_\theta^2 \frac{d}{dr} \ln(rf^{\text{open}}) + \frac{1}{r} \rho v_\phi^2 - \left(\frac{B_\theta^2}{8\pi} - \frac{B_\phi^2}{8\pi}\right) \frac{d}{dr} \ln f^{\text{open}}, \quad (4) \end{aligned} $$

$$ \begin{aligned} & \frac{\partial}{\partial t}(\rho v_\theta) + \frac{1}{r^2 f^{\text{open}}} \frac{\partial}{\partial r} \left[ \left( \rho v_r v_\theta - \frac{1}{4\pi} B_r B_\theta \right) r^2 f^{\text{open}} \right] \\ & = \left( \frac{B_r B_\theta}{4\pi} - \rho v_r v_\theta \right) \frac{d}{dr} \ln(rf^{\text{open}}) + \rho D_v^{\text{turb}}, \end{aligned} \quad (5) $$

$$ \begin{aligned} & \frac{\partial}{\partial t}(\rho v_\phi) + \frac{1}{r^2 f^{\text{open}}} \frac{\partial}{\partial r} \left[ \left( \rho v_r v_\phi - \frac{1}{4\pi} B_r B_\phi \right) r^2 f^{\text{open}} \right] \\ & = \left( \frac{B_r B_\phi}{4\pi} - \rho v_r v_\phi \right) / r, \end{aligned} \quad (6) $$

$$ \frac{1}{r^2 f^{\text{open}}} \frac{d}{dr} (B_r r^2 f^{\text{open}}) = 0, \quad (7) $$

$$ \begin{aligned} & \frac{\partial}{\partial t} B_\theta + \frac{1}{r^2 f^{\text{open}}} \frac{\partial}{\partial r} \left[ (v_r B_\theta - v_\theta B_r) r^2 f^{\text{open}} \right] \\ & = (v_r B_\theta - v_\theta B_r) \frac{d}{dr} \ln(rf^{\text{open}}) + \sqrt{4\pi\rho} D_b^{\text{turb}}, \end{aligned} \quad (8) $$

$$ \begin{aligned} & \frac{\partial}{\partial t} B_\phi + \frac{1}{r^2 f^{\text{open}}} \frac{\partial}{\partial r} \left[ (v_r B_\phi - v_\phi B_r) r^2 f^{\text{open}} \right] \\ & = (v_r B_\phi - v_\phi B_r) / r, \end{aligned} \quad (9) $$

$$ \begin{aligned} & \frac{\partial}{\partial t} e + \frac{1}{r^2 f^{\text{open}}} \frac{\partial}{\partial r} \left[ \left( (e + p_T) v_r - \frac{B_r}{4\pi} (v_\perp \cdot \mathbf{B}_\perp) + F_C \right) r^2 f^{\text{open}} \right] \\ & = -\rho v_r \frac{GM_*}{r^2} - Q_R, \end{aligned} \quad (10) $$

where

$$ e = \frac{p}{\gamma - 1} + \frac{1}{2} \rho v^2 + \frac{B_\perp^2}{8\pi}, \quad p_T = p + \frac{B_\perp^2}{8\pi}. \quad (11) $$

These are closed by the equation of state.

$$ p = c_1 \rho T, \quad c_1 = 1.36 \times 10^8 \text{ erg g}^{-1} \text{ K}^{-1} \quad (12) $$

where the value of $c_1$ is consistent with the equation of state of a fully ionized plasma that contains a few percent helium (alpha particle) by number. $D_v^{\text{turb}}$ and $D_b^{\text{turb}}$ represent the rate of turbulent dissipation of Alfvén wave per unit momentum (Shoda et al. 2018a), which will be discussed in Section 2.4. $F_C$ and $Q_R$ are the conductive flux and radiative cooling rate, respectively.

We employ a Spitzer–Härm type of thermal conductive flux (Spitzer & Härm 1953), with a quenching term that works at radial distances typically greater than $5R_\odot$:

$$ F_C = -\min \left( 1, \frac{\rho}{\rho_C} \right) \frac{B_r}{|\mathbf{B}|} \kappa_0 T^{5/2} \frac{dT}{dr}, \quad (13) $$

where we set $\kappa_0 = 10^{-6} \text{ erg cm}^{-1} \text{ s}^{-1} \text{ K}^{-7/2}$ and $\rho_C = 10^{-20} \text{ g cm}^{-3}$. The quenching term comes from the saturation of heat flux in interplanetary space (Salem et al. 2003; Bale et al. 2013). Although this quenching is known to be overestimated, it is unlikely to affect the numerical result as it generally occurs beyond the sonic point.

The radiative cooling is a combination of different types:

$$ Q_R = Q_{R,\text{thck}} \xi_1 + Q_{R,\text{thin}} (1 - \xi_1), \quad (14) $$

where $Q_{R,\text{thck}}$ and $Q_{R,\text{thin}}$ stand for the optically thick and thin radiative losses, respectively. The switching parameter $\xi_1$
mimics the optical depth, which takes $\xi_1 \approx 1$ in the photosphere and $\xi_1 \approx 0$ in the corona. Here we assume the following expression for $\xi_1$:

$$ \xi_1 = \max \left( 0, 1 - \frac{p_{chr}}{p} \right), \qquad (15) $$

where $p_{chr} = 1$ dyn cm$^{-2}$.

In the photosphere, where the optical depth is large, the balance between radiative heating and cooling keeps the temperature almost fixed. For this reason, following Gudiksen & Nordlund (2005), we approximate the optically thick cooling by an exponential cooling:

$$ Q_{R,thck} = \frac{1}{\tau_{thck}}(e_{int} - e_{int,ref}), \qquad (16) $$

$$ \tau_{thck} = 0.1 \left( \frac{\rho}{\bar{\rho}_*} \right)^{-5/3} s, \qquad (17) $$

where $\bar{\rho} = 10^{-7}$ g cm$^{-3}$ is the mean surface density, and $e_{int,ref}$ is the internal energy at a given reference temperature that mimics the radiation-balanced profile (e.g., Figure 4 in Cranmer & Winebarger 2019). $Q_{R,thck}$ only works near the surface because of the rapid increase of $\tau_{thck}$ with height.

Following Iijima (2016), the optically thin cooling function is composed of two different contributions. In the chromosphere, we employ the radiative cooling function given by Goodman & Judge (2012; $Q_{GJ}$), while in the corona, the optically thin cooling function taken from Rempel (2017) is used. These two functions are smoothly connected as a function of temperature using:

$$ Q_{R,thin} = Q_{GJ}(\rho, T)\xi_2 + n_p n_e \Lambda(T)(1 - \xi_2), \qquad (18) $$

$$ \xi_2 = \max \left( 0, \min \left( 1, \frac{T_{TR} - T}{\Delta T} \right) \right), \qquad (19) $$

where $T_{TR} = 15,000$ K and $\Delta T = 5000$ K.

## 2.3. Open-flux Filling Factor

There is a strong relationship between the open magnetic flux, $\Phi_{open}$, and the strength of magnetic braking (Vidotto et al. 2012, 2014a; Réville et al. 2015). The magnetic-flux conservation yields

$$ \Phi_{open} = 4\pi R_*^2 B_* f_*^{open} \approx 4\pi R_*^2 B_{eq,*} f_*^{open}, \qquad (20) $$

where $B_*$ is the characteristic field strength at the photosphere, and can be approximated by the equipartition value $B_{eq,*}$ that represents equal gas and magnetic pressures (Cranmer & Saar 2011). We assume that the stellar surface is divided into two areas: one with zero field and the other with equipartition field. Indeed, the Sun's photospheric magnetic field is observed to be spatially localized, exhibiting a nearly equipartition value (Tsuneta et al. 2008). Under this assumption, $f_*^{open}$ represents the fraction of the stellar surface covered by open magnetic flux (Saar 2001; Reiners et al. 2009). We note that the “open-flux filling factor” does not stand for the fraction of the open flux to the total flux $f_*^{open}/(f_*^{open} + f_*^{closed})$, where $f_*^{closed}$ is the fraction of the stellar surface covered by closed magnetic flux. We also note that, for the solar wind, $f_*^{open}$ (or, equivalently, the expansion factor) plays a role in determining the wind

speed (Wang & Sheeley 1990; Arge & Pizzo 2000; Fujiki et al. 2015).

Unfortunately, there is no established way to determine $f_*^{open}$ from the photospheric magnetic field, even for the Sun. For example, the widely used potential field source surface model (Schatten et al. 1969) consistently underestimates the open magnetic flux observed in situ spacecraft, which is referred to as the open-flux problem (Linker et al. 2017). However, it is thought that the dipolar magnetic field is the most significant contributor to the open magnetic flux (Réville et al. 2015; See et al. 2018). Recently, See et al. (2019) showed that, for most stars (especially with low Rossby numbers), the dipolar magnetic field is sufficient to determine the angular-momentum loss rate. For these reasons, we simply assume that the open magnetic flux is proportional to the surface-averaged unsigned dipolar magnetic field $\langle B_{dip} \rangle$:

$$ f_*^{open} \propto \frac{\langle B_{dip} \rangle}{B_{eq,*}}. \qquad (21) $$

In this work, we assume a power-law relation for $f_*^{open}$ as

$$ f_*^{open} = f_{\odot}^{open} \left( \frac{Ro}{Ro_{\odot}} \right)^{-1.2} \approx 10^{-3} \left( \frac{P_{rot}}{P_{rot,\odot}} \right)^{-1.2}, \qquad (22) $$

where we use $f_{\odot}^{open} \approx 10^{-3}$ as the the solar value of $f_*^{open}$. Note that we only consider main-sequence sunlike stars ($M_* = M_{\odot}$, $R_* = R_{\odot}$, $Z_* = Z_{\odot}$), and therefore, the Rossby number is a function of rotation rate only. Our implementation is in line with the derivation presented in See et al. (2019), who showed

$$ \langle B_{dip} \rangle \propto Ro^{-1.3 \pm 0.1}, \qquad (23) $$

based on a statistical analysis of Zeeman–Doppler imaging (ZDI) observations. The actual dependence of $\langle B_{dip} \rangle$ on Ro could be weaker because the ZDI observation tends to underestimate the magnetic field strength, especially when the field is weak (See et al. 2020). This supports our assumption that $\langle B_{dip} \rangle$ depends more weakly on Ro than the observation by See et al. (2019). In the future, a rotation-dependent correction factor to the ZDI-based $\langle B_{dip} \rangle$ values could be used. We summarize the input and output parameters of our simulations in Table 1.

We assume a two-step super-radial expansion of the magnetic field line: one expansion occurs in the chromosphere and the other in the corona (Cranmer & van Ballegooijen 2005). To implement such two-step expansion, we need to set the filling factor between the two expansion regions (at the coronal base), which we denote $f_{cor}^{open}$. Following Cranmer & Saar (2011), we simply assume a power-law relation between $f_*^{open}$ and $f_{cor}^{open}$ as

$$ f_{cor}^{open} = (f_{*}^{open})^{\theta_B}, \qquad (24) $$

where we use $\theta_B = 1/3$ as a reference value.

Once $f_*^{open}$ and $f_{cor}^{open}$ are given, we set the radial profile of $f_{open}(r)$ as

$$ f_{open}^{(r)} = f_{*}^{open} f_{1}^{\exp}(r) f_{2}^{\exp}(r), \qquad (25a) $$

$$ f_{1}^{\exp}(r) = \min \left[ f_{\text{cor}}^{\text{open}} / f_{*}^{\text{open}}, \exp \left( \frac{r - R_{*}}{2h_{\text{exp}}} \right) \right], \qquad (25b) $$
<table><caption>**Table 1**
Summary of the Input and Output Parameters of Our Simulations</caption><thead><tr><th>P<sub>rot</sub><br>(day)</th><th>f<sub>*</sub><sup>open</sup><br>(10<sup>-3</sup>)</th><th>M&#775;<sub>w</sub><br>(10<sup>-14</sup> M<sub>&odot;</sub> yr<sup>-1</sup>)</th><th>r<sub>&lambda;</sub><br>(R<sub>&odot;</sub>)</th><th>v<sub>r,A</sub><br>(10<sup>2</sup> km s<sup>-1</sup>)</th><th>&tau;<sub>w</sub><br>(10<sup>30</sup> erg)</th></tr></thead><tbody><tr><td>48</td><td>0.466</td><td>1.01</td><td>11.4</td><td>2.14</td><td>0.61</td></tr><tr><td>40</td><td>0.580</td><td>1.26</td><td>12.2</td><td>2.33</td><td>1.03</td></tr><tr><td>32</td><td>0.758</td><td>1.58</td><td>13.3</td><td>2.62</td><td>1.96</td></tr><tr><td>24</td><td>1.07</td><td>2.00</td><td>15.6</td><td>3.08</td><td>4.48</td></tr><tr><td>20</td><td>1.33</td><td>2.23</td><td>17.5</td><td>3.40</td><td>7.54</td></tr><tr><td>16</td><td>1.74</td><td>2.49</td><td>20.5</td><td>3.80</td><td>14.4</td></tr><tr><td>12</td><td>2.46</td><td>2.79</td><td>25.3</td><td>4.41</td><td>33.0</td></tr><tr><td>10</td><td>3.06</td><td>2.87</td><td>29.6</td><td>4.81</td><td>56.0</td></tr><tr><td>8</td><td>4.00</td><td>3.05</td><td>35.8</td><td>5.34</td><td>108</td></tr><tr><td>6</td><td>5.65</td><td>3.25</td><td>46.2</td><td>5.95</td><td>258</td></tr><tr><td>5</td><td>7.03</td><td>3.41</td><td>54.2</td><td>6.57</td><td>433</td></tr><tr><td>4</td><td>9.19</td><td>3.36</td><td>64.9</td><td>7.37</td><td>825</td></tr><tr><td>3</td><td>13.0</td><td>3.28</td><td>80.2</td><td>9.81</td><td>1650</td></tr><tr><td>2</td><td>21.1</td><td>3.09</td><td>107</td><td>16.5</td><td>3890</td></tr></tbody></table>

Note. The first two columns correspond to the input parameters (rotation period and open-flux filling factor), while the last four columns show the output parameters (mass-loss rate, Alfvén radius, Alfvén-point wind velocity, and angular-momentum loss rate).

$$f_2^{\text{exp}}(r) = \frac{\mathcal{F}(r) + f_{\text{cor}}^{\text{open}} + \mathcal{F}(R_*) (f_{\text{cor}}^{\text{open}} - 1)}{f_{\text{cor}}^{\text{open}} (\mathcal{F}(r) + 1)}, \quad (25c)$$

where $\mathcal{F}(r) = \exp\left(\frac{r - r_{\text{exp}}}{\sigma_{\text{exp}}}\right) f_1^{\text{exp}}$ and $f_2^{\text{exp}}$ represent the degree of flux-tube expansion in the chromosphere and corona, respectively. We assume that, in the stellar chromosphere, the flux tube expands so that the plasma beta is fixed until $f_1^{\text{exp}} = f_{\text{cor}}^{\text{open}}/f_{*}^{\text{open}}$ (Tsuneta et al. 2008). Such an expansion is approximately realized by setting the scale height $h_{\text{exp}}$ as

$$h_{\text{exp}} = a_{*}^{2}/g_{*}, \quad (26)$$

where $a_*$ and $g_*$ are the sound speed and gravitational acceleration at the stellar surface, respectively. Here we assume that the pressure scale height of the chromosphere is similar to the photospheric value. For the coronal expansion, we follow the formulation of Kopp & Holzer (1976), with $r_{\text{exp}}/R_* = 1.2$ and $\sigma_{\text{exp}}/R_* = 0.3$.

## 2.4. Alfvén-wave Turbulence

Broadband energy spectra observed in the solar wind indicate that the solar wind is at least partially heated by turbulence (Coleman 1968; Belcher & Davis 1971; Podesta et al. 2007; Chen et al. 2020). In fact, in the outer heliosphere, the observed turbulent dissipation accounts for the required heating rate of the solar wind (Carbone et al. 2009). Although it is still unclear how the solar wind is energized in and below the acceleration region, it is straightforward to assume that the heating process should be similar to what we observe in the distant solar wind; i.e., plasma is heated by turbulence in the solar atmosphere. Alfvén-wave turbulence is a promising candidate of such a heating mechanism. It is a type of MHD turbulence that is driven by the collision of bidirectional Alfvén waves or Elsässer variables (Kraichnan 1965; Dobrowolny et al. 1980; Howes & Nielson 2013). It is likely to develop in the stellar atmosphere (corona) and wind because the reflection of Alfvén waves therein naturally gives rise to wave-wave collisions (Matthaeus et al. 1999; Dmitruk et al. 2002). Alfvén-wave turbulence is now regarded as one of the most

dominant heating processes in coronal holes and the fast solar wind (Cranmer et al. 2007; Verdini & Velli 2007; Perez & Chandran 2013; van Ballegooijen & Asgari-Targhi 2016; Shoda et al. 2019), in coronal loops (van Ballegooijen et al. 2011; Verdini et al. 2012), and in the chromosphere (Verdini & Velli 2007; van Ballegooijen et al. 2011). Note, however, that other processes such as mode conversion (Moriyasu et al. 2004; Suzuki & Inutsuka 2005; Antolin et al. 2008), parametric decay instability (Suzuki & Inutsuka 2006; Tenerani & Velli 2013; Del Zanna et al. 2015; Réville et al. 2018; Shoda et al. 2018b), and phase mixing (Heyvaerts & Priest 1983; Magyar et al. 2017) are also likely to be important.

Without any additional terms, one-dimensional models cannot deal with Alfvén-wave turbulence, because it is a multidimensional process. To model the Alfvén-wave turbulence without expensive numerical cost, phenomenological treatments have been proposed (Hossain et al. 1995; Dmitruk et al. 2002). These models have been validated in previous solar wind simulations (e.g., van Ballegooijen & Asgari-Targhi 2016). Following Shoda et al. (2018a), we introduce a phenomenological model of turbulent dissipation as

$$D_{v}^{\text{turb}} = -\frac{c_{d}}{4\lambda_{\perp}}(|z_{\theta}^{+}| |z_{\theta}^{-}| + |z_{\theta}^{-}| |z_{\theta}^{+}|), \quad (27a)$$

$$D_{b}^{\text{turb}} = -\frac{c_{d}}{4\lambda_{\perp}}(|z_{\theta}^{+}| |z_{\theta}^{-}| - |z_{\theta}^{-}| |z_{\theta}^{+}|), \quad (27b)$$

where $\lambda_{\perp}$ is the perpendicular correlation length and $z_{\theta}^{\pm}$ are Elsässer variables (Elsässer 1950):

$$z_{\theta}^{\pm} = v_{\theta} \mp B_{\theta} / \sqrt{4\pi\rho}. \quad (28)$$

We assume that the correlation length increases with the flux-tube radius:

$$\lambda_{\perp} = \lambda_{\perp,*} \sqrt{\frac{B_{*}}{B_{r}}}. \quad (29)$$

In the photosphere, Alfvénic fluctuations are localized in the intergranular lanes where magnetic flux is concentrated (van Ballegooijen et al. 1998, 2011; Chitta et al. 2012). For this reason, we set the photospheric correlation length of Alfvén-wave turbulence as the typical width of intergranular lanes:

$$\lambda_{\perp,*} = 100 \text{ km}. \quad (30)$$

For the value of $c_d$ in Equations 27(a) and (b), following Shoda et al. (2018a), we set

$$c_d = 0.1, \quad (31)$$

which is supported by both a reduced-MHD simulation (van Ballegooijen & Asgari-Targhi 2017) and a shell-model calculation (Verdini et al. 2019). However, the best choice of $c_d$ remains controversial as one reduced-MHD calculation by Chandran & Perez (2019) shows $c_d \sim 1$. The uncertainty in $c_d$ is not a key issue in this work because the stellar-wind parameters appear to weakly depend on the value of $c_d$ (see Shoda et al. 2018a).

## 2.5. Simulation Domain and Boundary Condition

We solve the basic equations from the photosphere ($r = R_*$) to the outer boundary of the stellar wind ($r = r_{\text{out}}$). The extent of the simulation domain changes depending on $P_{\text{rot}}$, such that the $r_{\text{out}}$ is always beyond the fast-magnetosonic point. For
example, we set $r_{\text{out}}/R_* = 100$ when $P_{\text{rot}} = 24$ day and $r_{\text{out}}/R_* = 690$ when $P_{\text{rot}} = 2$ day.

The spatial resolution of the simulation domain is inhomogeneous. Below $r = 1.02R_*$, the grid size $\Delta r$ is fixed to $\Delta r = 20$ km independent of $\Omega_*$. $\Delta r$ increases with $r$ as a power law of $r$ above $r = 1.02R_*$ until it reaches the maximum value, $\Delta r_{\text{max}}$. To resolve Alfvén waves without large numerical cost, we increase $\Delta r_{\text{max}}$ with rotation rate, $\Omega_*$, because stellar-wind speed and Alfvén velocity are larger in faster rotators. Specifically, $\Delta r_{\text{max}} = 4 \times 10^3$ km for $P_{\text{rot}} = 48$ day and $\Delta r_{\text{max}} = 10^4$ km for $P_{\text{rot}} = 2$ day.

Beyond the outer boundary, $r_{\text{out}}$, a marginal simulation domain is set with gradually increasing grid size. Any numerical errors in the marginal region are unlikely to affect the simulation result since the outer boundary is always beyond the fast-magnetosonic point, where physical fluctuations cannot propagate back into the simulation domain.

Values evaluated at the inner boundary are denoted with the subscript *, and are given as follows. Fixed boundary conditions are imposed for $T$, $v_\phi$, $B_r$, and $B_\phi$:

$$ 
\begin{aligned}
T_* &= \frac{p_*}{c_1 \rho_*} = 6000 \text{ K}, & v_{\phi,*} &= R_* \Omega_* \\
B_{r,*} &= B_{\text{eq},*} = 1300 \text{ G}, & B_{\phi,*} &= 0,
\end{aligned}
\quad (32) 
$$

Note that $B_{\text{eq},*}$ is assumed to be constant with respect to $\Omega_*$. Because the photospheric motion exhibits a much smaller timescale than the rotation, the local property of the photosphere is independent from rotation rate.

To inject MHD waves at the photosphere, we impose time-dependent boundary conditions for density, velocity, and perpendicular magnetic field. Fluctuations of density and radial velocity are given as

$$ \rho_* = \bar{\rho}_* \left( 1 + \frac{v_{r,*}}{a_*} \right) \quad (33) $$

where

$$ \bar{\rho}_* = 10^{-7} \text{ g cm}^{-3}, \quad a_* = \sqrt{c_1 T_*}. \quad (34) $$

The time-dependent radial (vertical) velocity $v_{r,*}$ has a broadband spectrum of

$$ v_{r,*} \propto \sum_{N=0}^{10} \sin(2\pi f_N^l t + \phi_N^l) / \sqrt{2\pi f_N^l}, \quad (35) $$

where the $\phi_N^l$ is a random phase and the wave frequency $f_N^l$ ranges in $3.33 \times 10^{-3}$ Hz $\le f_N^l \le 3.33 \times 10^{-2}$ Hz. The lower limit of $f_N^l$ is set to be the cutoff frequency of acoustic waves at the stellar surface. The amplitude of $v_{r,*}$ is set so that the rms amplitude of upward acoustic waves is 0.9 km s$^{-1}$. Considering the downward wave contribution, the rms velocity at the surface is approximately 1.3 km s$^{-1}$, consistent with solar observations (Oba et al. 2017; Ishikawa et al. 2020).

The $\theta$-component of the velocity and magnetic field is given in terms of Elsässer variables from Equation (28). We impose a zero-derivative boundary condition on $z_\theta^-$ such that reflected Alfvén waves can be absorbed through the bottom boundary:

$$ \left. \frac{\partial}{\partial r} z_{\theta}^{-} \right|_{*} = 0. \quad (36) $$

Like $v_{r,*}$, the upward Elsässer variable $z_{\theta,*}^+$ is given with a broadband spectrum as follows:

$$ z_{\theta,*}^{+} \propto \sum_{N=0}^{20} \frac{\sin(2\pi f_N^l t + \phi_N^l)}{\sqrt{2\pi f_N^l}}, \quad (37) $$

where $\phi_N^l$ is a random phase and $f_N^l$ ranges in $1.00 \times 10^{-3}$ Hz $\le f_N^l \le 1.00 \times 10^{-2}$ Hz. The lower and upper limits of this frequency range approximate the turn-over timescale in granules and intergranular lanes (Hirzberger et al. 1999). The amplitude is tuned so that the rms amplitude of upward Alfvén waves is 1.2 km s$^{-1}$, which yields the rms photospheric transverse velocity of 1.7 km s$^{-1}$. This value is consistent with observations of the solar surface convection (de Wijn et al. 2008; Oba et al. 2020). There is evidence that the imposed spectrum of transverse waves may affect the dynamics of the resulting stellar wind (Shoda et al. 2018b). Therefore, in the future, simulations should be performed with a self-consistent convection zone (e.g., Rempel 2017) to remove the uncertainty in the wave generation process.

### 3. Trends in the Wind Simulations

#### 3.1. Overview of Rotation Dependence

In Figure 2, we show the time-averaged radial profiles of the simulated winds with various rotation rates: $P_{\text{rot}} = 24$ days (red solid line), 12 days (orange dashed-dotted line), 6 days (green dashed line), and 3 days (blue dotted line). To eliminate the influence of initial conditions, time averaging is conducted after the system reaches a quasi-steady state that is independent from the choice of initial condition. The averaging time is typically $1.2 \times 10^5$ s $\approx 1.39$ day. Note that we show the rms value for the wave amplitude $v_\phi$.

Panel (e) directly reflects the different rotation velocities used in our simulations; $v_\phi$ increases with rotation rate. Up to a certain height, the wind corotates with the stellar surface. The corotation breaks up below the Alfvén radius, and $v_\phi$ in turn begins to decrease. This behavior is consistent with the Weber-Davis solution, $v_{\phi,\text{WD}}$, which predicts

$$ 
\begin{aligned}
v_{\phi,\text{WD}} &\approx r \Omega_* && (r \ll r_A), \\
v_{\phi,\text{WD}} &\approx r_A^2 \Omega_* / r && (r \gg r_A).
\end{aligned}
\quad (38) 
$$

In fact, the radial profile of $v_\phi$ almost perfectly coincides with the Weber-Davis solution (thin black lines). Similarly, panel (f) shows that $-B_\phi/B_r$ increases with $r$, typically $-B_\phi/B_r \propto r$ far away from the star. This is also consistent with Weber-Davis solution that predicts the Parker-spiral relation (Parker 1958).

$$ -B_{\phi,\text{WD}}/B_r \approx r \Omega_* / v_{r,\infty} \quad (r \gg r_A). \quad (39) $$

Diamonds on each line indicate the Alfvén point. The Alfvén point occurs at larger radii as the rotation rate increases. This is a natural consequence of larger open magnetic flux and the larger coronal Alfvén speed of the faster rotators. As shown in Table 1, the largest Alfvén radius exceeds the mean orbital radius of Mercury ($\sim 83R_\odot$). This indicates that Mercury was possibly subject to magnetic star-planet interactions with the young Sun (see, e.g., Strugarek et al. 2014; Folsom et al. 2020).

Panels (a) and (b) show the weak dependence of density and temperature on $\Omega_*$. We can hardly see any differences between the four lines, especially near the coronal base. The
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison of the quasi-steady-state solutions after time averaging. The four lines correspond to $P_{\text{rot}} = 24$ days (red solid line), $P_{\text{rot}} = 12$ days (orange dashed-dotted line), $P_{\text{rot}} = 6$ days (green dashed line), and $P_{\text{rot}} = 3$ days (blue dotted line), respectively. Panel (a): mass density $\rho$, panel (b): temperature $T$, panel (c): radial velocity $v_r$, panel (d): rms wave amplitude $v_{\theta, \text{rms}}$, panel (e): rotation velocity $v_{\phi}$, and panel (f): field inclination $-B_{\phi}/B_r$. The thin black lines in panel (e) indicate the Weber-Davis solution and, for better visualization, the wind simulations are shown with dotted lines (only in panel (e)). The diamonds indicate the Alfvén point.</figcaption></figure>

rotation-rate dependences of stellar-wind parameters are shown more clearly in Figure 3. Panel (a) shows the dependence of the electron number density at the coronal base (where $p = 0.03$ dyne cm$^{-2}$), $n_{e,\text{cor}}$, on $\Omega_*$. In contrast to the observed relation (dashed line; Ivanova & Taam 2003), our simulations do not show any correlation. In panel (b), we show the coronal-base temperature ($T_{\text{cor}}$, circles) and the maximum temperature ($T_{\text{max}}$, diamonds) against $\Omega_*$, both of which are inconsistent with the relation $T_{\text{cor}} \propto \Omega_*^{0.45}$ found by O’Fionnagáin & Vidotto (2018).

These inconsistencies can be explained in the framework of our model. The observed scaling relations of $n_e$ and $T$ are deduced from X-ray emission that mostly comes from closed magnetic field loops, while the wind comes from the open-field region (e.g., Cranmer 2009). Therefore, the difference between our result and the observed trends indicates that the density and temperature scale in different ways for open and closed regions. To exemplify this, we also show a power-law fitting between $T_{\text{max}}$ and $\Omega_*$: $T_{\text{max}} \propto \Omega_*^{0.18}$, which is closer to the relation by O’Fionnagáin & Vidotto (2018).

The constant coronal temperature with respect to $\Omega_*$ is due to the constant Alfvén-wave energy flux transmitted into the corona. As Alfvén waves are the only source of coronal heating in our model, their constant energy flux leads to a constant coronal temperature. The coronal density is in general determined by the energy balance between radiative cooling and conductive heating (Hammer 1982; Withbroe 1988). Because the constant temperature yields constant basal conductive heating, the coronal density is also kept fixed with $\Omega_*$. Therefore, the constant coronal temperature and density are attributed to constant Alfvén-wave

energy flux in the corona, which is discussed in more detail in Section 4.

Panel (c) shows the $\Omega_*$ dependence of the wind terminal velocity. For $\Omega_*/\Omega_{\odot} \lesssim 7$, the wind velocity weakly depends on the rotation rate. However, when $\Omega_*/\Omega_{\odot} \gtrsim 7$, the wind velocity drastically increases with rotation. Beyond a critical rotation rate, the magneto-centrifugal force dominates the force balance in the wind acceleration, resulting in a strong acceleration of the wind (Belcher & MacGregor 1976). The critical point $\Omega_*/\Omega_{\odot} \approx 7$ turns out to be the regime-changing point in terms of energy budget. We will discuss this further in Section 4.

### 3.2. Angular-momentum Loss Rate (Torque)

One of the principal purposes of this work is to investigate whether the Alfvén-wave-driven magnetic rotator wind model can explain the observed spin-down of low-mass stars. Ignoring the core-envelope decoupling and internal-structure evolution, the stellar rotational evolution is described as

$$I_* \frac{d\Omega_*}{dt} = -\tau_w, \qquad (40)$$

where $I_*$ is the momentum of inertia of the star. If one assumes that the wind torque is approximated by $\tau_w \propto \Omega_*^{p+1}$, the solution of the rotational evolution yields

$$\Omega_* \propto t^{-1/p}, \qquad (41)$$

from which the Skumanich relation is reproduced when $p = 2$.

In the quasi-steady state, an analytical formulation of torque can be obtained (see, e.g., Lamers & Cassinelli 1999). The time-averaged mass conservation and magnetic-flux conservation are
<figure><img src="image_6.png" /><figcaption>**Figure 3.** $\Omega_*$ dependences of stellar-wind parameters. Panel (a): coronal-base electron number density, $n_{e,cor} \propto \Omega_*^{0.6}$ (Ivanova & Taam 2003). Panel (b): coronal-base temperature ($T_{cor}$, circles) and maximum temperature ($T_{max}$, diamonds). The dashed line is the observational single-power-law relation (O’Fionnagáin & Vidotto 2018). The dotted line is the power-law fitting between $T_{max}$ and $\Omega_*$. Panel (c): termination velocity of stellar wind.</figcaption></figure>

given by

$$ \dot{M}_w = 4\pi r^2 f^{\text{open}} \rho v_r = \text{const.}, \qquad (42) $$

$$ \Phi_{\text{open}} = 4\pi r^2 f^{\text{open}} B_r = \text{const.} \qquad (43) $$

Combining these equations gives the following identity:

$$ \frac{\Phi_{\text{open}}^2}{16\pi^2 \dot{M}_w} = r^2 f_{\text{open}}^2 \frac{v_A^2}{v_r} = \text{const.} = r_A^2 v_{r,A}, \qquad (44) $$

where we assume that the open-flux filling factor at the Alfvén point is unity. For simplicity, we assume that $r_A$ is spherically

symmetric. The torque is then given by

$$ \tau_w = \frac{2}{3} \dot{M}_w r_A^2 \Omega_* = \frac{2}{3} \frac{(B_{r,*} f_*^{\text{open}})^2}{v_{r,A}} R_*^4 \Omega_* \qquad (45) $$

where $v_{r,A}$ is the wind velocity at the Alfvén point. Substituting $B_{r,*} = 1300$ G, $R_* = 6.96 \times 10^{10}$ cm, and Equation (22), a semi-analytical expression of $\tau_w$ is obtained:

$$ \tau_w = 1.22 \times 10^{30} \left( \frac{v_{r,A}}{v_{g,\odot}} \right)^{-1} \left( \frac{\Omega_*}{\Omega_{\odot}} \right)^{3.4} \text{erg}, \qquad (46) $$

where $v_{g,\odot} = \sqrt{2GM_{\odot}/R_{\odot}} \approx 617$ km s$^{-1}$ is the escape velocity at the solar surface. To express $v_{r,A}/v_{g,\odot}$ as a function of $\Omega_*$, a numerical simulation is required.

The torques calculated from our numerical simulations are shown in Figure 4. Panel (a) shows how the wind torque varies with the stellar rotation rate. In the whole range of $\Omega_*$, the simulation results are well fitted by a single power law of

$$ \tau_w = 2.59 \times 10^{30} \left( \frac{\Omega_*}{\Omega_{\odot}} \right)^{2.82} \text{erg}, \qquad (47) $$

which yields

$$ \Omega_* \propto t^{-0.549}. \qquad (48) $$

We note that the data points in the fast-rotator regime deviate slightly from the fit line, which is possibly a result of the regime change (see Section 4). This spin-down law is consistent with the recent gyrochronology relation from Angus et al. (2015) who find $P_{\text{rot}} \propto t^{0.55}$. For $\Omega/\Omega_{\odot} = 1$, the calculated angular-momentum loss rate ($\tau_{w,\odot} = 2.59 \times 10^{30}$ erg) matches with the observed solar wind torque (Finley et al. 2018, 2019). However, this value is still smaller than the stellar-observation-based empirical value ($\tau_{w,\odot} = 6.3 \times 10^{30}$ erg; Matt et al. 2015) by a factor of 2.4. One possibility for this gap is that the Sun has a smaller amount of open magnetic flux than typical sunlike stars. We will discuss this point in more detail in Section 5.

Comparing Equations (46) and (47), one can tell that $v_{r,A}/v_{g,\odot}$ should depend on $\Omega_*$, specifically $v_{r,A}/v_{g,\odot} \propto \Omega_*^{0.58}$. This is directly confirmed in panel (b). If we simply assume that $v_{r,A}/v_{g,\odot} \approx 1$, which is not a bad approximation, the torque scales as $\tau_w \propto \Omega^{3.4}$, yielding slower spin-down than the observed one: $\Omega_* \propto t^{-0.417}$. In this respect, the $\Omega_*$ dependence of $v_{r,A}/v_{g,\odot}$ is also important in evaluating the spin-down law.

### 3.3. Alfvén Radius

Several works have produced semi-analytic scaling relations for $r_A$ (Kawaler 1988; Matt & Pudritz 2008; Matt et al. 2012). Here we compare two relations that are based on the open magnetic flux in the wind. These scaling relations are given in terms of the dimensionless wind-magnetization parameter $\Upsilon_{\text{open}}$, which is:

$$ \Upsilon_{\text{open}} = \frac{\Phi_{\text{open}}^2}{R_*^2 \dot{M}_w v_{g,\odot}}. \qquad (49) $$

For comparison, we convert the equatorial Alfvén radius from our simulations, $r_A$, to a latitudinally averaged value, $\langle r_A \rangle$,
<figure><img src="image_7.png" /><figcaption>Figure 4. Circles in each panel show the simulated panel (a): angular-momentum loss rates (torques) of the winds, $\tau_w$, and panel (b): stellar-wind velocities at the Alfvén point, $v_{r,A}$, normalized by the escape velocity at the surface $v_{g,\odot}$. Also shown by the dashed lines are power-law fittings to circles.</figcaption></figure>

based on the following formulation (Washimi & Shibata 1993):

$$ \tau_w = \frac{2}{3} \dot{M}_w r_A^2 \Omega_* = \dot{M}_w \langle r_A \rangle^2 \Omega_* \quad (50) $$

Therefore, $\langle r_A \rangle = \sqrt{2/3} r_A$, where spherical symmetry has been assumed.

We compare first to the scaling law given by Finley & Matt (2018):

$$ \langle r_A \rangle_{\text{FM18}} / R_* = \begin{cases} 0.33 \ \Upsilon_{\text{open}}^{0.371} & (\text{dipole}), \\ 0.46 \ \Upsilon_{\text{open}}^{0.329} & (\text{general}), \end{cases} \quad (51) $$

where the first case is fitted from simulations with only dipole fields and the second case corresponds to a fit using a range of simulations with combinations of dipole, quadrupole, and octupole geometries.

The second scaling relation is given by Réville et al. (2015) as follows:

$$ \langle r_A \rangle_{\text{R15}} / R_* = 0.64 \left[ \frac{\Upsilon_{\text{open}}}{\sqrt{1 + (f_{\text{bu}}/0.06)^2}} \right]^{0.31}, \quad (52) $$

where $f_{\text{bu}} = \Omega_* R_*^{3/2} (GM_*)^{-1/2}$ is the break-up fraction of the rotation speed.

In Figure 5 we compare our results with Finley & Matt (2018; top panel, dashed and dashed-dotted lines) and Réville et al. (2015; bottom panel, dashed line). In each panel, our results and a fitted power law are indicated by circles and the solid line. The results of both Réville et al. (2015) and Finley & Matt (2018) are consistent with our results, indicating that scaling relations for the Alfvén radius are robust regardless of simulation setting. Note that the open flux is an output of the simulations in Réville et al. (2015) and Finley & Matt (2018) while the mass-loss rate is mostly controlled by the coronal density and temperature imposed at the boundary condition. On the other hand, in our calculations, the mass-loss rate is an output while the open flux is an input. A more self-consistent treatment requires a full-sphere simulation with physics-based coronal heating and chromospheric evaporation.

As shown in Figure 5, our simulations yield a power-law relation between $\langle r_A \rangle / R_*$ and $\Upsilon_{\text{open}}$ of

$$ \langle r_A \rangle / R_* \propto \Upsilon_{\text{open}}^{0.36} \quad (53) $$

The origin of the exponent 0.36 is explained as follows. Rewriting Equation (44) in terms of $\Upsilon_{\text{open}}$, one obtains

$$ \langle r_A \rangle^2 / R_*^2 \propto \Upsilon_{\text{open}} \left( \frac{v_{r,A}}{v_{g,\odot}} \right)^{-1} \quad (54) $$

Suppose a power-law relation $v_{r,A}/v_{g,\odot} \propto (r_A/R_*)^q$ is satisfied, then

$$ \langle r_A \rangle / R_* \propto \Upsilon_{\text{open}}^{1/(2+q)} \quad (55) $$

It is evident from Figure 6 that the power-law relation $v_{r,A}/v_{g,\odot} \propto (r_A/R_*)^q$ is satisfied with $q = 0.77$, which yields the exponent in Equation (55) of $1/(2+q) = 0.361$. Also, our results are consistent with the scaling law by Pantolmos & Matt (2017) if we adopt the sound-to-escape velocity ratio $a/v_{g,\odot} = 0.23$. According to Pantolmos & Matt (2017), the $q$ value is sensitive to the coronal temperature. In our model, the coronal temperature is almost constant with respect to $\Omega_*$, and thus, all of our simulations are fitted by a unique $q$ value.

### 3.4. Mass-loss Rate

In this section, we discuss another interesting topic: the rotation dependence of mass-loss rates, $\dot{M}_w$. The top panel of Figure 7 shows how $\dot{M}_w$ varies with the stellar rotation rate (circles). Also shown by diamonds are the results with a fixed open-flux filling factor $f_{*}^{\text{open}} = 10^{-3}$ (see also Figure 10). $\dot{M}_w$ increases with $\Omega_*$ in the slow-rotator regime and saturates around $\dot{M}_w \sim 3.4 \times 10^{-14} M_{\odot} \text{ yr}^{-1}$ in the faster rotation cases. Observations of astrospheric line absorption show that the
<figure><img src="image_8.png" /><figcaption>Figure 5. Comparison of Alfvén-radius scaling laws from Finley & Matt (2018; top panel) and Réville et al. (2015; bottom panel).</figcaption></figure>

mass-loss rate tends to increase with X-ray flux (Wood et al. 2002, 2005, 2014), and thus with rotation rate (Güdel et al. 1997; Ribas et al. 2005; Wright et al. 2011; Magauda et al. 2020). However, we need to note that what is actually obtained by the astrospheric observation is the characteristic ram pressure, $P_w = 4\pi r_{\text{out}}^2 \rho_{\text{out}} v_{\text{r,out}}^2$ (Holzwarth & Jardine 2007), not the mass-loss rate, $\dot{M}_w = 4\pi r_{\text{out}}^2 \rho_{\text{out}} v_{\text{r,out}}$. Bearing this in mind, we henceforth focus on $P_w$ for comparison with observation.

The bottom panel of Figure 7 shows the relation between $P_w/P_{w,\odot}$ and $\Omega_*/\Omega_\odot$, which has a power-law relation of

$$P_w/P_{w,\odot} = (\Omega_*/\Omega_\odot)^{0.83}, \qquad (56)$$

<figure><img src="image_9.png" /><figcaption>Figure 6. $\langle r_A \rangle/R_*$ vs. $v_{r,A}/v_{g,\odot}$ (circles) and a power-law fit ($v_{r,A}/v_{g,\odot} \propto (\langle r_A \rangle/R_*)^{0.77}$; solid line). Also shown by the dashed line is the relation given by Pantolmos & Matt (2017) with $a/v_{g,\odot} = 0.23$.</figcaption></figure>

where we set $P_{w,\odot} = 5.0 \times 10^{19}$ dyne. Following Wright et al. (2011), we convert the rotation rate $\Omega_*$ to the X-ray flux $F_X$ as

$$F_X/F_{X,\odot} = (\Omega_*/\Omega_\odot)^{2.18}, \qquad (57)$$

where $F_{X,\odot} = 3 \times 10^4 \text{ erg cm}^{-2} \text{ s}^{-1}$. Note that (1) all of our simulation runs are in the unsaturated regime (in which stellar activities correlate with stellar rotation) and (2) there is a one-to-one relation between $F_X$ and $L_X/L_{\text{bol}}$ because the stellar radius and luminosity are fixed in our simulations. Combining Equations (56) and (57),

$$P_w/P_{w,\odot} = (F_X/F_{X,\odot})^{0.38}. \qquad (58)$$

Figure 8 shows $F_X-P_w$ for the astrospheric observations taken from Wood et al. (2014; symbols), the empirical relation from Wood et al. (2005; blue solid line), and our result (Equation (58); red dashed line). Our simulation result is consistent (within a factor of three) with the observations of 61 Vir, Sun, $\alpha$ Cen, $\varepsilon$ Ind, 61 Cyg A, $\xi$ Boo, Prox Cen, and EV Lac. A similar trend is found in the work of Holzwarth & Jardine (2007). Our model is able to explain a good fraction of the observations, although there exist non-negligible offsets for three K-dwarfs (36 Oph, 70 Oph, and $\varepsilon$ Eri). It is left for future work to test whether Equations (56) and (58) are valid for non-sunlike stars.

## 4. Wind Energetics

The physics of the stellar-wind heating and acceleration can be inferred by following the energy flow from the stellar surface to interplanetary space. For example, one can estimate the stellar-wind mass-loss rate analytically, based on wind energetics (Hansteen & Leer 1995; Cranmer & Saar 2011; Suzuki 2018). To understand what causes the saturation of the mass-loss rate, the energy budget in the stellar wind is discussed.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (Top panel) Mass-loss rate, $\dot{M}_w$, vs. rotation rate, $\Omega_*$: The circles show the fiducial cases ($f_*^{\text{open}} = 10^{-3}(\Omega_*/\Omega_\odot)^{1.2}$) and diamonds show the results with fixed $f_*^{\text{open}}$ ($f_*^{\text{open}} = 10^{-3}$; see also Figure 10 and discussions there). (Bottom panel) Characteristic wind ram pressure, $P_w = 4\pi r^2 f \rho v_\perp^2$, normalized by the solar value vs. rotation rate, $\Omega_*$: Shown by the dashed line is a power-law fit to the numerical results: $P_w/P_{w,\odot} = (\Omega_*/\Omega_\odot)^{0.83}$.</figcaption></figure>

### 4.1. Energy Conservation

After time averaging, the energy conservation law is written as follows:

$$ \frac{d}{dr}(L_K + L_E + L_A - L_C - L_G) = -4\pi r^2 f^{\text{open}} Q_{\text{rad}}, \quad (59) $$

where

$$ L_K = \frac{1}{2} \rho v_r^3 4\pi r^2 f^{\text{open}}, \quad (60a) $$

<figure><img src="image_11.png" /><figcaption>**Figure 8.** X-ray flux $F_X$ vs. characteristic ram pressure normalized by solar value $P_w/F_{w,\odot}$. The symbols indicate the astrophysical observations; diamonds, circles, and stars stand for G-type, K-type, and M-type stars, respectively, and binaries are doubly marked with corresponding symbols. Also, the red dashed and blue solid lines show our simulation result and the empirical relation proposed by Wood et al. (2005), respectively.</figcaption></figure>

$$ L_E = \frac{\gamma}{\gamma - 1} \rho v_r 4\pi r^2 f^{\text{open}}, \quad (60b) $$

$$ L_A = \left[ \left( \frac{1}{2} \rho v_{\perp}^2 + \frac{B_{\perp}^2}{4\pi} \right) v_r - \frac{B_r}{4\pi} (v_{\perp} \cdot B_{\perp}) \right] 4\pi r^2 f^{\text{open}}, \quad (60c) $$

$$ L_C = -F_C 4\pi r^2 f^{\text{open}}, \quad (60d) $$

$$ L_G = \rho v_r \frac{GM_*}{r} 4\pi r^2 f^{\text{open}} = \dot{M}_w \frac{GM_*}{r}. \quad (60e) $$

$L_A$, $L_E$, $L_A$, $L_C$, and $L_G$ correspond to the wind kinetic energy flux, enthalpy flux, Alfvén-wave energy flux, conductive flux, and gravitational energy flux, respectively. Bearing in mind that the $\theta$ and $\phi$ components stand for Alfvén waves and rotation, respectively, we can decompose $L_A$ as

$$ L_A = L_A^{\text{wav}} + L_A^{\text{rot}}, \quad (61a) $$

$$ L_A^{\text{wav}} = \left[ \left( \frac{1}{2} \rho v_{\theta}^2 + \frac{B_{\theta}^2}{4\pi} \right) v_r - \frac{B_r}{4\pi} v_{\theta} B_{\theta} \right] 4\pi r^2 f^{\text{open}}, \quad (61b) $$

$$ L_A^{\text{rot}} = \left[ \left( \frac{1}{2} \rho v_{\phi}^2 + \frac{B_{\phi}^2}{4\pi} \right) v_r - \frac{B_r}{4\pi} v_{\phi} B_{\phi} \right] 4\pi r^2 f^{\text{open}}, \quad (61c) $$

where $L_A^{\text{wav}}$ and $L_A^{\text{rot}}$ correspond to the luminosities of Alfvén waves and magneto-rotation, respectively.

### 4.2. Energetics in the Corona and Above

The energy conservation is more simply approximated above the coronal base, where the enthalpy flux and the radiative loss
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Luminosities in Equation (63) as functions of rotation rate $\Omega_*$ Shown are the wind kinetic energy flux at the outer boundary $L_{K,out}$ (red dashed line), rotational energy flux at the outer boundary $L_{A,out}^{rot}$ (blue dashed line), gravitational energy flux at the coronal base $L_{G-cor}$ (orange solid line), downward conductive energy flux at the coronal base $L_{C-cor}$ (green solid line), Alfvén-wave energy flux at the coronal base $L_{A-cor}^{wav}$ (cyan solid line), and rotational energy flux at the coronal base $L_{A-cor}^{rot}$ (blue solid line), respectively.</figcaption></figure>

are negligibly small.

$$ \frac{d}{dr}(L_K + L_A^{\text{wav}} + L_A^{\text{rot}} - L_C - L_G) \approx 0. \quad (62) $$

Integrating this equation from the coronal base to the distant stellar wind and ignoring minor components, one obtains

$$ L_{A-cor}^{\text{wav}} + L_{A-cor}^{\text{rot}} - L_{C-cor} - L_{G-cor} \approx L_{K,out} + L_{A-out}^{\text{rot}} \quad (63) $$

where $X_{cor}$ and $X_{out}$ denote $X$ measured at the coronal base and outer boundary, respectively.

Taking the coronal base as $r = 1.02R_*$, we have confirmed that this approximated energy balance relation is satisfied to within 2% error. In Figure 9, we plot each term in Equation (63) as a function of rotation rate $\Omega_*$. Figure 9 has several features:

1. In the slow-rotator regime ($\Omega_*/\Omega_{\odot} \lesssim 4$), the dominant coronal energy injection is by Alfvén waves: $L_{A-cor}^{\text{wav}} \gg L_{A-cor}^{\text{rot}}$. In this regime, the energy balance is approximated as $L_{A-cor}^{\text{wav}} \approx L_{G-cor} + L_{K,out}$, as assumed by Cranmer & Saar (2011).

2. In the fast-rotator regime ($\Omega_*/\Omega_{\odot} \gtrsim 10$), the rotation components becomes dominant: $L_{A-cor}^{\text{wav}} \ll L_{A-cor}^{\text{rot}}$. The energy balance relation is then $L_{A-cor}^{\text{rot}} \approx L_{K,out} + L_{A-out}^{\text{rot}}$.

3. The “regime change” from wave-driven wind ($L_{A-cor}^{\text{wav}} > L_{A-cor}^{\text{rot}}$) to rotation-driven wind ($L_{A-cor}^{\text{wav}} < L_{A-cor}^{\text{rot}}$) takes place around $\Omega_*/\Omega_{\odot} \approx 7$, or equivalently $P_{\text{rot}} \approx 3.6$ day. Note that the regime-changing period strongly depends on the filling factor of open flux.

An interesting behavior of the rotation-driven wind is that, in spite of the rapid increase of rotational energy injection $L_{A-cor}^{\text{rot}}$ with $\Omega_*$, the $L_{G-cor} \propto \dot{M}_w$ does not increase. This is because the

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Energy-flux dependence on the rotation rate $\Omega_*$ of the open-flux filling factor fixed to $f_*^{\text{open}} = 10^{-3}$. Definitions of the lines and markers are the same as in Figure 9.</figcaption></figure>

wind density is determined by the energy injected below the sonic (slow-magnetosonic) point (Hammer 1982; Leer et al. 1982; Hansteen & Leer 1995; Hansteen & Velli 2012). The magneto-rotational acceleration (magneto-centrifugal force) works in the supersonic region and works to accelerate the stellar wind without increasing the mass-loss rate. To show this, we perform a set of test simulations with different rotation rates ($P_{\text{rot}} = 24, 8, 2.4, 0.8$ day) and fixed open-flux filling factor to $f_*^{\text{open}} = 10^{-3}$ to see the purely rotational effect on the wind. Figure 10 shows $L_{K,out}$ (black), $L_{G-cor} \propto \dot{M}_w$ (blue), and $L_{A-cor}^{\text{rot}}$ (red) as functions of $\Omega_*$. The gravitational luminosity $L_{G-cor}$ (or the mass-loss rate) does not respond to the increase of rotational energy flux. Instead, the terminal wind kinetic energy flux increases with rotational energy. For this reason, enhanced rotation rate does not lead to enhanced mass-loss rate.

An approximated relation for the mass-loss rate is derived from our analysis. Consider the energy balance assumed in Cranmer & Saar (2011) of

$$ L_{K,out} + L_{G-cor} \approx L_{A-cor}^{\text{wav}} \quad (64) $$

We have already shown that this holds for the slow-rotator regime but breaks down in the presence of large rotational energy injection. However, Figure 9 shows that, in the whole range of $\Omega_*$, $L_{A-cor}^{\text{wav}}$ is well approximated as

$$ L_{A-cor}^{\text{wav}} \approx 2L_{G-cor} \quad (65) $$

which yields

$$ \dot{M}_w \approx \frac{L_{A-cor}^{\text{wav}}}{v_{g,\odot}^2} \quad (66) $$

Equation (66) is validated as follows. As we have already shown, the magneto-centrifugal force works to enhance the wind velocity but not to increase the mass-loss rate. In other words, the mass-loss rate remains unchanged even if the magneto-centrifugal force does not work. Therefore, the rotational terms in the energy budget in Equation (63) can be ignored in discussing the mass-loss rate. Ignoring the
<figure><img src="image_14.png" /><figcaption>**Figure 11.** $\Omega_*$ dependences of Poynting-flux luminosities. Alfvén-wave luminosities are shown by the solid lines for photospheric (red) and coronal-base (blue) values. Rotational luminosities are shown by the dashed lines for photospheric (red) and coronal-base (blue) values. A semi-analytical relation $L_{A,cor}^{\text{rot}} \propto \Omega_*^{3.83}$ is indicated by the dashed line.</figcaption></figure>

conductive flux, which is always minor $L_{C,cor}$, the energy conservation is reduced to

$$L_{K,out} + L_{G,cor} = \frac{1}{2} \dot{M}_w (v_{r,out}^2 + v_{g,\odot}^2) \approx L_{A,cor}^{\text{wav}} \quad (67)$$

We have confirmed by numerical simulations (not shown here) that, in the absence of magneto-centrifugal force, $v_{r,out}$ is always approximated by $v_{g,\odot}$. With $v_{r,out} \approx v_{g,\odot}$, Equation (66) is derived from Equation (67). Now, it turns out that the saturation of mass-loss rate comes from the saturation of $L_{A,cor}^{\text{wav}}$ (energy flux of Alfvén waves transmitted into the corona), which is further discussed in the following section.

### 4.3. Alfvén-wave Energetics in the Chromosphere

We further investigate the rotation dependence of Alfvén-wave energy flux (or luminosity) at the coronal base $L_{A,cor}^{\text{wav}}$. The solid lines in Figure 11 show the $\Omega_*$ dependences of the wave luminosities measured at the photosphere ($L_{A,pho}^{\text{wav}}$, red line), and the coronal base ($L_{A,cor}^{\text{wav}}$, blue line). The dashed lines represent the rotational luminosities: $L_{A,pho}^{\text{rot}}$ (red) and $L_{A,cor}^{\text{rot}}$ (blue). The photospheric value is measured 20 km above the stellar surface, eliminating the direct influence of the lower boundary condition. We can tell several interesting features from Figure 11:

1. Near the stellar surface, the wave energy flux is always larger than the rotation energy flux in the parameter range of our simulations. Our fastest rotating case has a surface rotation velocity of ~25 km s⁻¹, which is much larger than the wave amplitude of ~1.2 km s⁻¹. However, because the azimuthal magnetic field $B_\phi$ is small in the lower atmosphere (as the Weber–Davis solution predicts), the energy flux of magneto-rotation remains smaller than the wave energy flux.

2. $L_{A,cor}^{\text{wav}}$ is much smaller than $L_{A,pho}^{\text{wav}}$. This means that a large fraction of Alfvén waves dissipates between the photosphere and coronal base. Interestingly, the “energy transmission rate,” $L_{A,cor}^{\text{wav}}/L_{A,pho}^{\text{wav}}$, decreases with $\Omega_*$: This is why $L_{A,cor}^{\text{wav}}$ saturates with respect to $\Omega_*$ at rapid rotation in spite of the power-law dependence of $L_{A,pho}^{\text{wav}}$ on $\Omega_*$.

3. In contrast to the wave luminosity, the rotational luminosity hardly decreases between the photosphere and coronal base. The rotational luminosity steeply increases with $\Omega_*$, which is consistent with semi-analytical predictions ($L_{A,cor}^{\text{rot}} \propto \Omega_*^{3.83}$; see Appendix B). As a result, at the coronal base, the rotational energy flux overtakes the wave energy flux in the fast-rotator regime, typically from $P_{\text{rot}} \lesssim 4$ day. This transition is responsible for the wind regime change discussed in Section 4.2

To summarize, the stellar wind experiences a regime change at $P_{\text{rot}} \approx 4$ day. This results from a significant decrease in coronal Alfvén-wave energy flux, which is overtaken by the rotational energy flux at $P_{\text{rot}} \approx 4$ day. Note that the Alfvén-wave energy flux at the stellar surface is always larger than the rotation energy flux.

To determine the reason for the saturation of Alfvén-wave luminosity, we analyze the Alfvén-wave energy loss in the chromosphere. The conservation of Alfvén-wave energy flux (luminosity) is given as follows:

$$\frac{\partial}{\partial t}\left(\frac{1}{2} \rho v_{\theta}^{2}+\frac{B_{\theta}^{2}}{8 \pi}\right)+\frac{1}{4 \pi r^{2} f^{\text {open }}} \frac{\partial}{\partial r} L_{A}^{\text {wav }}=-\varepsilon_{r \leftrightarrow \theta}-Q_{\text {turb }} \quad (68)$$

where $Q_{\text{turb}}$ is the turbulent dissipation and $\varepsilon_{r \leftrightarrow \theta}$ represents the energy conversion between Alfvén wave and longitudinal motion:

$$Q_{\text{turb}}=c_{d} \rho \frac{\left|z_{\theta}^{+}\right| z_{\theta}^{- 2}+\left|z_{\theta}^{-}\right| z_{\theta}^{+ 2}}{4 \lambda_{\perp}} \quad (69)$$

$$\varepsilon_{r \leftrightarrow \theta}=-v_{r} \frac{\partial}{\partial r}\left(\frac{B_{\theta}^{2}}{8 \pi}\right)+v_{r}\left(\rho v_{\theta}^{2}-\frac{B_{\theta}^{2}}{4 \pi}\right) \frac{d}{d r} \ln (r f^{\text {open }}). \quad (70)$$

In this paper, we shall interpret $\varepsilon_{r \leftrightarrow \theta}$ as mode conversion, since the energy conversion between transverse and longitudinal waves is described by this term. Note that the mode conversion works efficiently in the chromosphere (Rosenthal et al. 2002; Bogdan et al. 2003), both transforming longitudinal waves to transverse waves (Schunker & Cally 2006; Shoda & Yokoyama 2018) and transverse waves to longitudinal waves (Hollweg et al. 1982; Kudoh & Shibata 1999; Matsumoto & Shibata 2010). We can define the loss fractions as

$$\Delta L_{A,\text{tot}}^{\text{wav}}=\Delta L_{A,\text{turb}}+\Delta L_{A, r \leftrightarrow \theta}, \quad (71)$$

$$\Delta L_{A, \text {turb }}=\int_{R_{*}}^{r} d r \ 4 \pi r^{2} f^{\text {open }} Q_{\text {turb }} \quad (72)$$

$$\Delta L_{A, r \leftrightarrow \theta}=\int_{R_{*}}^{r} d r \ 4 \pi r^{2} f^{\text {open }} \varepsilon_{r \leftrightarrow \theta}. \quad (73)$$

In Figure 12, we show the Alfvén-wave energy loss in the chromosphere. In panel (a), we show $L_{A}^{\text{wav}}$ as a function of height for different rotation periods. In panels (b)–(d), the relative fractions of energy loss $\Delta L_{A}^{\text{wav}}/L_{A,pho}^{\text{wav}}$ are plotted for each case (panel (b): $P_{\text{rot}} = 24$ day, panel (c): $P_{\text{rot}} = 12$ day, and panel (d): $P_{\text{rot}} = 6$ day), where total loss $\Delta L_{A,\text{tot}}^{\text{wav}}/L_{A,pho}^{\text{wav}}$ turbulence loss
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Alfvén-wave energy loss in the chromosphere. Panel (a): Alfvén-wave luminosities vs. height for different rotation rates: $P_{\text{rot}} = 24$ day (red solid line), $P_{\text{rot}} = 12$ day (green dashed–dotted line), and $P_{\text{rot}} = 6$ day (blue dashed line). Panel (b): normalized energy-loss fraction $\Delta L_{A,pho}^{\text{wav}}/L_{A,pho}^{\text{wav}}$ (black solid line) for $P_{\text{rot}} = 24$ day. Energy loss through turbulence and mode conversion are also shown with a red dashed line and a blue dashed–dotted line, respectively. Panels (c) and (d): same as panel (b) now for $P_{\text{rot}} = 12$ day and $P_{\text{rot}} = 6$ day, respectively.</figcaption></figure>

$\Delta L_{A,turb}/L_{A,pho}^{\text{wav}}$, and mode-conversion loss $\Delta L_{A,r \leftrightarrow \theta}/L_{A,pho}^{\text{wav}}$ are shown by the black solid, red dashed, and blue dashed–dotted lines, respectively. The black lines in panels (b)–(d) show that more than 90% of the Alfvén-wave energy flux is lost in the chromosphere. Moreover, the energy-loss fraction is larger for the faster-rotating case. Figure 13 shows the trend of the wave energy-loss fraction measured in the corona ($r - R_{\odot} = 10$ Mm), in which a larger loss fraction is clearly seen for the faster rotator. The aforementioned saturation of coronal Alfvén-wave energy flux $L_{A,cor}^{\text{wav}}$ is caused by this enhanced dissipation in fast rotators.

Comparing the red lines in panels (b)–(d), the increased energy loss is attributed to the increased turbulent loss. Given that the magnetic filling factor increases with rotation rate, the increased turbulent loss is a natural consequence for the following two reasons:

1. When the filling factor is large, magnetic-flux expansion is suppressed because one flux tube merges with the adjacent one after small expansion. The vortex size of turbulence (correlation length) is expected to expand with the flux tube, and thus remains small for fast rotators that are expected to have large open-flux filling factors. Therefore, Alfvén waves dissipate on smaller timescales, or equivalently, dissipate more quickly.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Fraction of Alfvén-wave energy loss in the chromosphere ($\Delta L_{A,pho}^{\text{wav}}/L_{A,pho}^{\text{wav}}$ measured at $r - R_{\odot} = 10$ Mm) as a function of rotation rate.</figcaption></figure>
2. The coronal magnetic field is stronger for cases with faster rotation. As a result, Alfvén-wave reflection at the transition region is enhanced. Therefore, Alfvén-wave turbulence (triggered by the Alfvén-wave reflection), should be promoted.

Although turbulent dissipation plays an important role in our model, we expect that our conclusion is not strongly affected by the amount of turbulent dissipation. Let us consider an extreme case with no turbulent dissipation. Alfvén waves propagate through the chromosphere with less dissipation and are strongly reflected at the transition region (Cranmer & van Ballegooijen 2005; Verdini & Velli 2007; Réville et al. 2018). Reflected Alfvén waves propagate backward without turbulent dissipation and reach the stellar surface. Since the net upward Alfvén-wave energy flux is reduced if more downward Alfvén waves are present, reduced turbulent dissipation leads to reduced $L_{A,pho}^{\text{wav}}$, which can also cause the saturation of the mass-loss rate. A detailed parameter survey on the turbulent correlation length should be done to test the above hypothesis.

We note that, because radiation dominates the energetics of the chromosphere, any dissipated energy is quickly radiated away. In this sense, our model is consistent with that of Suzuki et al. (2013), with the saturation of the mass-loss rate being compensated by an enhanced radiative loss.

To summarize our analysis, we have shown that a larger fraction of Alfvén-wave energy is dissipated in the chromosphere for faster-rotating cases. As a result, even though the photospheric Alfvén-wave luminosity follows a power-law relation in $\Omega_*$, the coronal-base Alfvén-wave luminosity saturates with increasing rotation, thus producing a saturation in the mass-loss rate.

# 5. Discussion

## 5.1. Comparison with Cranmer & Saar (2011)

A standard theoretical model of the stellar-wind mass-loss rate, for low-mass stars, is given by Cranmer & Saar (2011). Although both our model and that of Cranmer & Saar (2011) are based on Alfvén-wave heating and are calibrated by solar wind observations, the rotation dependence of the mass-loss rate is different. For example, when $P_{\text{rot}} = 2$ day, the Cranmer & Saar (2011) model yields a mass-loss rate that is 100 times larger than that predicted by our model. There are three factors that explain this difference.

1. Cranmer & Saar (2011) assumed a steeper dependence of $f_*^{\text{open}}$ on Ro with the exponent ranging between -2.5 and -3.4 in the unsaturated regime, while our model assumes much weaker dependence: $f_*^{\text{open}} \propto \text{Ro}^{-1.2}$. For example, when $\text{Ro}/\text{Ro}_\odot = 0.1$, $f_*^{\text{open}} = 0.36$ in Cranmer & Saar (2011), while $f_*^{\text{open}} = 0.016$ in this work. Given that the mass-loss rate approximately scales as $\propto (f_*^{\text{open}})^{5/7}$, this discrepancy yields a factor of 9.3 difference between Cranmer & Saar (2011) and our model. Since many observational aspects of stellar magnetism/winds are unresolved, theoretical mass-loss rates remain uncertain by around a factor of 10.

2. Cranmer & Saar (2011) employed a simplified model of Alfvén-wave propagation. Although their model also considers the turbulent dissipation of Alfvén waves, the difference between their wave equations and those used

in this work may lead to discrepancies in the resulting mass-loss rates. This hypothesis should be tested in the future by directly comparing the coronal wave energy between our model and that of Cranmer & Saar (2011).

3. Cranmer & Saar (2011) assumed that the wind speed is constant regardless of the open-flux filling factor. In reality, even without rotational acceleration, the wind speed tends to be higher for larger open-flux filling factors, which is explained as follows. Faster rotators exhibit larger coronal Alfvén speed that allows for more heat to be deposited beyond the sonic point. Given that the kinetic energy flux of the wind $(\propto \rho v_r^3)$ is constant, the mass-loss rate $(\propto \rho v_r)$ becomes smaller as the wind velocity $v_r$ increases in response to enhanced heating in the supersonic region.

## 5.2. Magnetic Transient Events

Our model assumes that the global magnetic structure is invariant on timescales of the stellar-wind acceleration. However, actual stellar magnetic fields can evolve in comparable or even shorter timescales than the wind acceleration. Specifically, the large-scale shuffling of magnetic field lines by super-granular motions is observed to cause magnetic reconnection and open closed magnetic features (Fisk et al. 1999; Antiochos et al. 2011; Moore et al. 2011; Higginson et al. 2017). Although the reconnection/loop-opening process is unlikely to be able to drive the majority of the solar wind (see, e.g., Cranmer & van Ballegooijen 2010; Lionello et al. 2016), it may play an important role for more active stars. Additionally, the reconnection/loop-opening process can work indirectly. For example, if the open magnetic field (carrying a quasi-steady wind) is rapidly connecting to closed loops with a high temperature, then the wind properties are determined by the closed loop temperature (O’Fionnagáin & Vidotto 2018). Such reconnection can also feed magnetohydrodynamic waves in addition to the surface granular motion (Cranmer 2018).

Eruptive processes such as CMEs can also be important in active stars (Aarnio et al. 2012; Drake et al. 2013). According to Cranmer (2017), CMEs could be a dominant source of mass loss for moderately faster rotators than the Sun. For much younger, much faster rotators, the centrifugally supported “slingshot prominences” are also expected to be present (Collier Cameron & Robinson 1989a, 1989b), and are likely to play a significant role in mass loss and magnetic braking. Recently Jardine & Collier Cameron (2019) have extended the $F_{\text{X}}-\dot{M}_{\text{w}}$ relation of Wood et al. (2014) to more active stars, based on mass-loss rates estimated from slingshot prominences. They show a significant mass loss through prominence ejection for such rapid rotators. The role of these eruptive processes should be taken into account in future works.

## 5.3. Comparison with Observations of Spin Evolution

In spite of successfully reproducing the stellar spin-down $\Omega_* \propto t^{-0.55}$, there exist several discrepancies between our model and stellar observations of spin evolution. As already mentioned in Section 3.2, the torque is smaller than the empirical value from stellar observation (Matt et al. 2015) but is consistent with solar wind observations (Finley et al. 2018). Since our model is calibrated by the solar wind, the deviation of our model from that of Matt et al. (2015) might be a result of
the solar magnetic field having an unusual character. However, reconstructions of the solar open magnetic flux from the last 9000 yr also recover the same solar wind torque as our model (see Finley et al. 2019).

Indeed, recent asteroseismic observations indicate that the solar dynamo could be in transition (Metcalfe et al. 2016). The deviation between our model and that of Matt et al. (2015) would be explained if the dynamo transition works to reduce the amount of open magnetic flux, by a factor of 2.4 from the canonical value. This hypothesis is consistent with the observed breakdown of gyrochronology (van Saders et al. 2016), and some spin-down models already take this effect into account (e.g., Garraffo et al. 2018). However, we must note that, from the perspective of dynamo simulations, large-scale field diminishing at $Ro > 1$ is not supported (Strugarek et al. 2018; Warnecke 2018; Guerrero et al. 2019). Thus, the small solar torque could be attributed to another mechanism.

### 5.4. Implications for Stars in the Saturated Regime

We have assumed that the filling factor of the open-flux regions monotonically increases with rotation rate, Equation (22), which is derived from stars in the slow-rotator regime of See et al. (2019). However, this power-law relation may be modified for fast rotators; the filling factor of the open regions may saturate at rapid rotation rates because a large fraction of the surface is expected to be covered by closed loops (these closed loops are thought to provide the observed coronal X-ray flux). This modification would affect the magnetic field strength in the chromosphere and the low corona, and changes the vertical profile of the Alfvén velocity there. A different profile of $v_A$ may enhance the transmitted fraction of Alfvénic waves through the transition region (Suzuki et al. 2013), which could increase the mass-loss rate in the fast-rotator regime.

In this work, we focused on the unsaturated regime of magnetic activity. However, many young low-mass stars (especially M dwarfs) lie in the saturated activity regime (e.g., Wright & Drake 2016; See et al. 2019). Given that the total open flux might be constant for stars in the saturated regime, our simulation results yield several implications for the winds of these stars. We expect that the mass-loss rate should be constant in the saturated regime. As shown in Figure 10, as long as the total open flux is fixed, increasing the rotation rate does not yield larger mass-loss rates. Instead, the increased rotational energy is used purely for wind acceleration. When the wind velocity is larger, $v_{r,A}$ (wind velocity at Alfvén point) should also be larger. According to the analytical expression of the torque in Equation (45), a larger $v_{r,A}$ yields a smaller $\tau_w/\Omega_*$, Thus, in the saturated regime, the torque could have a weaker-than-linear dependence on $\Omega_*$, In the future, this prediction could be directly tested by numerical simulations of stars in the saturated regime.

Numerical computations were carried out on Cray XC50 and PC cluster at Center for Computational Astrophysics, National Astronomical Observatory of Japan. M.S. is supported by Grant-in-Aid for Japan Society for the Promotion of Science (JSPS) Fellows and by the NINS program for cross-disciplinary study (grant Nos. 01321802 and 01311904) on Turbulence, Transport, and Heating Dynamics in Laboratory and Solar/ Astrophysical Plasmas: “SoLaBo-X.” T.K.S. is supported in part by Grants-in-Aid for Scientific Research from the MEXT of Japan, 17H01105. S.P.M., V.S., and A.J.F. are supported by the European Research Council, under the European Union’s

Horizon 2020 research and innovation program (agreement No. 682393, AWESoMeStars). A.A.V. is supported by the European Research Council, under the European Union’s Horizon 2020 research and innovation program (agreement No. 817540, ASTROFLOW). A.S. and A.S.B. acknowledge funding by ERC WHOLESUN 810218 grant, INSU/PNST, CNES-PLATO and CNES Solar Orbiter. A.S. acknowledges funding from the Programme National de Planétologie (PNP). V.R. is funded by the ERC SLOW_SOURCE project (SLOW_SOURCE—DLV-819189). This work benefited from discussions within the international team “The Solar and Stellar Wind Connection: Heating processes and angular-momentum loss,” supported by the International Space Science Institute (ISSI).

## Appendix A
### Derivation of Basic Equations

Except the nonideal terms (gravity, radiative loss, thermal conduction, and turbulent dissipation), our basic equations Equations (3)–(10) are derived from the typical ideal MHD equations as follows. Given the metrics of

$$h_r = 1, \quad h_\theta = rf^{\text{open}}, \quad h_\phi = r, \qquad (\text{A1})$$

and considering a one-dimensional system ($\partial/\partial\theta = \partial/\partial\phi = 0$), the nabla operators are expressed as follows:

$$\nabla \psi = \frac{\partial \psi}{\partial r} e_r, \qquad (\text{A2a})$$

$$\nabla \cdot \mathbf{A} = \frac{1}{r^2 f} \frac{\partial}{\partial r} (r^2 f A_r), \qquad (\text{A2b})$$

$$\nabla \times \mathbf{A} = -\mathbf{e}_{\theta} \frac{1}{r} \frac{\partial}{\partial r} (rA_{\phi}) + \mathbf{e}_{\phi} \frac{1}{rf} \frac{\partial}{\partial r} (rfA_{\theta}), \qquad (\text{A2c})$$

where, for simplicity, we denote $f^{\text{open}}$ as $f$. Using these expressions, each basic equation is derived in a straightforward manner. For example, the inertial and Lorentz forces in the equation of motion are written explicitly as

$$\begin{aligned} (\rho \mathbf{v} \cdot \nabla) \mathbf{v} &= \rho (\nabla \times \mathbf{v}) \times \mathbf{v} + \frac{1}{2} \rho \nabla (v^2) \\ &= e_r \left[ \rho v_r \frac{\partial}{\partial r} v_r - \rho v_\theta^2 \frac{d}{dr} \ln(rf) - \rho v_\phi^2 / r \right] \\ &\quad + \mathbf{e}_\theta \frac{\rho v_r}{rf} \frac{\partial}{\partial r} (rf v_\theta) + \mathbf{e}_\phi \frac{\rho v_r}{r} \frac{\partial}{\partial r} (r v_\phi), \end{aligned} \qquad (\text{A3})$$

$$\begin{aligned} \frac{1}{4\pi}(\nabla \times \mathbf{B}) \times \mathbf{B} = & e_r \left[ -\frac{B_\theta}{4\pi rf} \frac{\partial}{\partial r} (B_\theta rf) - \frac{B_\phi}{4\pi r} \frac{\partial}{\partial r} (rB_\phi) \right] \\ & + \mathbf{e}_\theta \frac{B_r}{4\pi rf} \frac{\partial}{\partial r} (rfB_\theta) + \mathbf{e}_\phi \frac{B_r}{4\pi r} \frac{\partial}{\partial r} (rB_\phi). \end{aligned} \qquad (\text{A4})$$

Similarly, the rotation of the electromotive force is

$$\begin{aligned} \nabla \times (\mathbf{v} \times \mathbf{B}) = & -\mathbf{e}_{\theta} \frac{1}{r} \frac{\partial}{\partial r} [r(v_r B_{\theta} - v_{\theta} B_r)] \\ & - \mathbf{e}_{\phi} \frac{1}{rf} \frac{\partial}{\partial r} [rf(v_r B_{\phi} - v_{\phi} B_r)]. \end{aligned} \qquad (\text{A5})$$

After these calculations, one can obtain the basic equations after rewriting in a conservation form.
## Appendix B

# An Analytical Formulation of Coronal Rotational Luminosity

The coronal rotational luminosity $L_{A,cor}^{rot}$ is obtained analytically based on Weber–Davis solution. We begin with the analytical expression of (time-averaged) $v_{\phi}$ and $B_{\phi}$:

$$v_{\phi} = r\Omega_{*}\frac{M_{A}^{2}L/(r^{2}\Omega_{*}) - 1}{M_{A}^{2} - 1}, \quad B_{\phi} = \frac{B_{r}}{v_{r}}(v_{\phi} - r\Omega_{*}). \quad (B1)$$

Near the coronal base where the wind velocity is negligibly small, we can approximate $v_{\phi}$ and $B_{\phi}$ to the first order of $M_{A}^{2}$ as

$$v_{\phi} \approx r\Omega_{*}\left(1 - \frac{M_{A}^{2}r_{A}^{2}}{r^{2}}\right), \quad B_{\phi} \approx -r\Omega_{*}\frac{B_{r}M_{A}^{2}r_{A}^{2}}{v_{r}r^{2}}. \quad (B2)$$

To the first order of $M_{A}^{2}$, the rotational luminosity at the coronal base is given as

$$L_{A,cor}^{rot} = -4\pi r^2 f^{open} \cdot \frac{B_r}{4\pi} v_{\phi} B_{\phi} \bigg|_{cor} \approx r_A^2 \Omega_* \dot{M}_w = \tau_w \Omega_* \quad (B3)$$

It is interesting to see that the rotational luminosity is approximated by $\tau_{w}\Omega_{*}$. Using Equation (44), the above formulation is further simplified as

$$L_{A,cor}^{rot} \approx \frac{\Phi_{open}^{2}}{16\pi^{2}v_{r,A}} \Omega_{*}^{2} \propto \Omega_{*}^{3.83}, \quad (B4)$$

where we have used $\Phi_{open} \propto \Omega_{*}^{1.2}$ and $v_{r,A} \propto \Omega_{*}^{0.57}$.

## ORCID iDs

Munehito Shoda https://orcid.org/0000-0002-7136-8190

Takeru K. Suzuki https://orcid.org/0000-0001-9734-9601

Sean P. Matt https://orcid.org/0000-0001-9590-2274

Steven R. Cranmer https://orcid.org/0000-0002-3699-3134

Aline A. Vidotto https://orcid.org/0000-0001-5371-2675

Antoine Strugarek https://orcid.org/0000-0002-9630-6463

Victor See https://orcid.org/0000-0001-5986-3423

Victor Réville https://orcid.org/0000-0002-2916-3837

Adam J. Finley https://orcid.org/0000-0002-3020-9409

Allan Sacha Brun https://orcid.org/0000-0002-1729-8267

## References

* Aarnio, A. N., Matt, S. P., & Stassun, K. G. 2012, ApJ, 760, 9

* Ahuir, J., Brun, A. S., & Strugarek, A. 2020, A&A, 635, A170

* Airapetian, V. S., Barnes, R., Cohen, O., et al. 2020, IJASB, 19, 136

* Alfvén, H. 1947, MNRAS, 107, 211

* Allan, A., & Vidotto, A. A. 2019, MNRAS, 490, 3760

* Angus, R., Aigrain, S., Foreman-Mackey, D., & McQuillan, A. 2015, MNRAS, 450, 1787

* Antiochos, S. K., Mikić, Z., Titov, V. S., Lionello, R., & Linker, J. A. 2011, ApJ, 731, 112

* Antolin, P., Shibata, K., Kudoh, T., Shiota, D., & Brooks, D. 2008, ApJ, 688, 669

* Arge, C. N., & Pizzo, V. J. 2000, JGR, 105, 10465

* Argiroffi, C., Reale, F., Drake, J. J., et al. 2019, NatAs, 3, 742

* Bale, S. D., Pulupa, M., Salem, C., Chen, C. H. K., & Quataert, E. 2013, ApJL, 769, L22

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

* Barnes, S. A. 2007, ApJ, 669, 1167

* Barnes, S. A. 2010, ApJ, 722, 222

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

* Belcher, J. W., & MacGregor, K. B. 1976, ApJ, 210, 498

* Bogdan, T. J., Carlsson, M., Hansteen, V. H., et al. 2003, ApJ, 599, 626

* Brun, A. S., & Browning, M. K. 2017, LRSP, 14, 4

* Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073

* Carbone, V., Marino, R., Sorriso-Valvo, L., Noullez, A., & Bruno, R. 2009, PhRvL, 103, 061102

* Chandran, B. D. G., & Perez, J. C. 2019, JPIPh, 85, 905850409

* Chen, C. H. K., Bale, S. D., Bonnell, J. W., et al. 2020, ApJS, 246, 53

* Chitta, L. P., van Ballegooijen, A. A., Rouppe van der Voort, L., DeLuca, E. E., & Kariyappa, R. 2012, ApJ, 752, 48

* Coleman, P. J., Jr. 1968, ApJ, 153, 371

* Collier Cameron, A., & Robinson, R. D. 1989a, MNRAS, 236, 57

* Collier Cameron, A., & Robinson, R. D. 1989b, MNRAS, 238, 657

* Cranmer, S. R. 2009, LRSP, 6, 3

* Cranmer, S. R. 2017, ApJ, 840, 114

* Cranmer, S. R. 2018, ApJ, 862, 6

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

* Cranmer, S. R., & van Ballegooijen, A. A. 2005, ApJS, 156, 265

* Cranmer, S. R., & van Ballegooijen, A. A. 2010, ApJ, 720, 824

* Cranmer, S. R., van Ballegooijen, A. A., & Edgar, R. J. 2007, ApJS, 171, 520

* Cranmer, S. R., & Winebarger, A. R. 2019, ARA&A, 57, 157

* Davies, G. R., Chaplin, W. J., Farr, W. M., et al. 2015, MNRAS, 446, 2959

* de Wijn, A. G., Lites, B. W., Berger, T. E., et al. 2008, ApJ, 684, 1469

* Del Zanna, L., Matteini, L., Landi, S., Verdini, A., & Velli, M. 2015, JPIPh, 81, 325810102

* Delorme, P., Collier Cameron, A., Hebb, L., et al. 2011, MNRAS, 413, 2218

* Dmitruk, P., Matthaeus, W. H., Milano, L. J., et al. 2002, ApJ, 575, 571

* Dobrowolny, M., Mangeney, A., & Veltri, P. 1980, PhRvL, 45, 144

* Drake, J. J., Cohen, O., Yashiro, S., & Gopalswamy, N. 2013, ApJ, 764, 170

* Elsässer, W. M. 1950, PhRv, 79, 183

* Finley, A. J., Deshmukh, S., Matt, S. P., Owens, M., & Wu, C.-J. 2019, ApJ, 883, 67

* Finley, A. J., & Matt, S. P. 2017, ApJ, 845, 46

* Finley, A. J., & Matt, S. P. 2018, ApJ, 854, 78

* Finley, A. J., Matt, S. P., & See, V. 2018, ApJ, 864, 125

* Fisk, L. A., Schwadron, N. A., & Zurbuchen, T. H. 1999, JGR, 104, 19765

* Folsom, C. P., Fionnagáin, D. Ó, Fossati, L., et al. 2020, A&A, 633, A48

* Fujiki, K., Tokumaru, M., Iju, T., Hakamada, K., & Kojima, M. 2015, SoPh, 290, 2491

* 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. 2016, ApJL, 833, L4

* Garraffo, C., Drake, J. J., Dotter, A., et al. 2018, ApJ, 862, 90

* Goodman, M. L., & Judge, P. G. 2012, ApJ, 751, 75

* Güdel, M. 2007, LRSP, 4, 3

* Güdel, M., Guinan, E. F., & Skinner, S. L. 1997, ApJ, 483, 947

* Gudiksen, B. V., & Nordlund, Å. 2005, ApJ, 618, 1020

* Guerrero, G., Zaire, B., Smolarkiewicz, P. K., et al. 2019, ApJ, 880, 6

* Hammer, R. 1982, ApJ, 259, 779

* Hansteen, V. H., & Leer, E. 1995, JGR, 100, 21577

* Hansteen, V. H., & Velli, M. 2012, SSRv, 172, 89

* Heyvaerts, J., & Priest, E. R. 1983, A&A, 117, 220

* Higginson, A. K., Antiochos, S. K., DeVore, C. R., Wyper, P. F., & Zurbuchen, T. H. 2017, ApJ, 837, 113

* Hirzberger, J., Bonet, J. A., Vázquez, M., & Hanslmeier, A. 1999, ApJ, 515, 441

* Hollweg, J. V., Jackson, S., & Galloway, D. 1982, SoPh, 75, 35

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

* Hossain, M., Gray, P. C., Pontius, D. H., Jr., Matthaeus, W. H., & Oughton, S. 1995, PhFl, 7, 2886

* Hotta, H., Rempel, M., & Yokoyama, T. 2016, Sci, 351, 1427

* Howes, G. G., & Nielson, K. D. 2013, PhPl, 20, 072302

* Iijima, H. 2016, PhD thesis, The Univ. Tokyo

* Irwin, J., & Bouvier, J. 2009, in IAU Symp. 258, The Ages of Stars, ed. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse (Cambridge: Cambridge Univ. Press), 363

* Ishikawa, R. T., Katsukawa, Y., Oba, T., et al. 2020, ApJ, 890, 138

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

* Jardine, M., & Collier Cameron, A. 2019, MNRAS, 482, 2853

* Johnstone, C. P. 2017, A&A, 598, A24

* Johnstone, C. P., Güdel, M., Brott, I., & Lüftinger, T. 2015a, A&A, 577, A28

* Johnstone, C. P., Güdel, M., Stökl, A., et al. 2015b, ApJL, 815, L12

* Johnstone, C. P., Khodachenko, M. L., Lüftinger, T., et al. 2019, A&A, 624, L10

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

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

* Kraft, R. P. 1967, ApJ, 150, 551
Kraichnan, R. H. 1965, PhFl, 8, 1385

Kudoh, T., & Shibata, K. 1999, ApJ, 514, 493

Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge: Cambridge Univ. Press), 452

Lammer, H., Selsis, F., Chassefieèle, E., et al. 2010, AsBio, 10, 45

Leer, E., Holzer, T. E., & Fla, T. 1982, SSRv, 33, 161

Leighton, R. B. 1969, ApJ, 156, 1

Linker, J. A., Caplan, R. M., Downs, C., et al. 2017, ApJ, 848, 70

Lionello, R., Török, T., Titov, V. S., et al. 2016, ApJL, 831, L2

Magaudda, E., Stelzer, B., Covey, K. R., et al. 2020, arXiv:2004.02904

Magyar, N., VanDoorsselaere, T., & Goossens, M. 2017, NatSR, 7, 14820

Matsumoto, T., & Shibata, K. 2010, ApJ, 710, 1857

Matt, S., & Pudritz, R. E. 2008, ApJ, 678, 1109

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

Matt, S. P., MacGregor, K. B., Pinsonneault, M. H., & Greene, T. P. 2012, ApJL, 754, L26

Matthaeus, W. H., Zank, G. P., Oughton, S., Mullan, D. J., & Dmitruk, P. 1999, ApJL, 523, L93

Meibom, S., Barnes, S. A., Platais, I., et al. 2015, Natur, 517, 589

Metcalfe, T. S., Egeland, R., & van Saders, J. 2016, ApJL, 826, L2

Moore, R. L., Sterling, A. C., Cirtain, J. W., & Falconer, D. A. 2011, ApJL, 731, L18

Moriyasu, S., Kudoh, T., Yokoyama, T., & Shibata, K. 2004, ApJL, 601, L107

Notsu, Y., Maehara, H., Honda, S., et al. 2019, ApJ, 876, 58

O’Fionnagáin, D., & Vidotto, A. A. 2018, MNRAS, 476, 2465

Oba, T., Iida, Y., & Shimizu, T. 2020, ApJ, 890, 141

Oba, T., Riethmüller, T. L., Solanki, S. K., et al. 2017, ApJ, 849, 7

Osterbrock, D. E. 1961, ApJ, 134, 347

Pantolmos, G., & Matt, S. P. 2017, ApJ, 849, 83

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

Parker, E. N. 1988, ApJ, 330, 474

Perez, J. C., & Chandran, B. D. G. 2013, ApJ, 776, 124

Pinto, R. F., & Rouillard, A. P. 2017, ApJ, 838, 89

Podesta, J. J., Roberts, D. A., & Goldstein, M. L. 2007, ApJ, 664, 543

Radick, R. R., Thompson, D. T., Lockwood, G. W., Duncan, D. K., & Baggett, W. E. 1987, ApJ, 321, 459

Rappazzo, A. F., Velli, M., Einaudi, G., & Dahlburg, R. B. 2008, ApJ, 677, 1348

Reiners, A., Basri, G., & Browning, M. 2009, ApJ, 692, 538

Rempel, M. 2017, ApJ, 834, 10

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

Réville, V., Folsom, C. P., Strugarek, A., & Brun, A. S. 2016, ApJ, 832, 145

Réville, V., Tenerani, A., & Velli, M. 2018, ApJ, 866, 38

Réville, V., Velli, M., Panasenco, O., et al. 2020, ApJS, 246, 24

Ribas, I., Guinan, E. F., Güdel, M., & Audard, M. 2005, ApJ, 622, 680

Rosenthal, C. S., Bogdan, T. J., Carlsson, M., et al. 2002, ApJ, 564, 508

Saar, S. H. 2001, in ASP Conf. Ser. 223, Recent Measurements of (and Inferences About) Magnetic Fields on K and M Stars (CD-ROM Directory: Contribs/ Saar1), ed. R. J. Garcia Lopez, R. Rebolo, & M. R. Zapaterio Osorio (San Francisco, CA: ASP), 292

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

Salem, C., Hubert, D., Lacombe, C., et al. 2003, ApJ, 585, 1147

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

Schatzman, E. 1962, AnAp, 25, 18

Schröder, K.-P., & Cuntz, M. 2005, ApJL, 630, L73

Schunker, H., & Cally, P. S. 2006, MNRAS, 372, 551

See, V., Jardine, M., Vidotto, A., et al. 2018, MNRAS, 474, 536

See, V., Lehmann, L., Matt, S. P., & Finley, A. J. 2020, ApJ, 894, 69

See, V., Matt, S. P., Finley, A. J., et al. 2019, ApJ, 886, 120

Shoda, M., Suzuki, T. K., Asgari-Targhi, M., & Yokoyama, T. 2019, ApJL, 880, L2

Shoda, M., & Yokoyama, T. 2018, ApJ, 854, 9

Shoda, M., Yokoyama, T., & Suzuki, T. K. 2018a, ApJ, 853, 190

Shoda, M., Yokoyama, T., & Suzuki, T. K. 2018b, ApJ, 860, 17

Skumanich, A. 1972, ApJ, 171, 565

Spitzer, L., & Härm, R. 1953, PhRv, 89, 977

Strugarek, A., Beaudoin, P., Charbonneau, P., & Brun, A. S. 2018, ApJ, 863, 35

Strugarek, A., Brun, A. S., Matt, S. P., & Réville, V. 2014, ApJ, 795, 86

Suzuki, T. K. 2007, ApJ, 659, 1592

Suzuki, T. K. 2018, PASJ, 70, 34

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

Suzuki, T. K., & Inutsuka, S.-i. 2005, ApJL, 632, L49

Suzuki, T. K., & Inutsuka, S.-I. 2006, JGRA, 111, 6101

Tenerani, A., & Velli, M. 2013, JGRA, 118, 7507

Toriumi, S., & Wang, H. 2019, LRSP, 16, 3

Tsuneta, S., Ichimoto, K., Katsukawa, Y., et al. 2008, ApJ, 688, 1374

van Ballegooijen, A. A., & Asgari-Targhi, M. 2016, ApJ, 821, 106

van Ballegooijen, A. A., & Asgari-Targhi, M. 2017, ApJ, 835, 10

van Ballegooijen, A. A., Asgari-Targhi, M., Cranmer, S. R., & DeLuca, E. E. 2011, ApJ, 736, 3

van Ballegooijen, A. A., Nisenson, P., Noyes, R. W., et al. 1998, ApJ, 509, 435

van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ, 782, 81

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

Velli, M. 1994, ApJL, 432, L55

Verdini, A., Grappin, R., & Montagud-Camps, V. 2019, SoPh, 294, 65

Verdini, A., Grappin, R., Pinto, R., & Velli, M. 2012, ApJL, 750, L33

Verdini, A., & Velli, M. 2007, ApJ, 662, 669

Vidotto, A. A., & Cleary, A. 2020, MNRAS, 494, 2417

Vidotto, A. A., Fares, R., Jardine, M., et al. 2012, MNRAS, 423, 3285

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

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

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

Warnecke, J. 2018, A&A, 616, A72

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

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

Withbroe, G. L. 1988, ApJ, 325, 442

Wood, B. E., Müller, H.-R., Redfield, S., & Edelman, E. 2014, ApJL, 781, L33

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. 2016, Natur, 535, 526

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