# Effect of Differential Rotation on the Magnetic Braking of Low-mass and Solar-like Stars: A Proof-of-concept Study

Lewis G. Ireland, Sean P. Matt, Charlie R. Davey, Owain L. Harris, Tobias

W. Slade-Harajda, Adam J. Finley, Claudio Zanni

► To cite this version:

Lewis G. Ireland, Sean P. Matt, Charlie R. Davey, Owain L. Harris, Tobias W. Slade-Harajda, et al.  
Effect of Differential Rotation on the Magnetic Braking of Low-mass and Solar-like Stars: A Proof-of-  
concept Study. The Astrophysical Journal, 2022, 925, 10.3847/1538-4357/ac3a71.  
insu-03745339

HAL Id: insu-03745339

https://insu.hal.science/insu-03745339

Submitted on 4 Aug 2022

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

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

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

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

# Effect of Differential Rotation on the Magnetic Braking of Low-mass and Solar-like Stars: A Proof-of-concept Study

Lewis G. Ireland<fnref n="1" />, Sean P. Matt<fnref n="1" />, Charlie R. Davey<fnref n="1" />, Owain L. Harris<fnref n="1" />, Tobias W. Slade-Harajda<fnref n="1" />,
Adam J. Finley<fnref n="2" />, and Claudio Zanni<fnref n="3" />

Received 2021 August 27; revised 2021 November 15; accepted 2021 November 15; published 2022 January 27

## Abstract

On the main sequence, low-mass and solar-like stars are observed to spin down over time, and magnetized stellar winds are thought to be predominantly responsible for this significant angular momentum loss. Previous studies have demonstrated that the wind torque can be predicted via formulations dependent on stellar properties, such as magnetic field strength and geometry, stellar radius and mass, wind mass-loss rate, and stellar rotation rate. Although these stars are observed to experience surface differential rotation, torque formulations so far have assumed solid-body rotation. Surface differential rotation is expected to affect the rotation of the wind and thus the angular momentum loss. To investigate how differential rotation affects the torque, we use the PLUTO code to perform 2.5D magnetohydrodynamic, axisymmetric simulations of stellar winds, using a colatitude-dependent surface differential rotation profile that is solar-like (i.e., rotation is slower at the poles than the equator). We demonstrate that the torque is determined by the average rotation rate in the wind so that the net torque is less than that predicted by assuming solid-body rotation at the equatorial rate. The magnitude of the effect is essentially proportional to the magnitude of the surface differential rotation, for example, resulting in a torque for the Sun that is ~20% smaller than predicted by the solid-body assumption. We derive and fit a semianalytic formulation that predicts the torque as a function of the equatorial spin rate, magnitude of differential rotation, and wind magnetization (depending on the dipolar magnetic field strength and mass-loss rate, combined).

*Unified Astronomy Thesaurus concepts: Stellar winds (1636); Stellar evolution (1599); Stellar magnetic fields (1610); Stellar rotation (1629); Magnetohydrodynamical simulations (1666); Main sequence stars (1000); Low mass stars (2050); Solar mass (1506)*

## 1. Introduction

The evolution of rotation rates of low-mass and solar-like stars on the main sequence (MS) is a complex function of their mass and age, as illustrated, for example, by rotation-period-mass diagrams from observations of large stellar clusters (see, e.g., Barnes 2003; Irwin & Bouvier 2009; Meibom et al. 2011; Bouvier et al. 2014, pp. 433–450; Stauffer et al. 2016; Davenport 2017). Generally, these stars are observed to spin down with time (Skumanich 1972), due to the extraction of angular momentum by magnetized stellar winds (Parker 1958; Schatzman 1962; Weber & Davis 1967; Mestel 1968). Models for this rotational evolution (e.g., Kawaler 1988; MacGregor & Brenner 1991; Bouvier et al. 1997; Bouvier 2008; Denissenkov et al. 2010; Matt et al. 2012; Reiners & Mohanty 2012; Gallet & Bouvier 2015; van Saders & Pinsonneault 2013; Gallet & Bouvier 2013; Lanzafame & Spada 2015; Matt et al. 2015; Blackman & Owen 2016; Gondoin 2017; Garraffo et al. 2018; Amard et al. 2019; Breimann et al. 2021) rely on theoretical formulations for stellar-wind torques, which depend on a variety of stellar properties, such as masses and radii, magnetic field strengths and topology, coronal temperatures, and mass-loss rates. Due to the importance of magnetized stellar winds for understanding MS rotational evolution, there is a need for the

development of formulations that predict the torque and include all of the important physical processes in stellar winds.

Early and analytic models of stellar winds and angular momentum loss (e.g., Parker 1958; Schatzman 1962; Parker 1965; Weber & Davis 1967; Mestel 1968; Kawaler 1988; Sauty et al. 2011) studied thermal-pressure-driven winds in the presence of magnetic fields, demonstrating that initially subsonic stellar winds accelerate beyond the sound speed and beyond the magnetic Alfvén wave speed (at the Alfvén radius). The angular momentum loss is enhanced by the presence of a magnetic field, and magnetocentrifugal processes in rapid rotators can even enhance the overall wind speed. Weber & Davis (1967) also demonstrated that the stellar-wind torque depends on the square of the Alfvén radius, which by mechanical analogy is often referred to as the “magnetic lever arm” length. More recently, stellar-wind dynamics have been studied numerically via magnetohydrodynamic (MHD) simulations, allowing for the comparison between self-consistently computed values of stellar-wind torque and those predicted via parameterizations of the Alfvén radius. Matt & Pudritz (2008) investigated how the torque depends on magnetic dipole field strength and mass-loss rate around slow rotators (see also Washimi & Shibata 1993; Cohen & Drake 2014). Subsequent works explored how the torques depend on magnetocentrifugal effects around fast rotators (Matt et al. 2012), magnetic fields with more complex magnetic geometries (Matt & Pudritz 2008; Réville et al. 2015; Garraffo et al. 2016; Finley & Matt 2017, 2018), coronal temperature (Pantolmos & Matt 2017), and magnetic cycle variations (Finley et al. 2018; Perri et al. 2018; Finley et al. 2019). So far, studies that have been

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK; L.G.Ireland@exeter.ac.uk</footnote>

<footnote n="2">²Department of Astrophysics-AIM, University of Paris-Saclay and University of Paris, CEA, CNRS, Gif-sur-Yvette Cedex F-91191, France</footnote>

<footnote n="3">³INAF—Osservatorio Astrofisico di Torino, Strada Osservatorio 20, I-10025 Pino Torinese, Italy</footnote>
used to derive formulations for stellar-wind torques have assumed that the stellar surface (or more specifically, the base of the wind-emitting region) rotates as a solid body.

Stars behave as rotating fluids, and so the surface rotation rate of a star is not always uniform, but instead can vary as a function of colatitude (see, e.g., Kitchatinov & Rüdiger 1999; Küker et al. 2011). Common trackers of surface differential rotation for the Sun include sunspots and Doppler measurements (see, e.g., Beck 2000 for further discussion and comparison between measurements). Typical measured values show the poles to be spinning ~30% slower than the equator (Snodgrass & Ulrich 1990). For other F, G, and K MS stars, the surface differential rotation has been inferred from, for example, Doppler imaging (see, e.g., Collier Cameron et al. 2002), spectral line Fourier transforms (see, e.g., Reiners & Schmitt 2002a), and photometric measurements of stellar periods (see, e.g., Reinhold et al. 2013). A commonly used prescription to describe the surface rotation rate as a function of colatitude, $\theta$, is $\Omega(\theta) = \Omega_{*,eq}(1 - \alpha \cos^2 \theta)$, where $\Omega_{*,eq}$ is the angular rotation rate at the equator, and the relative differential rotation rate is defined as $\alpha = \Delta\Omega/\Omega_{*,eq} \equiv (\Omega_{*,eq} - \Omega_{*,pole})/\Omega_{*,eq}$, where $\Omega_{*,pole}$ is the angular rotation rate at the poles. Thus, for the Sun, $\alpha \approx 0.3$. This relative differential rotation rate is observed to be approximately correlated with rotation period, with the vast majority of F, G, and K MS stars having $\alpha < 0.3$, and most observed samples of rapidly rotating stars (that rotate faster, on average, than the Sun) showing $\alpha \to 0$ (see, e.g., Collier Cameron 2002; Reiners & Schmitt 2002b; Barnes et al. 2005; Reiners 2006; Collier Cameron 2007; Balona & Abedigamba 2016). Most stars appear to be undergoing solar-like rotation, where the equator spins faster than the poles ($\alpha > 0$). However, it is possible that some stars are undergoing antisolar rotation, where the poles are faster than the equator ($\alpha < 0$), or more complex profiles such as cylindrical-banded rotation, which consists of alternating zonal jets (see, e.g., Brun et al. 2017 and references therein for a discussion of solar-like and antisolar differential rotation).

We expect a priori that differential rotation will influence the rotation of the wind, hence the stellar torque. In order to quantify and be able to predict the effects of differential rotation on the torque, we perform 74 2.5D axisymmetric stellar-wind simulations, changing $\Omega_{*,eq}$, $\alpha$ (restricted to solar-like cases), and the dipolar magnetic field strength. We show that the stellar-wind torque is simply predicted in terms of the average ("effective") rotation rate of the wind material. Because the wind emanates from a region surrounding the poles, the "effective" rotation rate is slower than that of the stellar equator, and so the net torque is smaller than would be predicted assuming a solid-body rotation at the rate of the equator. We find a formulation for predicting this "effective" rotation rate (and thus the torque) as a function of the differential rotation rate, and other stellar and wind parameters. In this study, we purely investigate how differential rotation can modify the rotation rate at the base of the wind and how this affects global angular momentum loss. We are concerned only about the rotation rate of open magnetic field lines, independent of the effects on the coronal dynamics and small-scale magnetic field structure.

In Section 2, we describe the numerical setup for our simulations, the initial conditions and final solutions, boundary conditions, normalizations, and parameter space explored in this study. In Section 3, we list our simulations, present their qualitative behavior, and introduce the global quantities we use,

such as the mass flow rate, the torque, and the unsigned magnetic flux. In Section 4, we introduce a torque formulation for the stellar wind that accounts for the differential rotation profile investigated in this parameter study. In Section 5, we discuss and conclude our findings.

## 2. Numerical Method

### 2.1. Numerical Setup

The simulations in this study are solved numerically governed by the following MHD equations:

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) &= 0, \\
\frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot \left[ \rho \mathbf{v} \mathbf{v} + \left( p + \frac{\mathbf{B} \cdot \mathbf{B}}{8\pi} \right) \mathbf{I} - \frac{\mathbf{B}\mathbf{B}}{4\pi} \right] &= \rho \mathbf{g}, \\
\frac{\partial E}{\partial t} + \nabla \cdot \left[ \left( E + p + \frac{\mathbf{B} \cdot \mathbf{B}}{8\pi} \right) \mathbf{v} - \frac{(\mathbf{v} \cdot \mathbf{B})\mathbf{B}}{4\pi} \right] &= \rho \mathbf{g} \cdot \mathbf{v}, \\
\frac{\partial \mathbf{B}}{\partial t} + \nabla \times (\mathbf{B} \times \mathbf{v}) &= 0,
\end{aligned}
\quad (1) $$

which represent the mass continuity, momentum, energy, and magnetic induction equation, respectively. Here, $\rho$ represents the mass density, $\mathbf{v}$ the velocity field, $p$ the thermal pressure, $\mathbf{B}$ the magnetic field, $\mathbf{I}$ the identity matrix, and $\mathbf{g} = -(GM_*/R^2)\hat{\mathbf{R}}$ the gravitational acceleration ($G$ is the gravitational constant, $M_*$ is the stellar mass, and $R$ is the spherical radius). The total energy density is written as

$$ E = \rho u + \rho \frac{\mathbf{v} \cdot \mathbf{v}}{2} + \frac{\mathbf{B} \cdot \mathbf{B}}{8\pi}, \quad (2) $$

where $u$ is specific internal energy (per unit mass). We use a polytropic wind, with an equation of state of the form $\rho u = p/(\gamma - 1)$, where $\gamma$ is the adiabatic index. We numerically solve Equations (1)–(2) using the PLUTO code (Mignone et al. 2007, 2012), using a finite-volume Godunov scheme and a linearized ROE Riemann solver (Roe 1981). We use the constrained transport method to control the divergence-free constraint $\nabla \cdot \mathbf{B} = 0$ (see Tóth et al. 2005).

All simulations are 2.5D (2D computational domain, with 3D vector components) and adopt the $(R, \theta)$ spherical coordinate system. Axisymmetry is assumed about the stellar rotation axis. We define the cylindrical radius using $r = R \sin \theta$, where $R$ is the spherical radius. We use a logarithmic grid for the radial direction, consisting of 256 grid cells and covering [1, 60]$R_*$, where $R_*$ is the stellar radius. We use a uniform grid for the $\theta$ direction, consisting of 512 grid cells and covering [0, $\pi$].

### 2.2. Initial Conditions and Final Solutions

We initialize the computational domain with a stellar corona and a magnetic field configuration. The stellar corona is initialized with a 1D spherically symmetric, isotropic, polytropic Parker wind solution, which is defined by the ratio of the stellar surface sound speed $c_{s,*} = (\gamma p_*/\rho_*)^{1/2}$ to the stellar surface escape velocity $v_{esc} = (2GM_*/R_*)^{1/2}$, where $p_*$ and $\rho_*$ are pressure and density at the stellar surface, respectively. We set $\gamma = 1.05$ (near isothermal), so the wind is heated during expansion without requiring explicit heating terms in the energy equation from
Equation (1) (see the discussion in Finley & Matt 2018 and references therein). The magnetic field configuration is initially purely dipolar and aligned with the axis of rotation. The R and θ components are described by

$$B_R(R, \theta) = B_* \left(\frac{R_*}{R}\right)^3 \cos \theta \quad (3)$$

and

$$B_{\theta}(R, \theta) = \frac{1}{2} B_{*} \left(\frac{R_{*}}{R}\right)^{3} \sin \theta, \quad (4)$$

respectively, where $B_*$ represents the polar stellar surface magnetic field strength.

We run each simulation until it reaches a steady state or, in some cases, an oscillating state (as described in Pantolmos & Matt 2017), which we regard as the final wind solution. The winds are characterized by a continuous and ubiquitous flow from the lower radial boundary (the “star”), through the grid, and out of the outer radial boundary, at speeds that exceed all information-carrying-wave speeds. Thus, the final solutions are insensitive to the initial conditions but are instead essentially entirely determined by the lower radial boundary conditions (and somewhat influenced by the axial boundary, at $\theta = 0$ and $\pi$), described in the following section.

## 2.3. Boundary Conditions

The computational domain is enclosed by a boundary filled with ghost cells, allowing for the implementation of boundary conditions. At the inner boundary, the density and thermal-pressure profiles from the Parker wind solution are kept fixed. The R magnetic field component is kept fixed to conserve the total stellar flux, but the $\theta/\phi$ magnetic field components are free to change via a linear extrapolation. The poloidal velocity $v_p$ is imposed to be parallel to the poloidal magnetic field $B_p$, as well as the continuity of the mass flux per magnetic flux (an axisymmetric MHD invariant),

$$\kappa = \frac{\rho v_p \cdot B_p}{|B_p|^2}, \quad (5)$$

along magnetic field lines, which ensures smooth inflow of the stellar wind (see, e.g., Ustyugova et al. 1999; Keppens & Goedbloed 2000). The stellar rotation at the inner boundary is set via the toroidal velocity,

$$v_{\phi} = r \Omega_{*}(\theta) + v_{p} \frac{B_{\phi}}{B_{p}}, \quad (6)$$

where $\Omega_*(\theta)$ is the stellar rotation rate at a given colatitude $\theta$. By using Equation (6), we impose the differential rotation to the rate of rotation of the magnetic surfaces. The toroidal speed of the plasma (which determines, e.g., the Doppler shifts) slightly deviates from this and therefore actually contains some very small degree of differential rotation (even for a solid rotator). On the other hand, because $v_p$ (the injection speed of the wind) is smaller than $r\Omega_*$ (at least for faster rotators) and $B_\phi/B_p$ should be also reasonably small, the two speeds (plasma and field) almost coincide. We adopt the following simple differential rotation profile for the stellar surface rotation rate:

$$\Omega_{*}(\theta) = \Omega_{*,eq}(1 - \alpha \cos^{2}\theta), \quad (7)$$

where $\Omega_{*,eq}$ is the angular rotation rate at the stellar equator, and $\alpha = \Delta\Omega/\Omega_{*,eq} \equiv (\Omega_{*,eq} - \Omega_{*,pole})/\Omega_{*,eq}$, where $\Omega_{*,pole}$ is the angular rotation rate at the poles. This rotation profile is kept fixed throughout the computation. At the outer boundary, all quantities have vanishing derivatives, e.g., $dv_R/dR = 0$, allowing them to flow outward from the computational domain.

## 2.4. Simulation Parameters

We investigate a parameter space where we systematically vary the following quantities:

1. Surface polar magnetic field strength, $B_*$; in practice, this is controlled via the input parameter $v_A/v_{esc}$, i.e., the ratio of the surface polar Alfvén velocity $v_A = B_*/(4\pi\rho_*)^{1/2}$ and the stellar escape velocity $v_{esc} = (2GM_*/R_*)^{1/2}$.

2. Stellar equatorial rotation rate, expressed as a fraction of the breakup rate, $f_{eq} = \Omega_{*,eq}R_*/v_{K,*}$.

3. Relative differential rotation rate, $\alpha = (\Omega_{*,eq} - \Omega_{*,pole})/\Omega_{*,eq}$.

Where $\Omega_{*,eq}$ and $\Omega_{*,pole}$ are the equatorial and polar stellar rotation rates, respectively. We fix $\gamma = 1.05$ and the stellar-wind sound speed at the stellar surface to $c_{s,*} = 0.25v_{esc}$ for all simulations.

## 2.5. Units and Normalization

We perform simulations in dimensionless units; here, we list normalization factors required to convert quantities into physical units representative of different types of stars. We express length in units of $R_*$, density in units of its value at the base of the corona $\rho_*$, and velocities in units of the stellar surface Keplerian velocity, $v_{K,*} = (GM_*/R_*)^{1/2}$ (where $M_*$ is the stellar mass). The following normalizations can then be derived from the aforementioned base dimensionalizations: time in units of $t_0 = R_*/v_{K,*}$, magnetic field strength in units of $B_0 = (4\pi\rho_*(v_{K,*}^2)^{1/2})$, mass flow rate in units of $\dot{M}_0 = \rho_* R_*^2 v_{K,*}$, torque in units of $\dot{J}_0 = \rho_* R_*^3 v_{K,*}^2$, and magnetic flux in units of $\Phi_0 = (4\pi\rho_* R_*^4 v_{K,*}^2)^{1/2}$.

By adopting solar values for radius, $R_\odot = 6.96 \times 10^{10}$ cm, mass $M_\odot = 1.99 \times 10^{33}$ g, and density at the base of the corona $\rho_\odot = 2.46 \times 10^{-16}$ g cm$^{-3}$ (see section 3), we can make direct comparisons with main-sequence solar-like and low-mass stars using the following normalizations:

$$\begin{align*}
v_{K,*} &= 437 \left(\frac{M_*}{M_\odot}\right)^{1/2} \left(\frac{R_*}{R_\odot}\right)^{-1/2} \text{ km s}^{-1}, \\
B_0 &= 2.43 \left(\frac{\rho_*}{\rho_\odot}\right)^{1/2} \left(\frac{M_*}{M_\odot}\right)^{1/2} \left(\frac{R_*}{R_\odot}\right)^{-1/2} \text{ G}, \\
t_0 &= 0.0184 \left(\frac{M_*}{M_\odot}\right)^{-1/2} \left(\frac{R_*}{R_\odot}\right)^{3/2} \text{ days}, \\
\dot{M}_0 &= 8.26 \times 10^{-13} \left(\frac{\rho_*}{\rho_\odot}\right) \left(\frac{M_*}{M_\odot}\right)^{1/2} \left(\frac{R_*}{R_\odot}\right)^{3/2} \text{ M}_\odot \text{ yr}^{-1}, \\
\dot{J}_0 &= 1.58 \times 10^{32} \left(\frac{\rho_*}{\rho_\odot}\right) \left(\frac{M_*}{M_\odot}\right) \left(\frac{R_*}{R_\odot}\right)^2 \text{ erg}, \\
\Phi_0 &= 1.18 \times 10^{22} \left(\frac{\rho_*}{\rho_\odot}\right)^{1/2} \left(\frac{M_*}{M_\odot}\right)^{1/2} \left(\frac{R_*}{R_\odot}\right)^{3/2} \text{ Mx.} \tag{8}
\end{align*}$$
Furthermore, assuming the solar fraction of the equatorial breakup rate $f_{\odot,eq} = 4.73 \times 10^{-3}$ (using a typical measured value of the solar equatorial sidereal rotation period, 24.47 days; Snodgrass & Ulrich 1990), the equatorial stellar rotation period ($P_{*,eq} = 2\pi/\Omega_{*,eq} \equiv 2\pi R_*/(f_{eq}v_{K,\star})$) can be expressed as

$$P_{*,eq} = 24.47 \left( \frac{f_{eq}}{f_{\odot,eq}} \right)^{-1} \left( \frac{M_{\star}}{M_{\odot}} \right)^{-1/2} \left( \frac{R_{\star}}{R_{\odot}} \right)^{3/2} \text{ days.} \quad (9)$$

### 3. Simulations of Stellar Winds

This study consists of 74 simulations in total, exploring the parameter space outlined in Section 2.4. Each simulation is ran for a period roughly corresponding to 16 stellar rotation periods for $f_{eq} = 0.05$ ($t_0 = 2000$). Typically, this corresponds to weeks or months in simulated time, where changes in the stellar rotation rate are negligible, justifying our use of fixed $f_{eq}$. All simulations evolve to a quasi-steady state in this simulation time period. Simulation input and output parameters can be found in Table 1.

In Table 2, we make direct comparisons between our simulations and examples of solar-like and low-mass stars by dimensionalizing our numerical simulation units using normalizations shown in Section 2.5. For $M_*/M_\odot = (0.10, 0.25, 0.50, 1.00)$, we list the range of surface magnetic field strengths and stellar rotation periods that are covered by our parameter space. For each stellar mass, we also list the ranges of mass-loss rate, torque, and open magnetic flux that result from these simulations. The stellar radius for each mass is taken from the grid of low-mass stellar models produced by Amard et al. (2019) (at 5 Gyr). For the $M_*=M_\odot$ case, we pick a value of $\rho_*$ so that the mass-loss rate of model 1 (the closest to a representative solar case) is $\dot{M} = \dot{M}_\odot \approx 1.78 \times 10^{-14} M_\odot \text{yr}^{-1}$, which is a reasonable value for the Sun.<fnref n="4" /> For $M_*=0.1 M_\odot$, we pick a value of $\rho_*$ so that the range of $B_*$ in our simulations ranges from 100–800 G, which is within the range bracketed by Zeeman-Doppler magnetogram observations (see, e.g., Morin et al. 2010; See et al. 2019a, 2019b).<fnref n="5" /> For the intermediate masses, $M_* = 0.25–0.50 M_\odot$, we choose interpolated values of $\rho_*$ assuming $\rho_* \propto M_*^{-3}$ (given the three orders of magnitude decrease in $\rho_*$ for one order of magnitude increase in stellar mass) for simplicity.

There are relatively few observational constraints on the mass-loss rate of low-mass and solar-like stars. However, for example, Wood et al. (2021) show that stars within our mass range in Table 2 have mass-loss rates of the order of the solar value, with a potential scatter of a few orders of magnitude. We find that the range of $\dot{M}$ shown in Table 2 is bracketed by the observational range suggested by Wood et al. (2021). Regarding the rotation rates, we chose a range to cover a parameter space that is well within the slowly rotating regime and also pushing into the rapidly rotating regime. Consequently, the range of rotation periods in Table 2 for the solar case includes the range typically observed (see, e.g., Bouvier et al. 2014, pp. 433–450) but it also extends to much longer periods than what is applicable for solar-like stars. On the other

hand, for the low-mass star cases, our parameter space includes the relatively fast-rotation side of the observed rotation rate distributions, but it does not extend to rotational periods as long as observed in the oldest stars (Bouvier et al. 2014). However, stars with longer rotation periods than those covered by our parameter space should still be compatible with our torque formulations because they are in they are in the slowly rotating regime, where changes in angular momentum are linear with stellar rotation rate.

#### 3.1. Morphology of Stellar-wind Systems

In this section, we demonstrate the qualitative properties of these stellar-wind simulations, focusing on the effects differential rotation have on the dynamical behavior of the stellar wind. Figure 1 shows the poloidal velocity distribution for the computational domains of three stellar-wind simulations: model 46 ($f_{eq} = 0.05$, $\alpha = 0$), model 14 ($f_{eq} = 0.001$, $\alpha = 0$), and model 34 ($f_{eq} = 0.05$, $\alpha = 0.75$). These are taken at simulation times corresponding to the middle of their respective time-median domains. The figures show the magnetic field lines (white), the ejection of stellar wind (green) along open field lines, the sonic Mach surface (black), and the Alfvén Mach surface (blue). For the rapid solid-body rotator (model 46), centrifugal effects become increasingly noticeable in the stellar wind, enhancing wind speeds at mid-to-high colatitudes and causing the wind to self-collimate and bend toward the poles; this results in an elongated Alfvén surface (see, e.g., Sakurai 1985; Washimi & Shibata 1993; Matt & Balick 2004). For the slow solid-body rotator (model 14), centrifugal effects are negligible, therefore the dynamics of the wind are predominantly determined by the stellar field strength and the mass-loss rate in the wind; this results in a more spherical Alfvén surface. For the differential rotator (model 34), even though the equator rotates as quickly as the rapid solid-body case, the rotation of the wind itself (originating from low colatitudes) is much slower, reducing centrifugal effects and resulting in dynamical behavior similar to that of the slow solid-body case.

In addition, we probe differences in the rotation of the wind between solid-body and differential rotation models, by calculating the effective rotation rate across the computational domain:

$$\Omega_{\text{eff}} = \frac{1}{R \sin\theta} \left( v_{\phi} - \frac{B_{\phi} \kappa}{\rho} \right). \quad (10)$$

In order to calculate accurate stellar torque values, this quantity is required to be conserved within some tolerance along magnetic field lines (as shown in Zanni & Ferreira 2009). In practice, we check its convergence by calculating $\Omega_{\text{eff}}$ in each grid cell and tracing the field line at this point back to its footprint on the stellar surface. We then compare this quantity to the stellar surface rotation rate at that given footprint’s colatitude, $\theta_f$, i.e., $\Omega_*(\theta_f)$; full conservation of this quantity is satisfied when $\Omega_{\text{eff}}(R, \theta) = \Omega_*(\theta_f)$. In Figure 2, we plot $\Omega_{\text{eff}}$ (normalized by the equatorial spin rate, $\Omega_{*,eq}$) as a function of the angle on the stellar surface between the footprint colatitude and its closest pole, for the same rapid solid-body (model 46) and differential rotation (model 34) cases from Figure 1. We plot the differential rotation profile (expressed in Equation (7)) as a function of $\theta_f$ (normalized by $\Omega_{*,eq}$) for comparison, where the

<footnote n="4">⁴ This solar mass-loss rate is the mean of 27 day averages via measurements of the solar-wind speed and density from the Advanced Composition Explorer (ACE) spacecraft (see Finley et al. 2018).</footnote>

<footnote n="5">⁵ See et al. (2019b) give a range of dipolar field strengths averaged over the stellar surface, $\langle B_d \rangle$; therefore, values of $B_*$ are likely to be within a factor of 2.</footnote>
<table><caption>**Table 1**
Variable Input Parameters and Outputed Global Variables for All Simulations</caption><thead><tr><th>Model</th><th>f<sub>eq</sub></th><th>α</th><th>v<sub>A</sub>/v<sub>esc</sub></th><th>B<sub>*</sub>a/B<sub>0</sub></th><th>M&#775;/M<sub>0</sub>(10<sup>-2</sup>)</th><th>j/J<sub>0</sub>(10<sup>-2</sup>)</th><th>Φ<sub>open</sub>/Φ<sub>0</sub></th><th>⟨r<sub>A,eq</sub>⟩/R<sub>*</sub></th><th>Υ<sub>*</sub>(10<sup>4</sup>)</th><th>ω<sub>emp</sub></th><th>Model</th><th>f<sub>eq</sub></th><th>α</th><th>v<sub>A</sub>/v<sub>esc</sub></th><th>B<sub>*</sub>a/B<sub>0</sub></th><th>M&#775;/M<sub>0</sub>(10<sup>-2</sup>)</th><th>j/J<sub>0</sub>(10<sup>-2</sup>)</th><th>Φ<sub>open</sub>/Φ<sub>0</sub></th><th>⟨r<sub>A,eq</sub>⟩/R<sub>*</sub></th><th>Υ<sub>*</sub>(10<sup>4</sup>)</th><th>ω<sub>emp</sub></th></tr></thead><tbody><tr><td>1</td><td>0.001</td><td>0</td><td>1</td><td>1.41</td><td>2.15</td><td>0.102</td><td>2.87</td><td>6.89</td><td>0.26</td><td>0.946</td><td>38</td><td>0.05</td><td>0.25</td><td>2</td><td>2.83</td><td>2.1</td><td>8.95</td><td>5.2</td><td>9.23</td><td>1.06</td><td>0.827</td></tr><tr><td>2</td><td>0.001</td><td>0</td><td>2</td><td>2.83</td><td>1.67</td><td>0.184</td><td>4.44</td><td>10.5</td><td>1.34</td><td>0.943</td><td>39</td><td>0.05</td><td>0.25</td><td>4</td><td>5.66</td><td>1.69</td><td>17.1</td><td>8.3</td><td>14.2</td><td>5.3</td><td>0.859</td></tr><tr><td>3</td><td>0.001</td><td>0</td><td>4</td><td>5.66</td><td>1.21</td><td>0.321</td><td>6.44</td><td>16.3</td><td>7.41</td><td>0.946</td><td>40</td><td>0.05</td><td>0.25</td><td>8</td><td>11.3</td><td>1.43</td><td>32.3</td><td>13.1</td><td>21.3</td><td>25</td><td>0.865</td></tr><tr><td>4</td><td>0.001</td><td>0</td><td>8</td><td>11.3</td><td>0.886</td><td>0.628</td><td>9.72</td><td>26.6</td><td>40.3</td><td>1.05</td><td>41</td><td>0.05</td><td>0.5</td><td>1</td><td>1.41</td><td>2.37</td><td>3.36</td><td>3.07</td><td>5.33</td><td>0.236</td><td>0.598</td></tr><tr><td>5</td><td>0.001</td><td>0.25</td><td>1</td><td>1.41</td><td>2.15</td><td>0.0816</td><td>2.87</td><td>6.16</td><td>0.26</td><td>0.757</td><td>42</td><td>0.05</td><td>0.5</td><td>2</td><td>2.83</td><td>1.95</td><td>6.36</td><td>4.97</td><td>8.07</td><td>1.14</td><td>0.607</td></tr><tr><td>6</td><td>0.001</td><td>0.25</td><td>2</td><td>2.83</td><td>1.67</td><td>0.145</td><td>4.44</td><td>9.32</td><td>1.34</td><td>0.744</td><td>43</td><td>0.05</td><td>0.5</td><td>4</td><td>5.66</td><td>1.49</td><td>11.5</td><td>7.68</td><td>12.4</td><td>5.99</td><td>0.609</td></tr><tr><td>7</td><td>0.001</td><td>0.25</td><td>4</td><td>5.66</td><td>1.21</td><td>0.25</td><td>6.44</td><td>14.4</td><td>7.41</td><td>0.736</td><td>44</td><td>0.05</td><td>0.5</td><td>8</td><td>11.3</td><td>1.23</td><td>22.4</td><td>12.2</td><td>19.1</td><td>29.1</td><td>0.642</td></tr><tr><td>8</td><td>0.001</td><td>0.25</td><td>8</td><td>11.3</td><td>0.886</td><td>0.484</td><td>9.74</td><td>23.4</td><td>40.3</td><td>0.808</td><td>45</td><td>0.05</td><td>0.75</td><td>1</td><td>1.41</td><td>2.27</td><td>2.2</td><td>2.98</td><td>4.41</td><td>0.246</td><td>0.399</td></tr><tr><td>9</td><td>0.001</td><td>0.5</td><td>1</td><td>1.41</td><td>2.15</td><td>0.0612</td><td>2.87</td><td>5.34</td><td>0.26</td><td>0.568</td><td>46</td><td>0.05</td><td>0.75</td><td>2</td><td>2.83</td><td>1.82</td><td>3.84</td><td>4.74</td><td>6.49</td><td>1.23</td><td>0.378</td></tr><tr><td>10</td><td>0.001</td><td>0.5</td><td>2</td><td>2.83</td><td>1.67</td><td>0.106</td><td>4.44</td><td>7.97</td><td>1.34</td><td>0.544</td><td>47</td><td>0.05</td><td>0.75</td><td>4</td><td>5.66</td><td>1.35</td><td>6.47</td><td>7.15</td><td>9.78</td><td>6.62</td><td>0.36</td></tr><tr><td>11</td><td>0.001</td><td>0.5</td><td>4</td><td>5.66</td><td>1.21</td><td>0.178</td><td>6.43</td><td>12.2</td><td>7.41</td><td>0.526</td><td>48</td><td>0.05</td><td>0.75</td><td>8</td><td>11.3</td><td>1.05</td><td>11.9</td><td>11</td><td>15.1</td><td>34.2</td><td>0.366</td></tr><tr><td>12</td><td>0.001</td><td>0.5</td><td>8</td><td>11.3</td><td>0.886</td><td>0.34</td><td>9.75</td><td>19.6</td><td>40.3</td><td>0.568</td><td>49</td><td>0.1</td><td>0</td><td>1</td><td>1.41</td><td>4.12</td><td>14.6</td><td>4.06</td><td>5.95</td><td>0.135</td><td>1.02</td></tr><tr><td>13</td><td>0.001</td><td>0.75</td><td>1</td><td>1.41</td><td>2.15</td><td>0.0409</td><td>2.87</td><td>4.36</td><td>0.26</td><td>0.38</td><td>50</td><td>0.1</td><td>0</td><td>2</td><td>2.83</td><td>3.75</td><td>29.3</td><td>6.95</td><td>8.84</td><td>0.596</td><td>1.05</td></tr><tr><td>14</td><td>0.001</td><td>0.75</td><td>2</td><td>2.83</td><td>1.67</td><td>0.0671</td><td>4.43</td><td>6.34</td><td>1.34</td><td>0.344</td><td>51</td><td>0.1</td><td>0</td><td>4</td><td>5.66</td><td>3.13</td><td>53.6</td><td>11.1</td><td>13.1</td><td>2.85</td><td>1.03</td></tr><tr><td>15</td><td>0.001</td><td>0.75</td><td>4</td><td>5.66</td><td>1.21</td><td>0.107</td><td>6.44</td><td>9.44</td><td>7.41</td><td>0.315</td><td>52</td><td>0.1</td><td>0</td><td>8</td><td>11.3</td><td>2.57</td><td>90.3</td><td>17</td><td>18.6</td><td>13.9</td><td>0.912</td></tr><tr><td>16</td><td>0.001</td><td>0.75</td><td>8</td><td>11.3</td><td>0.885</td><td>0.196</td><td>9.74</td><td>14.9</td><td>40.4</td><td>0.327</td><td>53</td><td>0.1</td><td>0.25</td><td>1</td><td>1.41</td><td>3.55</td><td>11.5</td><td>3.85</td><td>5.68</td><td>0.157</td><td>0.856</td></tr><tr><td>17</td><td>0.01</td><td>0</td><td>1</td><td>1.41</td><td>2.17</td><td>1.03</td><td>2.88</td><td>6.88</td><td>0.257</td><td>0.948</td><td>54</td><td>0.1</td><td>0.25</td><td>2</td><td>2.83</td><td>3.21</td><td>23</td><td>6.57</td><td>8.48</td><td>0.696</td><td>0.887</td></tr><tr><td>18</td><td>0.01</td><td>0</td><td>2</td><td>2.83</td><td>1.71</td><td>1.91</td><td>4.53</td><td>10.6</td><td>1.3</td><td>0.969</td><td>55</td><td>0.1</td><td>0.25</td><td>4</td><td>5.66</td><td>2.59</td><td>40.5</td><td>10.4</td><td>12.5</td><td>3.46</td><td>0.846</td></tr><tr><td>19</td><td>0.01</td><td>0</td><td>4</td><td>5.66</td><td>1.24</td><td>3.31</td><td>6.58</td><td>16.4</td><td>7.22</td><td>0.961</td><td>56</td><td>0.1</td><td>0.25</td><td>8</td><td>11.3</td><td>2.28</td><td>74.4</td><td>16.5</td><td>18.1</td><td>15.7</td><td>0.811</td></tr><tr><td>20</td><td>0.01</td><td>0</td><td>8</td><td>11.3</td><td>0.941</td><td>6.67</td><td>10.2</td><td>26.6</td><td>38</td><td>1.08</td><td>57</td><td>0.1</td><td>0.5</td><td>1</td><td>1.41</td><td>3.04</td><td>8.23</td><td>3.59</td><td>5.2</td><td>0.184</td><td>0.656</td></tr><tr><td>21</td><td>0.01</td><td>0.25</td><td>1</td><td>1.41</td><td>2.16</td><td>0.821</td><td>2.88</td><td>6.16</td><td>0.258</td><td>0.76</td><td>58</td><td>0.1</td><td>0.5</td><td>2</td><td>2.83</td><td>2.66</td><td>16.1</td><td>6.02</td><td>7.79</td><td>0.839</td><td>0.672</td></tr><tr><td>22</td><td>0.01</td><td>0.25</td><td>2</td><td>2.83</td><td>1.71</td><td>1.5</td><td>4.52</td><td>9.39</td><td>1.31</td><td>0.763</td><td>59</td><td>0.1</td><td>0.5</td><td>4</td><td>5.66</td><td>2.12</td><td>28.4</td><td>9.43</td><td>11.6</td><td>4.21</td><td>0.643</td></tr><tr><td>23</td><td>0.01</td><td>0.25</td><td>4</td><td>5.66</td><td>1.24</td><td>2.65</td><td>6.7</td><td>14.7</td><td>7.22</td><td>0.773</td><td>60</td><td>0.1</td><td>0.5</td><td>8</td><td>11.3</td><td>1.82</td><td>52.6</td><td>15.2</td><td>17</td><td>19.6</td><td>0.63</td></tr><tr><td>24</td><td>0.01</td><td>0.25</td><td>8</td><td>11.3</td><td>0.915</td><td>4.94</td><td>10</td><td>23.4</td><td>39</td><td>0.82</td><td>61</td><td>0.1</td><td>0.75</td><td>1</td><td>1.41</td><td>2.63</td><td>5.25</td><td>3.33</td><td>4.46</td><td>0.213</td><td>0.444</td></tr><tr><td>25</td><td>0.01</td><td>0.5</td><td>1</td><td>1.41</td><td>2.16</td><td>0.615</td><td>2.87</td><td>5.34</td><td>0.259</td><td>0.57</td><td>62</td><td>0.1</td><td>0.75</td><td>2</td><td>2.83</td><td>2.2</td><td>9.54</td><td>5.41</td><td>6.58</td><td>1.01</td><td>0.431</td></tr><tr><td>26</td><td>0.01</td><td>0.5</td><td>2</td><td>2.83</td><td>1.7</td><td>1.1</td><td>4.5</td><td>8.04</td><td>1.32</td><td>0.558</td><td>63</td><td>0.1</td><td>0.75</td><td>4</td><td>5.66</td><td>1.73</td><td>16.9</td><td>8.53</td><td>9.87</td><td>5.15</td><td>0.419</td></tr><tr><td>27</td><td>0.01</td><td>0.5</td><td>4</td><td>5.66</td><td>1.23</td><td>1.89</td><td>6.67</td><td>12.4</td><td>7.27</td><td>0.551</td><td>64</td><td>0.1</td><td>0.75</td><td>8</td><td>11.3</td><td>1.39</td><td>28</td><td>13</td><td>14.2</td><td>25.6</td><td>0.376</td></tr><tr><td>28</td><td>0.01</td><td>0.5</td><td>8</td><td>11.3</td><td>0.905</td><td>3.45</td><td>9.89</td><td>19.6</td><td>39.5</td><td>0.574</td><td>65</td><td>0.2</td><td>0</td><td>1</td><td>1.41</td><td>10.4</td><td>42.4</td><td>5.52</td><td>4.52</td><td>0.0538</td><td>1.04</td></tr><tr><td>29</td><td>0.01</td><td>0.75</td><td>1</td><td>1.41</td><td>2.15</td><td>0.411</td><td>2.87</td><td>4.37</td><td>0.259</td><td>0.381</td><td>66</td><td>0.2</td><td>0</td><td>2</td><td>2.83</td><td>9.77</td><td>85.9</td><td>9.63</td><td>6.63</td><td>0.229</td><td>1.07</td></tr><tr><td>30</td><td>0.01</td><td>0.75</td><td>2</td><td>2.83</td><td>1.68</td><td>0.681</td><td>4.46</td><td>6.37</td><td>1.33</td><td>0.348</td><td>67</td><td>0.2</td><td>0</td><td>4</td><td>5.66</td><td>7.97</td><td>147</td><td>15</td><td>9.64</td><td>1.12</td><td>0.978</td></tr><tr><td>31</td><td>0.01</td><td>0.75</td><td>4</td><td>5.66</td><td>1.21</td><td>1.11</td><td>6.56</td><td>9.57</td><td>7.39</td><td>0.325</td><td>68</td><td>0.2</td><td>0</td><td>8</td><td>11.3</td><td>6.26</td><td>238</td><td>22.9</td><td>13.9</td><td>5.71</td><td>0.853</td></tr><tr><td>32</td><td>0.01</td><td>0.75</td><td>8</td><td>11.3</td><td>0.892</td><td>1.82</td><td>9.25</td><td>14.3</td><td>40</td><td>0.302</td><td>69</td><td>0.2</td><td>0.5</td><td>1</td><td>1.41</td><td>6.11</td><td>24.9</td><td>4.89</td><td>4.51</td><td>0.0914</td><td>0.748</td></tr><tr><td>33</td><td>0.05</td><td>0</td><td>1</td><td>1.41</td><td>2.66</td><td>5.85</td><td>3.28</td><td>6.62</td><td>0.21</td><td>0.986</td><td>70</td><td>0.2</td><td>0.5</td><td>2</td><td>2.83</td><td>5.5</td><td>47.7</td><td>8.29</td><td>6.59</td><td>0.406</td><td>0.739</td></tr><tr><td>34</td><td>0.05</td><td>0</td><td>2</td><td>2.83</td><td>2.29</td><td>11.8</td><td>5.47</td><td>10.1</td><td>0.974</td><td>1.05</td><td>71</td><td>0.2</td><td>0.5</td><td>4</td><td>5.66</td><td>4.2</td><td>75.9</td><td>12.7</td><td>9.53</td><td>2.13</td><td>0.649</td></tr><tr><td>35</td><td>0.05</td><td>0</td><td>4</td><td>5.66</td><td>1.87</td><td>22.7</td><td>8.85</td><td>15.5</td><td>4.76</td><td>1.09</td><td>72</td><td>0.25</td><td>0</td><td>2</td><td>2.83</td><td>14.8</td><td>126</td><td>10.8</td><td>5.83</td><td>0.151</td><td>1.08</td></tr><tr><td>36</td><td>0.05</td><td>0</td><td>8</td><td>11.3</td><td>1.61</td><td>42.6</td><td>14.1</td><td>23</td><td>22.2</td><td>1.08</td><td>73</td><td>0.25</td><td>0</td><td>4</td><td>5.66</td><td>11.9</td><td>211</td><td>16.9</td><td>8.46</td><td>0.749</td><td>0.967</td></tr><tr><td>37</td><td>0.05</td><td>0.25</td><td>1</td><td>1.41</td><td>2.5</td><td>4.58</td><td>3.17</td><td>6.05</td><td>0.223</td><td>0.794</td><td>74</td><td>0.3</td><td>0</td><td>4</td><td>5.66</td><td>17.8</td><td>299</td><td>18.7</td><td>7.54</td><td>0.502</td><td>1</td></tr></tbody></table>

**Note.**

<footnote>ª B<sub>*</sub> is not a fundamental input parameter, but it is simply derived from v<sub>A</sub>/v<sub>esc</sub> and tabulated here, for convenience.</footnote>
<table><caption>**Table 2**
Example Ranges of Global Properties in Physical Units (for Our Parameter Space) for $M_{\star} = 0.1 - 1 M_{\odot}$ (5 Gyr)</caption><thead><tr><th>M<sub>\*</sub><br>(M<sub>☉</sub>)</th><th>R<sub>\*</sub><br>(R<sub>☉</sub>)</th><th>v<sub>K,\*</sub><br>(km s<sup>-1</sup>)</th><th>ρ<sub>\*</sub><br>(g cm<sup>-3</sup>)</th><th>B<sub>\*</sub><br>(G)</th><th>P<sub>\*,eq</sub><br>(days)</th><th>M&#775;<sub>\*</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>j<br>(erg)</th><th>Φ<sub>open</sub><br>(Mx)</th></tr></thead><tbody><tr><td>0.10</td><td>0.117</td><td>404</td><td>2.46 × 10<sup>-13</sup></td><td>100 - 800</td><td>0.0489 - 14.7</td><td>9.23 × 10<sup>-14</sup> - 1.85 × 10<sup>-12</sup></td><td>8.84 × 10<sup>28</sup> - 6.46 × 10<sup>32</sup></td><td>1.35 × 10<sup>22</sup> - 1.08 × 10<sup>23</sup></td></tr><tr><td>0.25</td><td>0.255</td><td>432</td><td>1.57 × 10<sup>-14</sup></td><td>27.2 - 218</td><td>0.0995 - 29.8</td><td>3.02 × 10<sup>-14</sup> - 6.06 × 10<sup>-13</sup></td><td>6.75 × 10<sup>28</sup> - 4.93 × 10<sup>32</sup></td><td>1.74 × 10<sup>22</sup> - 1.39 × 10<sup>23</sup></td></tr><tr><td>0.50</td><td>0.455</td><td>458</td><td>1.97 × 10<sup>-15</sup></td><td>10.2 - 81.5</td><td>0.168 - 50.3</td><td>1.27 × 10<sup>-14</sup> - 2.55 × 10<sup>-13</sup></td><td>5.37 × 10<sup>28</sup> - 3.92 × 10<sup>32</sup></td><td>2.07 × 10<sup>22</sup> - 1.65 × 10<sup>23</sup></td></tr><tr><td>1.00</td><td>1.00</td><td>437</td><td>2.46 × 10<sup>-16</sup></td><td>3.44 - 27.5</td><td>0.386 - 116</td><td>7.32 × 10<sup>-15</sup> - 1.47 × 10<sup>-13</sup></td><td>6.48 × 10<sup>28</sup> - 4.73 × 10<sup>32</sup></td><td>3.37 × 10<sup>22</sup> - 2.69 × 10<sup>23</sup></td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Poloidal velocity color maps, showing (left) model 46 ($f_{\text{eq}} = 0.05$, $\alpha = 0$), (center) model 14 ($f_{\text{eq}} = 0.001$, $\alpha = 0$), and (right) model 34 ($f_{\text{eq}} = 0.05$, $\alpha = 0.75$). All simulations above have a surface polar field strength $B_*/B_0 = 2.83$. Magnetic field lines (white), velocity vectors normalized by $v_{K,*}$ (green), the sonic surface (black), and the Alfvén surface (blue) are included. We include a secondary color bar that gives poloidal velocity in physical units (km s⁻¹) for the 0.1–1 M☉ cases in Table 2. As $v_{K,*}$ for these cases only varies by a maximum of ≈7% from the mean value $\langle v_{K,*} \rangle = 432.75$ km s⁻¹, we adopt this to produce a single color bar that approximately represents the entire mass range.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** $\Omega_{\text{eff}}/\Omega_{*,\text{eq}}$ as a function of the angle on the stellar surface between the footprint colatitude and its closest pole, showing (left) model 46 ($f_{\text{eq}} = 0.05$, $\alpha = 0$) and (right) model 34 ($f_{\text{eq}} = 0.05$, $\alpha = 0.75$). Both simulations have a field strength $B_*/B_0 = 2.83$. The red dotted line shows the stellar differential rotation profile $\Omega_*(\theta_f)$ adopted at the inner boundary of the computational domain (equation (7)), normalized by $\Omega_{*,\text{eq}}$. The color map illustrates the divergence of $\Omega_{\text{eff}}$ from the stellar differential rotation profile $\Omega_*(\theta_f)$. The white (gray) area illustrates the region on the stellar surface in which the magnetic field is open (closed).</figcaption></figure>

color map illustrates the divergence of $\Omega_{\text{eff}}$ from this profile.<fnref n="6" /> In both cases, the white area represents the area on the stellar surface in which the stellar wind is being ejected, i.e., where the magnetic field lines are open. It is evident that the stellar-wind region decreases with increasing differential rotation. For solid-body rotation, $\Omega_{\text{eff}}$ appears to be well conserved in the region of

the stellar wind (where $\Omega_{\text{eff}}/\Omega_*(\theta_f) \approx 1$), and the rotation of the wind is roughly constant. For differential rotation, $\Omega_{\text{eff}}$ is also well conserved in the wind, but the rotation of the wind decreases with decreasing colatitude (toward the poles).

## 3.2. Global Flow Quantities and Efficiencies

The main aim of this paper is to investigate how the stellar wind interacts with a star that is undergoing differential rotation. By integrating over a spherical shell of surface S

<footnote n="6">⁶ Divergence between $\Omega_{\text{eff}}$ and $\Omega_*(\theta_f)$ typically occurs at the transition between open and closed field lines, due to increasing numerical diffusion.</footnote>
perpendicular to the poloidal flow at a given radius R, we express the radial profiles of the stellar-wind mass flow rate, torque, and unsigned magnetic flux as

$$ \dot{M} = \int_S \rho v_p \cdot dS, \qquad (11) $$

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

and

$$ \Phi = \int_S |\mathbf{B} \cdot d\mathbf{S}|, \qquad (13) $$

respectively, where

$$ \Lambda = R \sin \theta \left( v_{\phi} - \frac{B_{\phi}}{\kappa} \right) \qquad (14) $$

is a quantity related to the angular momentum flux $F_{AM} = \Lambda \rho v$. In practice, we determine global values of $\dot{M}$ and $\dot{j}$ for all simulations by taking the median values over all spherical shells at $R > 7.68R_*$ (corresponding to the outer half of the logarithmic radial grid) in order to avoid numerical effects near the inner boundary. The unsigned magnetic flux for a system with dipolar topology initially falls as $1/R$, but reaches a constant value that represents the open flux in the stellar wind, $\Phi_{open}$, outside the Alfvén surface, i.e., when the thermal and hydrodynamic wind pressure exceeds the magnetic pressure. In practice, we determine global values of $\Phi_{open}$ for all simulations by taking the median values over all spherical shells at $R > 40R_*$. Values of $\dot{M}$, $\dot{j}$, and $\Phi_{open}$ for all simulations can be found in Table 1.

We define a “wind magnetization” parameter that encapsulates the magnetic and mass-loss stellar-wind properties, based on properties at the stellar surface:

$$ \Upsilon_* = \frac{\Phi_*^2}{4\pi R_*^2 \dot{M} v_{\text{esc}}} \qquad (15) $$

where $\Phi_* = \alpha \pi R_*^2 B_*$ is the total stellar flux ($\alpha = 2$ for our dipolar surface configuration). Values of $\Upsilon_*$ for all simulations can be found in Table 1. In Figure 3, we illustrate the extent of the parameter space explored by our simulations, in terms of $f_{eq}$ and $\Upsilon_*$. Symbol shapes, colors, and borders demonstrate how each parameter in our study varies. Increasing $B_*$ noticeably increases $\Upsilon_*$ (at fixed $f_{eq}$ and $\alpha$). For slower rotation rates ($f_{eq} \le 0.01$), changes in $f_{eq}$ and/or $\alpha$ have a negligible effect on $\Upsilon_*$ (at fixed $B_*$); therefore, simulations with changing $\alpha$ (for fixed $B_*$ and $f_{eq}$) are seen to overlap in this figure. For more rapid rotation ($f_{eq} \ge 0.05$), increasing $f_{eq}$ and/or decreasing $\alpha$ decreases $\Upsilon_*$ (at fixed $B_*$).

Equation (15) demonstrates that these simulations are degenerate in $B_*^2/\dot{M}$. This means that for a real star, the range of $\Upsilon_*$ shown in Figure 3 can be taken to represent a range of $B_*$ for some fixed value of $\dot{M}$, a range of $\dot{M}$ for some fixed $B_*$, or for any combination of the two. Table 2 shows the range of $B_*$ and $\dot{M}$ in physical units for a few representative MS star masses and choices of coronal-density normalization. Similarly, the range of $f_{eq}$ shown in Figure 3 represents a range of equatorial rotation periods, where example physical values are given in Table 2.

<figure><img src="image_6.png" /><figcaption>Figure 3. $f_{eq}$ as a function of $\Upsilon_*$ for all simulations. To differentiate between different input parameters, we use differing symbol styles for $f_{eq}$, marker colors for $\alpha$, and marker edge colors for $B_*$.</figcaption></figure>

## 4. Torque Formulation for Simulations of Stellar Winds

Through dimensional analysis, a general parameterization for the stellar-wind torque can be written as

$$ \dot{j} = \dot{M} \langle \Omega_{\text{wind}} \rangle \langle r_A \rangle^2, \qquad (16) $$

where $\langle \Omega_{\text{wind}} \rangle$ represents an “effective” stellar rotation rate in the stellar-wind region, and $\langle r_A \rangle$ is a characteristic length scale that represents the “effective magnetic lever arm” and quantifies the efficiency of the stellar-wind torque (see, e.g., Weber & Davis 1967). The Alfvén radius (normalized by $R_*$) is determined by rearranging Equation (16):

$$ \frac{\langle r_A \rangle}{R_*} = \left( \frac{\dot{j}}{\dot{M} \langle \Omega_{\text{wind}} \rangle R_*^2} \right)^{1/2} \qquad (17) $$

(following Washimi & Shibata 1993; Matt & Pudritz 2008). However, $\langle \Omega_{\text{wind}} \rangle$ is not known a priori, so we instead normalize the angular momentum loss in the stellar wind using the equatorial rotation rate $\Omega_{\star,eq}$, giving

$$ \frac{\langle r_{A,eq} \rangle}{R_*} = \left( \frac{\dot{j}}{\dot{M} \Omega_{\star,eq} R_*^2} \right)^{1/2} \equiv \omega^{1/2} \frac{\langle r_A \rangle}{R_*} \qquad (18) $$

where we define

$$ \omega = \frac{\langle \Omega_{\text{wind}} \rangle}{\Omega_{\star,eq}} \equiv \frac{\langle f_{\text{wind}} \rangle}{f_{\text{eq}}} \qquad (19) $$

where $\langle f_{\text{wind}} \rangle$ represents an “effective” stellar breakup fraction in the stellar wind. Values of $\langle r_{A,eq} \rangle/R_*$ for all simulations can be found in Table 1. Therefore, we rewrite Equation (16) in terms of simulation parameters:

$$ \dot{j} = \dot{M} \Omega_{\star,eq} \langle r_{A,eq} \rangle^2, \qquad (20) $$

We adopt the stellar-wind torque formulation developed by Matt et al. (2012):

$$ \frac{\langle r_A \rangle}{R_*} = K_{\star,1} \left( \frac{\Upsilon_*}{\beta_*} \right)^{m_*} \qquad (21) $$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** $\langle r_{A,eq} \rangle / R_*$ as a function of $\Upsilon_*/\beta_*$ (fixing $\omega = 1$) for all simulations. The black line shows the solid-body fit $\langle r_{A,eq} \rangle / R_* = 0.932(\Upsilon_*/\beta_*)^{0.258}$ (with $K_{*,2} = 0.263$).</figcaption></figure>

where $\Upsilon_*$ is the “magnetization parameter” of the stellar wind based on properties at the stellar surface (expressed in Equation (15)),

$$ \beta_* = \left[ 1 + \left( \frac{\omega f_{eq}}{K_{*,2}} \right)^2 \right]^{1/2} \qquad (22) $$

is a “magnetocentrifugal correction” term for the wind efficiency (expressed in the form adopted by Réville et al. 2015), and $K_{*,1}$, $K_{*,2}$, and $m_*$ are best-fit dimensionless parameters. The Alfvén radius is inversely proportional to the stellar-wind mass-loss rate. Using Equation (18), we rewrite Equation (21) in terms of $\langle r_{A,eq} \rangle / R_*$:

$$ \frac{\langle r_{A,eq} \rangle}{R_*} = K_{*,1} \left( \frac{\Upsilon_*}{\beta_*} \right)^{m_*} \omega^{1/2}. \qquad (23) $$

In Figure 4, we plot $\langle r_{A,eq} \rangle / R_*$ as a function of $\Upsilon_*/\beta_*$, (fixing $\omega = 1$ for the $\beta_*$ term). By solely fitting the solid-body rotating cases, we find the best-fit dimensionless parameters to be $K_{*,1} = 0.932$, $K_{*,2} = 0.263$, and $m_* = 0.258$. For convenience, all of the best-fit dimensionless parameters presented in Section 4 can be found in Table 3. The cases with solid-body rotation are well fit by the power law represented by Equation (21), illustrated by the small amount of scatter of those cases around the best-fit line. However, the cases with differential rotation exhibit a systematic deviation of the points toward smaller $\langle r_{A,eq} \rangle / R_*$ values for larger amounts of differential rotation (larger $\alpha$). It is clear that any stellar-wind torque formulation that assumes solid-body rotation (i.e., one that uses the equatorial rotation rate and assumes $\omega = 1$) will increasingly overestimate the efficiency of the stellar-wind torque for increasing solar-like differential rotation. We can define an empirical value for the normalized “effective” rotation rate in the wind, $\omega_{\text{emp}}$, as being the deviation of the simulated values from the simple power-law relationship of the solid-body-rotator case. That is, the “effective” rotation rate in the wind can be measured using Equation (23) and the fit

<table><caption>**Table 3**

Best-fit Dimensionless Parameters from Scaling Laws in Section 4</caption><thead><tr><th>Parameter</th><th>Value</th><th>Equations</th></tr></thead><tbody><tr><td>$K_{*,1}$</td><td>0.932</td><td>(21)</td></tr><tr><td>$K_{*,2}$</td><td>0.263</td><td></td></tr><tr><td>$m_*$</td><td>0.258</td><td></td></tr><tr><td>$K_\theta$</td><td>1.06</td><td>(31)</td></tr><tr><td>$K_{\Phi,1}$</td><td>1.16</td><td></td></tr><tr><td>$K_{\Phi,2}$</td><td>0.0111</td><td></td></tr><tr><td>$m_\phi$</td><td>-0.165</td><td></td></tr></tbody></table>

parameters ($K_{*,1}$, $K_{*,2}$, and $m_*$):

$$ \omega_{\text{emp}} = \left( \frac{\langle r_{A,eq} \rangle / R_*}{K_{*,1}(\Upsilon_*/\beta_*)^{m_*}} \right)^2. \qquad (24) $$

Equation (24) is a transcendental equation (due to the $\omega$ factor present in $\beta_*$), therefore, we solve this iteratively. Values of $\omega_{\text{emp}}$ for all simulations can be found in Table 1.

In order to predict the stellar-wind torque as a function of the stellar surface properties, we must adopt a formulation for $\omega$. The “effective” rotation rate in the wind, as defined by $\omega$, is a torque-averaged rotation rate over all open field lines carrying wind material. Because our adopted differential rotation profile is a simple function of colatitude on the stellar surface, we can define the colatitude on the star that happens to rotate at the same rate as the “effective” rotation of the wind, $\theta_{\text{eff}}$. Equations (7) and (19) thus define the “effective” rotation rate in terms of this colatitude:

$$ \omega_{\text{th}} = 1 - \alpha \cos^2 \theta_{\text{eff}}. \qquad (25) $$

Because the rotation rate increases from pole to equator, and because the integrated torque (and thus “effective” rotation in the wind) is dominated by flow from the highest colatitudes, we expect that $\theta_{\text{eff}}$ will have a similar value and similar scaling to the maximum colatitude on the stellar surface that emits a wind (as opposed to the higher colatitudes where there is a dead zone), which we call $\theta_{\text{open}}$. In practice, we explored a few different possible relationships between $\theta_{\text{eff}}$ and $\theta_{\text{open}}$. Given the functional form of Equation (25), we found that a relatively simple form and acceptable fit to the data is obtained by assuming that $\cos^2 \theta_{\text{eff}}$ is proportional to $\cos^2 \theta_{\text{open}}$, leading to

$$ \omega_{\text{th}} = 1 - \alpha K_{\theta} \cos^2 \theta_{\text{open}}, \qquad (26) $$

where $K_\theta$ is a best-fit dimensionless parameter. By evaluating Equation (13) at $R_*$ across the average area enclosing the stellar wind, i.e., for $0 \le \theta \le \theta_{\text{open}}$ and $(\pi - \theta_{\text{open}}) \le \theta \le \pi$, using the dipolar definition for $B_R$ (Equation (3)), one can write

$$ \cos^2 \theta_{\text{open}} = 1 - \frac{\Phi_{\text{open}}}{\Phi_*}, \qquad (27) $$

giving

$$ \omega_{\text{th}} = 1 - \alpha K_{\theta} \left( 1 - \frac{\Phi_{\text{open}}}{\Phi_{*}} \right) \qquad (28) $$

for our dipolar configuration. Therefore, when $\alpha \ne 0$, the “effective” rotation rate of the wind (normalized by the equatorial stellar rotation rate) is a function of $\alpha$ and the fraction of open flux in the stellar wind.
In previous stellar-wind simulations, it has been shown that the fraction of open flux increases with the wind mass-loss rate (at fixed $f_{eq}/B_\star$); more specifically, $\Phi_{open}/\Phi_\star$ scales inversely with $\Upsilon_\star$ (see, e.g., Réville et al. 2015; Finley & Matt 2017, 2018; Pantolmos & Matt 2017). Therefore, we parameterize $\Phi_{open}/\Phi_\star$ using the following functional fit:

$$ \frac{\Phi_{\text{open}}}{\Phi_{\star}} = K_{\Phi,1} \left( \frac{\Upsilon_{\star}}{\beta_{\Phi}} \right)^{m_{\Phi}}, \qquad (29) $$

where

$$ \beta_{\Phi} = \left\{ 1 + \left[ \frac{(1 - \alpha)f_{\text{eq}}}{K_{\Phi,2}} \right]^2 \right\}^{1/2} \qquad (30) $$

is the “magnetocentrugal correction” term for the fraction of open flux in the stellar wind, and $K_{\Phi,1}$, $K_{\Phi,2}$, and $m_\Phi$ are best-fit dimensionless parameters. By solely fitting the solid-body rotation cases, we find the best-fit dimensionless parameters to be $K_{\Phi,1} = 1.16$, $K_{\Phi,2} = 0.0111$, and $m_\Phi = -0.165$.<fnref n="7" /> In Equation (29), we assume $\omega \approx (1 - \alpha)$ in the $\beta_\Phi$ term for simplicity to avoid producing a transcendental equation for the full parameterization of $\omega$; this approximation gives simulation points that lie within $\lesssim 10\%$ of the formulation if $\omega$ was used. Therefore, one can rewrite Equation (28) as

$$ \omega_{\text{th}} = 1 - \alpha K_{\theta} \left[ 1 - K_{\Phi,1} \left( \frac{\Upsilon_{\star}}{\beta_{\Phi}} \right)^{m_{\Phi}} \right]. \qquad (31) $$

We determine the best-fit dimensionless parameter $K_\theta = 1.06$ by fitting the $\omega_{\text{emp}}$ values with the expression for $\omega_{\text{th}}$ (equation (31)). Conceptually, given the differential profile used in this study, $\theta_{\text{open}}$ corresponds to the colatitude of the field line with the maximum stellar rotation rate in the stellar wind. Our value of $K_\theta$ demonstrates that $\langle\Omega_{\text{wind}}\rangle$ acts along a field line connected at a colatitude that is near, but lower than, $\theta_{\text{open}}$. In other words, the location on the stellar surface that rotates at the same rate as the “effective” rotation rate of the wind is located within the wind-emitting region and near the last open field line, as expected.

In Figure 5, we plot $\omega_{\text{emp}}$ as a function of $\Upsilon_\star/\beta_\Phi$. When the differential rotation increases, i.e., as $\alpha$ increases, the “effective” rotation rate in the wind is smaller. For a given differential rotation profile, the trend with $\Upsilon_\star/\beta_\Phi$ illustrates that the “effective” rotation rate in the wind is sensitive to the fraction of open flux in the wind. Specifically, if there is more fractional open flux (which occurs for smaller $\Upsilon_\star/\beta_\Phi$), the “effective” rotation of the wind is faster because it includes more flow from higher colatitudes. We superimpose lines representing the fully parameterized $\omega_{\text{th}}$ for each $\alpha$ explored in our study (Equation (31)). The scatter present for $\alpha = 0$ cases is solely due to deviations in the relationship between $\langle r_{A,eq} \rangle/R_\star$ and $\Upsilon_\star/\beta_\star$, whereas for $\alpha > 0$ cases, the scatter is likely a combination of this and from our parameterization of $\omega$.

In Figure 6, we plot $(\langle r_{A,eq} \rangle/R_\star)\omega_{\text{th}}^{-1/2}$ as a function of $\Upsilon_\star/\beta_\star$ (using $\omega_{\text{th}}$ in the expression of $\beta_\star$). We find all simulation points to lie within $\lesssim 10\%$ of the fitting function. Therefore, our adoption (and parameterization) of $\omega$ appears to strongly account for the stellar differential rotation effects on

<figure><img src="image_8.png" /><figcaption>Figure 5. $\omega_{\text{emp}}$ as a function of $\Upsilon_\star/\beta_\Phi$ for all simulations. Dotted lines represent the theoretical solution $\omega_{\text{th}}$, for a given $\alpha$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6. $(\langle r_{A,eq} \rangle/R_\star)\omega_{\text{th}}^{-1/2}$ as a function of $\Upsilon_\star/\beta_\star$ for all simulations. The black line shows the fit $(\langle r_{A,eq} \rangle/R_\star)\omega_{\text{th}}^{-1/2} = 0.932(\Upsilon_\star/\beta_\star)^{0.258}$ (with $K_{\star,2} = 0.263$).</figcaption></figure>

the efficiency of the stellar-wind torque. Finally, by combining Equations (20), (23), and (31), we express a fully parameterized estimate for the stellar-wind torque, for a wide range of differential rotation rates and stellar surface field strengths:

$$ j = K_{\star,1}^2 \dot{M} \Omega_{\star,eq} R_{\star}^2 \left( \frac{\Upsilon_{\star}}{\beta_{\star}} \right)^{2m_{\star}} \times \left\{ 1 - \alpha K_{\theta} \left[ 1 - K_{\Phi,1} \left( \frac{\Upsilon_{\star}}{\beta_{\Phi}} \right)^{m_{\Phi}} \right] \right\}. \qquad (32) $$

## 5. Discussion and Conclusions

In this work, we perform 74 2.5D axisymmetric MHD stellar-wind simulations to investigate how the stellar-wind torque experienced by the star is impacted by differential rotation and the dipolar stellar surface magnetic field strength. We create a torque formulation that encapsulates these effects, allowing for the investigation of the rotational evolution of MS stars using 1D stellar evolution codes. We chose a simplified setup, including a purely dipolar magnetic field, thermal/

<footnote n="7">⁷ Note, $K_{\Phi,2} \neq K_{\star,2}$. The $K_{\star,2}$ term captures the combined effects of the wind acceleration and the fact that the open flux is changing, whereas $K_{\Phi,2}$ captures only the latter.</footnote>
polytropic wind driving, and a monotonic differential rotation profile.

## 5.1. How Differential Rotation Affects Stellar-wind Torques

We determine a stellar-wind torque formulation through the parameterization of the effective Alfvén radius, based on the formulation employed by Matt et al. (2012). Our solid-body rotation simulations follow a $\langle r_A \rangle - \Upsilon_* - f$ relationship (Equation (21)) with a power-law index of $m_* = 0.258$, which is comparable to that found in previous stellar-wind studies (with an initial dipolar configuration) within systematic uncertainties (see, e.g., Matt et al. 2012; Réville et al. 2015; Finley & Matt 2017; Pantolmos & Matt 2017; Pantolmos et al. 2020). In addition, our magnetocentrifugal correction term for the wind efficiency, $\beta_*$, has a large fitting constant of $K_{*,2} = 0.263$, suggesting that rotation has an impact on the wind acceleration only for our most-rapid cases. This value is slightly larger than that determined by Réville et al. (2015) (~0.2), but much larger than that determined by Matt et al. (2012) (0.07); as discussed in the former, these discrepancies might be due to differences in boundary conditions.

For cases undergoing surface differential rotation, the effective Alfvén radius now depends on the “effective” stellar rotation rate in the wind, $\langle\Omega_{\text{wind}}\rangle$, which is not known a priori (Equation (17)). We renormalize the Alfvén radius using the equatorial stellar rotation rate, $\Omega_{*,eq}$, introducing the term $\omega = \langle\Omega_{\text{wind}}\rangle/\Omega_{*,eq}$ into the formulation (Equation (18)). We find that $\omega$ decreases with the magnitude of surface differential rotation (measured via $\alpha$). Furthermore, $\omega$ increases with the fraction of open flux in the wind (for a given $\alpha$); this corresponds to the stellar wind occupying a larger area of the domain, allowing for the inclusion of faster flow from higher colatitudes. The fraction of open flux in the wind itself inversely scales with the wind magnetization, $\Upsilon_*$; for rapidly rotating cases, it is also sensitive to magnetocentrifugal effects on the wind acceleration and the open flux combined, $\beta_\Phi$. Our simulations follow a $\Phi_{\text{open}}/\Phi_* - \Upsilon_* - f$ relationship with a power-law index of $m_\Phi = -0.165$, which is again comparable to previous open flux stellar-wind studies (see, e.g., Réville et al. 2015; Finley & Matt 2017; Pantolmos & Matt 2017; Pantolmos et al. 2020). Our $\beta_\Phi$ term has a proportionality constant of $K_{\Phi,2} = 0.0111$, suggesting that the open flux is more sensitive to rotation compared to the effects on wind acceleration solely (encapsulated by our $\beta_*$ term).

Overall, for our solar-like differential rotation cases, i.e., $\omega < 1$, we find the rotation in the wind to be decreased compared to a solid-body case with identical stellar properties (rotating at $\Omega_{*,eq}$ with a fixed mass-loss rate), decreasing the stellar-wind torque by approximately a factor of $\omega$ (ignoring centrifugal effects).<fnref n="8" /> We note that the parameterization of $\omega$ used in our formulation is specific to the differential rotation profile adopted; however, it is likely that one can derive $\omega$ for different profiles and follow the framework of this paper to predict the torque.

Differential rotation, in general, can shear small-scale magnetic field structures and affect the coronal dynamics (see, e.g., Lionello et al. 2006; Morgan 2011; Lionello et al. 2020). However, in this study, we solely investigate the effects

of differential rotation on the rotation rate at the base of the wind, and how this affects global angular momentum loss, independent of coronal dynamics. In this case, we are concerned only about the rotation rate of open magnetic field lines (i.e., single footprints), rather than connected footprints undergoing shear. In a 3D domain, differential rotation can redistribute the open magnetic flux of a nonaxisymmetric field (e.g., a tilted dipole); hence, the assumed initial magnetic flux distribution changes in time and the solution cannot be stationary. On the other hand, the timescales to establish a wind solution are much shorter than the differential rotation and magnetic field reorganization timescales. Therefore, our axisymmetric results represent the solution for a wind in which the timescale of reorganization is much longer than the time to establish the (quasi-)steady-state conditions in the wind. We also note that assuming a dipolar magnetic configuration is a simplification, and real solar-like and low-mass stars exhibit smaller-scale magnetic field structures, in addition to their dipolar components (e.g., Donati & Landstreet 2009). However, various works have suggested that the open flux, hence the angular momentum loss, is dominated by the dipolar magnetic field component (Jardine et al. 2016; See et al. 2016, 2017; Finley & Matt 2018). Our simulation results should apply to any real star, under the condition that the differential rotation of the open field lines is similar to the profiles we explored. Future studies should consider more realistic magnetic field geometries and differential rotation profiles, in order to explore whether magnetic topology, closed-loop shearing, and other coronal dynamics affect the differential rotation in the wind.

## 5.2. Variations during a Solar Magnetic Cycle

We use our modified prescription to illustrate the potential effect differential rotation could have on predicting the wind torque of the Sun, using time-varying solar properties (over solar cycles 23 to 24) calculated in Finley et al. (2018). Dipolar magnetic field strength is calculated using synoptic magnetograms of the surface field strength (from solar cycles 23 to 24) from two sources: the Michelson Doppler Imager at the Solar and Heliospheric Observatory (SOHO/MDI), and the Helioseismic and Magnetic Imager at the Solar Dynamic Observatory (SDO/HMI). The estimated dipolar field strength (at the pole) varies in the range 0.108–3.14 G (average: 1.54 G). The mass-loss rate is calculated from 27 day averages via measurements of the solar-wind speed and density from the ACE spacecraft, i.e., $\dot{M}_\odot = 4\pi \langle R^2 v_R(R) \rho(R) \rangle_{27\text{day}}$. The estimated mass-loss rate varies in the range (0.697 – 4.30) $\times 10^{-14} \text{ M}_\odot \text{ yr}^{-1}$ (average: $1.78 \times 10^{-14} \text{ M}_\odot \text{ yr}^{-1}$). A typical measured value of the solar equatorial sidereal rotation period is 24.47 days (Snodgrass & Ulrich 1990); this corresponds to a rotation rate of $\Omega_{\odot,eq} = 2.97 \times 10^{-6} \text{ rad s}^{-1}$, giving $f_{\odot,eq} = 4.73 \times 10^{-3}$. We adopt a relative differential rotation rate of $\alpha_\odot = 0.3$.

Using these observations, we determine corresponding values of $\omega_\odot$ (via Equation (31)) as a function of time. In the top panel of Figure 7, we illustrate the “effective” rotation of the wind’s sensitivity to the stellar-wind area by plotting the colatitude of the opening field line (determined with $\omega_\odot$ via Equation (26)) as a function of time for solar cycles 23 and 24 (black line). Gray lines are interpolated values, due to missing data for these periods. We superimpose a color map of $\omega_\odot(\theta_{\text{open},\odot} = \theta)$ (determined by substituting $\theta$ into Equation (26)), which demonstrates the distribution of its

<footnote n="8">⁸ For fixed mass-loss rate, one can show $\dot{J}/\dot{J}_{\text{SB}} = (\beta_{*,\text{SB}}/\beta_*)^{2m_*} \omega$, where the subscript “SB” represents the solid-body value. For $f_{\text{eq}} \ll K_{*,2}$, i.e., when centrifugal effects are negligible, $\dot{J}/\dot{J}_{\text{SB}} \approx \omega$.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Estimated opening colatitude of the solar wind (top), and estimated percentage change in predicted solar torque (bottom), as a function of time (in years) over solar cycles 23 and 24.</figcaption></figure>

theoretical value for a given colatitude of the opening field line. At the solar minima, where the stellar-wind region is at its smallest (corresponding to lower opening colatitudes), $\omega_{\odot} \approx 0.8$, i.e., the “effective” rotation of the wind is roughly 80% of the equatorial rotation rate. At the most extreme snapshot of the solar maxima, where the stellar-wind region is at its largest (corresponding to the highest opening colatitude), $\omega_{\odot} \approx 0.99$, i.e., the “effective” rotation of the wind is determined predominately by the most-rapid flows close to the equatorial value.

We determine the change in the predicted solar torque using our modified prescription, compared to the solid-body prediction. For the solid-body approximation, we adopt the Carrington sidereal rotation period of 25.4 days, which corresponds to the period it takes for fixed features on the solar surface at a colatitude of 64° (typical for sunspots and other solar activity) to rotate to the same apparent position when viewed from Earth (Carrington 1863, pp. 221, 244); this corresponds to a rotation rate of $\Omega_{\odot,SB} = 2.86 \times 10^{-6}$ rad s$^{-1}$, giving $f_{\odot,SB} = 4.56 \times 10^{-3}$. Using Equation (32) and assuming magnetocentrifugal effects on the solar-wind efficiency are negligible ($\omega_{\odot}f_{\odot,eq} \sim f_{\odot,SB} \ll K_{*,2}$), the ratio of our torque prediction and the solid-body prediction for the solar case can be interpreted as

$$ \frac{\dot{J}_{\odot}}{\dot{J}_{\odot,SB}} = \frac{\Omega_{\odot,eq}}{\Omega_{\odot,SB}} \omega_{\odot}. \qquad (33) $$

In the bottom panel of Figure 7, we plot the percentage change in the predicted solar torque as a function of time. At the solar minima, our formulation predicts a solar torque that is roughly 18% smaller than the solid-body prediction. Therefore, for the majority of the solar cycle, the open magnetic field extends over a region that is rotating slower than the solid-body rotation rate, resulting in a weaker torque. However, at the most extreme snapshot of the solar maxima, our formulation actually predicts a solar torque that is roughly 2.5% larger, but this is

primarily due to the fact that $\Omega_{\odot,eq} > \Omega_{\odot,SB}$. In this extreme case, where a large fraction of surface field is opened into the wind, the open region engulfs surface regions that are rotating faster than the solid-body value, and thus the torque is increased.

This analysis makes use of only the axisymmetric dipole field, which may be a reasonable approximation for most of the solar cycle except for at solar maximum, where global magnetic field configurations are often nonaxisymmetric and multipolar (DeRosa et al. 2012). However, this work could be a foundation for a more detailed analysis in which the average rotation rate of the open field regions is calculated during the solar cycle, i.e., accounting for nonaxisymmetry and multipolar magnetic fields. Some observations and theoretical work suggest that while the solar photosphere rotates as described, the solar corona (and thus the base of the solar wind) may in fact be rotating differently, or even potentially as a solid body (see, e.g., Insley et al. 1995; Giordano & Mancuso 2008; Bagashvili et al. 2017; Pinto et al. 2021). It is likely that the rotation of the solar corona determines the rotation rate of the wind, therefore it is possible that these systems are much more complex than anticipated by these formulations and that surface differential rotation may not be noticeably responsible for any torque discrepancies between observations and theory.

### 5.3. Significance of Differential Rotation for Main-sequence Stellar Spin-down

Most F, G, and K MS stars are observed to have a relative differential rotation rate similar to or smaller than the Sun, i.e., $\alpha < 0.3$, and most rapidly rotating stars are observed to be rotating almost as a solid body, i.e., $\alpha \to 0$ (see, e.g., Collier Cameron 2002; Reiners & Schmitt 2002b; Barnes et al. 2005; Reiners 2006; Collier Cameron 2007; Balona & Abedigamba 2016). For the cases similar to the Sun, differential rotation could affect the stellar torque on the order of tens of percent according to our prescription (as discussed in Section 5.2). The effect of differential rotation on wind torques is expected to be even smaller for more rapidly rotating stars, as $\alpha$ is observed to decrease with rotation period. In the case of antisolar rotation ($\alpha < 0$), which we did not explicitly explore in our parameter study, our formulation shows $\omega$ to be inversely proportional to the fractional open flux (i.e., proportional to $\Upsilon_{*}/\beta_{\Phi}$), suggesting that MS antisolar rotators actually spin down more efficiently. Finally, depending on the method used for determining rotation rates, it is likely that the observed rates are not generally equal to the maximum/equatorial rate, but rather to some intermediate rate. The observed rotation rates can thus be closer to the “effective” rotation of the wind, which would reduce the error in using solid-body stellar torque predictions.

Generally, depending on the precision required, models using solid-body rotation (and indeed using observations of rotation rates where the colatitude of the rotation rate observed is not known) are probably acceptable for studying stellar spin-down, especially because other uncertainties in the modeling (e.g., knowing the mass-loss rates or magnetic field properties) are likely comparable or larger. However, for stars demonstrating relative differential rotation comparable to, or stronger than, that of the Sun, the results presented in this paper should be considered.

L.G.I. and S.P.M. acknowledge support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement
No. 682393; sAWESoMeStars: Accretion, Winds, and Evolution of Spins and Magnetism of Stars; http://empslocal.ex.ac.uk/AWESoMeStars). A.J.F. is supported by the ERC Synergy grant "Whole Sun," No. 810218. C.Z. acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 742095; SPIDI: Star-Planets-Inner Disk-Interactions; http://spidi-eu.org).

The authors would like to acknowledge the use of the University of Exeter High-Performance Computing (HPC) facility in carrying out this work.

We thank Andrea Mignone and others for the development and maintenance of the PLUTO code. Figures within this work are produced using the Python package Matplotlib (Hunter 2007).

*Software:* Matplotlib (Hunter 2007), PLUTO (Mignone et al. 2007, 2012).

## ORCID iDs

Lewis G. Ireland https://orcid.org/0000-0002-8833-1204

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

Charlie R. Davey https://orcid.org/0000-0001-6734-0104

Owain L. Harris https://orcid.org/0000-0002-4198-2499

Tobias W. Slade-Harajda https://orcid.org/0000-0002-4993-0757

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

Claudio Zanni https://orcid.org/0000-0003-0204-8190

## References

Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77

Bagashvili, S. R., Shergelashvili, B. M., Japaridze, D. R., et al. 2017, A&A, 603, A134

Balona, L. A., & Abedigamba, O. P. 2016, MNRAS, 461, 497

Barnes, J. R., Cameron, A. C., Donati, J. F., et al. 2005, MNRAS, 357, L1

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

Beck, J. G. 2000, SoPh, 191, 47

Blackman, E. G., & Owen, J. E. 2016, MNRAS, 458, 1548

Bouvier, J. 2008, A&A, 489, L53

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

Bouvier, J., Matt, S. P., Mohanty, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Arizona Univ. Press)

Breimann, A. A., Matt, S. P., & Naylor, T. 2021, ApJ, 913, 75

Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192

Carrington, R. C. 1863, Observations of the Spots on the Sun from November 9, 1853, to March 24, 1861, Made at Redhill (London: Williams & Norgate)

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

Collier Cameron, A. 2002, AN, 323, 336

Collier Cameron, A. 2007, AN, 328, 1030

Collier Cameron, A., Donati, J. F., & Semel, M. 2002, MNRAS, 330, 699

Davenport, J. R. A. 2017, ApJ, 835, 16

Denissenkov, P. A., Pinsonneault, M., Terndrup, D. M., & Newsham, G. 2010, ApJ, 716, 1269

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

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

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

Finley, A. J., See, V., & Matt, S. P. 2019, ApJ, 876, 44

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, A&A, 595, A110

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

Giordano, S., & Mancuso, S. 2008, ApJ, 688, 656

Gondoin, P. 2017, A&A, 599, A122

Hunter, J. D. 2007, CSE, 9, 90

Insley, J. E., Moore, V., & Harrison, R. A. 1995, SoPh, 160, 1

Irwin, J., & Bouvier, J. 2009, in IAU Symp. 258, The Ages of Stars, ed. R. M. Ros et al. (Cambridge: Cambridge Univ. Press), 363

Jardine, M., Vidotto, A. A., & See, V. 2016, MNRAS, 465, L25

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

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

Kitchatinov, L. L., & Rüdiger, G. 1999, A&A, 344, 911

Küker, M., Rüdiger, G., & Kitchatinov, L. L. 2011, A&A, 530, A48

Lanzafame, A. C., & Spada, F. 2015, A&A, 584, A30

Lionello, R., Linker, J. A., Mikić, Z., & Riley, P. 2006, ApJL, 642, L69

Lionello, R., Titov, V. S., Mikić, Z., & Linker, J. A. 2020, ApJ, 891, 14

MacGregor, K. B., & Brenner, M. 1991, ApJ, 376, 204

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

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

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

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

Meibom, S., Mathieu, R. D., Stassun, K. G., Liebesny, P., & Saar, S. H. 2011, ApJ, 733, 115

Mestel, L. 1968, MNRAS, 138, 359

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

Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7

Morgan, H. 2011, ApJ, 738, 189

Morin, J., Donati, J. F., Petit, P., et al. 2010, MNRAS, 407, 2269

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

Pantolmos, G., Zanni, C., & Bouvier, J. 2020, A&A, 643, A129

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

Parker, E. N. 1965, SSRv, 4, 666

Perri, B., Brun, A. S., Réville, V., & Strugarek, A. 2018, JPIPh, 84, 765840501

Pinto, R. F., Poirier, N., Rouillard, A. P., et al. 2021, A&A, 653, A92

Reiners, A. 2006, A&A, 446, 267

Reiners, A., & Mohanty, S. 2012, ApJ, 746, 43

Reiners, A., & Schmitt, J. H. M. M. 2002a, A&A, 384, 155

Reiners, A., & Schmitt, J. H. M. M. 2002b, A&A, 393, L77

Reinhold, T., Reiners, A., & Basri, G. 2013, A&A, 560, A4

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

Roe, P. L. 1981, JCoPh, 43, 357

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

Sauty, C., Meliani, Z., Lima, J. J. G., et al. 2011, A&A, 533, A46

Schatzman, E. 1962, AnAp, 25, 18

See, V., Jardine, M., Vidotto, A. A., et al. 2016, MNRAS, 466, 1542

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

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

See, V., Matt, S. P., Folsom, C. P., et al. 2019a, ApJ, 876, 118

Skumanich, A. 1972, ApJ, 171, 565

Snodgrass, H. B., & Ulrich, R. K. 1990, ApJ, 351, 309

Stauffer, J., Rebull, L., Bouvier, J., et al. 2016, AJ, 152, 115

Tóth, G., Sokolov, I. V., Gombosi, T. I., et al. 2005, JGRA, 110, A12226

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

van Saders, J. L., & Pinsonneault, M. H. 2013, ApJ, 776, 67

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

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

Wood, B. E., Müller, H.-R., Redfield, S., et al. 2021, ApJ, 915, 37

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