# Numerical quantification of the wind properties of cool main sequence stars

Judy J. Chebly 1,2, Julián D. Alvarado-Gómez 1, Katja Poppenhäger1,2, Cecilia
Garraffo3,
1Leibniz Institute for Astrophysics, An der Sternwarte 16, 14482, Potsdam,
Germany
2Institute of Physics and Astronomy, University of Potsdam, Potsdam-Golm,
14476, Germany
3Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA
02138, USA E-mail<EMAIL_ADDRESS>0000-0003-0695-6487 0000-0001-5052-3473
0000-0003-1231-2194 0000-0002-8791-6286

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

As a cool star evolves, it loses mass and angular momentum due to magnetized
stellar winds which affect its rotational evolution. This change has
consequences that range from the alteration of its activity to influences over
the atmosphere of any orbiting planet. Despite their importance, observations
constraining the properties of stellar winds in cool stars are extremely
limited. Therefore, numerical simulations provide a valuable way to understand
the structure and properties of these winds. In this work, we simulate the
magnetized winds of 21 cool main-sequence stars (F-type to M-dwarfs), using a
state-of-the-art 3D MHD code driven by observed large-scale magnetic field
distributions. We perform a qualitative and quantitative characterization of
our solutions, analyzing the dependencies between the driving conditions
(e.g., spectral type, rotation, magnetic field strength) and the resulting
stellar wind parameters (e.g., Alfvén surface size, mass loss rate, angular
momentum loss rate, stellar wind speeds). We compare our models with the
current observational knowledge on stellar winds in cool stars and explore the
behaviour of the mass loss rate as a function of the Rossby number.
Furthermore, our 3D models encompass the entire classical Habitable Zones (HZ)
of all the stars in our sample. This allows us to provide the stellar wind
dynamic pressure at both edges of the HZ and analyze the variations of this
parameter across spectral type and orbital inclination. The results here
presented could serve to inform future studies of stellar wind-magnetosphere
interactions and stellar wind erosion of planetary atmospheres via ion escape
processes.

###### keywords:

exoplanets – stars: atmospheres – stars: magnetic fields – stars: mass-loss –
stars: winds, outflows

††pubyear: 2023††pagerange: Numerical quantification of the wind properties of
cool main sequence stars–A

##  1 Introduction

For many decades, scientists have known that the Sun has a mass outflow, which
is most visible in the behavior of comet tails (e.g., Biermann 1957). It has
also been established that solar wind is a natural byproduct of the heating
processes that produce the hot solar corona ($T~{}\sim~{}10^{6}$ K). As a
result, all cool main-sequence stars ($M_{\bigstar}~{}\leqslant$ 1.3 M⊙) with
analogous hot coronae, evidenced from their measured X-ray properties (Schmidt
et al. 1995; Pizzolato et al. 2000; Wright et al. 2011a), should have similar
winds (Parker 1958). Magnetic fields are thought to play a key role as an
energy source for the corona and the expanding solar atmosphere (e.g.,
Aschwanden 2006; Klimchuk 2015; Velli et al. 2015). Recent theories have shown
that in addition to magnetic fields, wave dissipation (via turbulence) and
magnetic reconnection could also play a role in energizing and shaping the
spatial properties of the solar wind (see, Ofman 2010; Cranmer 2012; Hansteen
& Velli 2012; Cranmer et al. 2015).

Winds, even if relatively weak, play an important role in stellar evolution
for stars of different spectral types causing the star to lose angular
momentum and slow its rotation over time (Weber & Davis 1967; Skumanich 1972;
Matt et al. 2012; Gallet & Bouvier 2013, 2015; Johnstone et al. 2015a, b;
Ahuir et al. 2020). As a result, the magnetic activities that constitute the
space weather (i.e., stellar winds, flares, coronal mass ejections) will
decrease with age in low-mass stars (Skumanich 1972; Güdel et al. 1997; Ribas
et al. 2005; Vidotto et al. 2014a). These changes in the host star will also
affect the evolution of planetary atmospheres and habitability (Tian 2008;
Kislyakova et al. 2014; Airapetian et al. 2017).

Direct measurements of the solar wind by spacecraft such as the Advanced
Composition Explorer (ACE, Stone et al. 1998; McComas et al. 1998), Ulysses
(McComas et al., 2003), and Parker solar probe (Kasper et al., 2021) have
improved our knowledge and understanding of its properties. On the other hand,
detecting a solar-like wind emitted by another star has proven extremely
challenging. This is not surprising, given how difficult it is to observe the
solar wind remotely. The latter carries a very low mass loss rate
($\dot{M}_{\odot}=\rm{2\times 10^{-14}~{}\dot{M_{\odot}}~{}{\rm yr^{-1}}}$,
see Feldman 1977; Wood 2004), which implies relatively low densities (near the
heliopause: $\sim$ 0.002 cm-3, Gurnett & Kurth 2019). Similarly, its high
temperature and elevated ionization state, make it difficult to detect with
simple imaging or spectroscopic techniques. As a result, properties such as
the associated mass loss rates, angular momentum loss rates, and terminal
velocities, crucial to understand stellar winds in low-mass stars, remain
poorly constrained.

Attempts to directly detect thermal radio emission from the plasma stream in
cool stars have not yet led to any discovery (Lim & White 1996; Drake et al.
1993; van den Oord & Doyle 1997; Gaidos et al. 2000; Villadsen et al. 2014;
Fichtinger et al. 2017). Current radio telescopes are not optimized for this
method; they can only detect winds much stronger than those from the Sun.
Moreover, the coronae of these active stars are also radio sources, making it
difficult to determine the exact source of the emission. Nevertheless, this
method has been able to establish upper limits for solar analogs of 1.3
$\times$ $10^{\rm-10}$ $\dot{M}_{\rm\odot}$ yr-1 (Gaidos et al. 2000;
Fichtinger et al. 2017). Another proposed method for direct detection is to
look for X-ray emission from nearby stars. As the star’s winds propagate, they
collide with the Local Interstellar Medium (ISM), forming "astrospheres"
similar to the Sun’s heliosphere (Wood, 2004). The charge exchange between the
highly ionized stellar wind and the ISM produces X-ray photons with energies
ranging from 453 to 701 eV. However, this method was unable to detect
circumstellar charge exchange X-ray emission even from the nearest star,
Proxima Centauri (Wargelin & Drake, 2002).

Similar to the charge exchange X-ray emission method, the Ly-$\rm\alpha$
absorption technique assumes the presence of the charge exchange phenomenon.
In this case, however, we are interested in the neutral hydrogen wall formed
at the astrospherical outer boundary by the interaction between the stellar
wind and the ISM. This exchange has been detected as excess HI Ly-$\rm\alpha$
absorption in Hubble Space Telescope UV stellar spectra (Linsky & Wood, 2014).
With nearly 30 measurements to date, spectroscopic analyses of the stellar HI
lines have proven to be the best method to unambiguously detect and measure
weak solar-like winds as well as some evolved cool stars (Wood et al., 2021).

Using this method,Wood et al. (2005b) found evidence for some increase in
$\dot{M}$ with magnetic activity, corresponding to a power-law relation in the
form $\dot{M}\propto F_{\rm X}^{\rm 1.34\pm 0.18}$ with $F_{\rm X}<10^{\rm 6}$
$\rm erg~{}cm^{\rm-2}~{}s^{\rm-1}$. However, this relation does not seem to
hold anymore for more active stars ($F_{\rm X}>10^{\rm 6}$ $\rm erg$ cm-2
s-1), mainly M-dwarfs (Wood et al., 2005b, 2014). Recently, Wood et al. (2021)
established a power law ($\dot{M}\propto F_{\rm X}^{\rm 0.77\pm 0.04}$)
between the $\dot{M}$ per unit surface area and the X-ray surface flux for
coronal winds for a broader selection of stars, including G, K, and new
$\dot{M}$ estimates for M-dwarfs. They found that the relation breaks even for
stars with $F_{\rm X}<10^{\rm 6}$ $\rm erg$ cm-2 s-1 (e.g., GJ 436, which has
$F_{\rm X}$ = 4.9 $\times 10^{\rm 4}$ erg cm-2 s-1, where the $\dot{M}$ was
estimated by using the planet as a probe for the stellar wind Vidotto &
Bourrier 2017) with the magnetic topology being a possible factor for the
scatter.

While extremely useful, the search for astrospherical absorption is influenced
by a number of critical factors. For instance, this method is strongly
dependent on the relative velocity of the stellar rest frame and the ISM flow
velocity ($V_{\rm ISM}$). As well as on the angle, $\theta$, between the
upwind direction of the ISM flow and the line-of-sight to the star (Wood et
al., 2021). It also requires prior knowledge of the properties of the ISM such
as the density and its ionization state (Wood et al. 2005a; Redfield & Linsky
2008a). Finally, its applicability is limited to relatively nearby stars
($\lesssim 15$ pc) due to the absorption of the ISM.

Due to the scarcity of observational data and associated limitations,
numerical simulations can be used to improve our understanding of stellar
winds. Models based on Alfvén waves are more commonly used to simulate the
stellar wind from stars other than the Sun (Suzuki, 2006). This is because
these waves are considered to be key mechanism for heating and accelerating
the solar wind (van der Holst et al. 2014; Van Doorsselaere et al. 2020).

In this study, we present a detailed numerical characterization of the stellar
wind properties of cool main-sequence stars (early F to M-dwarfs) covering a
range of rotation rates and magnetic field strengths. We compute steady-state
stellar wind solutions using a state-of-the-art 3D MHD model and provide
consistent qualitative and quantitative comparisons. Our goal is to better
understand the different stellar wind properties as a function of the driving
parameters, allowing us to explore the expected stellar wind conditions in the
circumstellar region around planet-hosting stars.

This paper is organized as follows: Section 2 describes the numerical model
and properties of the selected stellar sample. In Sect. 3, we present our
numerical results, discuss the derived trends in the stellar wind properties,
and compare our results with observations. This information is then used to
quantify the stellar wind conditions and explore their implications in the
context of the classical habitable zone (HZ) around cool main-sequence stars.
Conclusions and summary are provided in Sect. 4.

Figure 1: Simulated stellar wind environment for HD 179949 including multi-
domain models. The bottom panel contains the surface field distribution (in
Gauss) provided by the ZDI map and used to drive the AWSoM solution within the
stellar corona (SC) domain (middle panel). The blue-green color bar represents
the radial magnetic field strength on the stellar surface. Within SC, the gray
iso-surface corresponds to the Alfvén surface of the stellar wind (see Sect.
3.1). Selected magnetic field lines are shown in white. The steady-state
solution is propagated from the coupling region (62 - 67R★) to the entire
Inner Astrosphere (IA) domain (1200 R★ in each cartesian direction; upper
panel). The central gray sphere in the top panel denotes the boundary of IA
with the SC domain at 67 $R_{\bigstar}$. This domain contains the inner and
outer edges of the habitable zone (gray circles). Color-coded (top and middle
panel) is the wind dynamic pressure ($P_{\rm dyn}$ = $\rho U^{\rm 2}$),
normalized to the nominal Sun-Earth value ($\simeq$ 1.5 nPa), visualized on
the equatorial plane of both domains. The z-axis indicates the assumed stellar
rotation axis of the star. Figure 2: Examples of surface field distribution
(in Gauss) of our sample stars retrieved from ZDI maps. The rows represent
different spectral types going from late F-(top) to M-dwarfs (bottom) as
indicated. The color-code represents the normalized radial magnetic field for
a given row. The z-axis indicates the assumed stellar rotation axis for all
the stars in the panel. The slowest rotation in our sample of 21 stars is HD
219134 (K3V, $P_{\rm rot}$ = 42.2 d), and the fastest is GJ 1245 B (M6V,
$P_{\rm rot}$ = 0.71 d). The radial magnetic field strength ranges from 5 G to
1.5 kG.  Table 1: Fundamental parameters of our sample. Columns 1–8 list the
star number, name, spectral type (SpT), stellar mass ($M_{\bigstar}$), stellar
radius ($R_{\bigstar}$), rotation period ($P_{\rm\scriptscriptstyle rot}$),
effective temperature ($T_{\rm\scriptscriptstyle eff}$), and luminosity,
respectively (See et al. 2019a and references therein).

ID number |  Star |  SpT |  $M_{\bigstar}$ [$M\odot$] |  $R_{\bigstar}$ [$R\odot$] |  $P_{\rm\scriptscriptstyle rot}$ [d] |  $T_{\rm eff}$ [K] |  $L_{\rm\bigstar}$ [$L_{\rm\odot}$]
---|---|---|---|---|---|---|---
$1$ |  $\tau$ Boo |  F7V |  1.34 |  1.46 |  3 |  6387 |  3.0
$2$ |  HD 179949 |  F8V |  1.21 |  1.19 |  7.6 |  6168 |  1.80
$3$ |  HD 35296 |  F8V |  1.06 |  1.1 |  3.48 |  6202 |  1.60
$4$ |  HN Peg |  G0V |  1.1 |  1.04 |  4.55 |  5974 |  1.20
$5$ |  HD 190771 |  G2V |  1.06 |  1.01 |  8.8 |  5834 $\pm$ 50 |  0.99
$6$ |  TYC 1987-509-1 |  G7V |  0.9 |  0.83 |  9.43 |  5550† |  0.52 $\pm$ 0.03
$7$ |  HD 73256 |  G8V |  1.05 |  0.89 |  14 |  5480† |  0.72
$8$ |  HD 130322 |  K0V |  0.79 |  0.83 |  26.1 |  5400 † |  0.5
$9$ |  HD 6569 |  K1V |  0.85 |  0.76 |  7.13 |  5170 |  0.36 $\pm$ 0.01
$10$ |  $\epsilon$ Eri |  K2V |  0.85 |  0.72 |  11 |  5125 $\pm$ 87 |  0.3 $\pm$ 0.06
$11$ |  HD 189733 |  K2V |  0.82 |  0.76 |  12.5 |  4939 |  0.34
$12$ |  HD 219134 |  K3V |  0.81 |  0.78 |  42.2 |  4835 † |  0.27
$13$ |  TYC 6878-0195-1 |  K4V |  0.65‡ |  0.64‡ |  5.72 |  4600† |  0.8 $\pm$ 0.32
$14$ |  61 Cyg A |  K5V |  0.66 |  0.62 |  34.2 |  4655 † |  0.15
$15$ |  HIP 12545 |  K6V |  0.58‡ |  0.57‡ |  4.83 |  4300 † |  0.4 $\pm$ 0.06
$16$ |  TYC 6349-0200-1 |  K7V |  0.54‡ |  0.54‡ |  3.39 |  4100† |  0.3 $\pm$ 0.02
$17$ |  DT Vir |  M0V |  0.59 |  0.53 |  2.85 |  3850† |  0.055
$18$ |  GJ 205 |  M1.5V |  0.63 |  0.55 |  33.6 |  3690† |  0.061$\pm$ 0.006
$19$ |  EV Lac |  M3.5V |  0.32 |  0.3 |  4.37 |  3267 |  0.013
$20$ |  YZ CMi |  M4.5V |  0.32 |  0.29 |  2.77 |  3129 |  0.012
$21$ |  GJ 1245 B |  M6V |  0.12 |  0.14 |  0.71 |  3030§ |  0.0016

‡ Reference from Vidotto et al. (2014b)
† Reference from Pecaut & Mamajek (2013)
§ Reference from Vidotto et al. (2013)

##  2 Model description

We simulate stellar winds in cool main-sequence stars using the state-of-the-
art Space Weather Modeling Framework (SWMF; Sokolov et al. 2013; van der Holst
et al. 2014; Gombosi et al. 2018). The SWMF is a set of physics-based models
(from the solar corona to the outer edge of the heliosphere) that can be run
independently or in conjunction with each other (Tóth et al., 2012). This
model uses the numerical schemes of the Block Adaptive Tree Solar Roe-Type
Upwind Scheme (BATS-R-US; Powell et al. 1999) MHD solver. For a detailed
description of the model, see Gombosi et al. (2021). The multi-domain solution
starts with a calculation using the Solar/Stellar Corona (SC) module which
incorporates the Alfvén Wave Solar Model (AWSoM; van der Holst et al. 2014).
This module provides a description of the coronal structure and the stellar
wind acceleration region. The simulation is then coupled to a second module
known as the Inner Heliosphere/Astrosphere111This module is formally labeled
IH within the SWMF, but since we are working with low-mass main sequence
stars, we will refer to it as the Inner Astrosphere (IA) domain. (IA). In this
way, it is possible to propagate the stellar wind solution up to Earth’s orbit
and beyond. The model has been extensively validated and updated employing
remote sensing as well as in-situ solar data (e.g., Oran et al. 2017; Sachdeva
et al. 2019; van der Holst et al. 2019).

AWSoM is driven by photospheric magnetic field data, which is normally
available for the Sun in the form of synoptic magnetograms (Riley et al.,
2014). A potential field source surface method is used to calculate the
initial magnetic field (more details in the following section). This
information is used by AWSoM to account for heating and radiative cooling
effects, as well as the Poynting flux entering the corona, and empirical
turbulent dissipation length scales. With the interplay between the magnetic
field distribution, the extrapolation of the potential field, and the
thermodynamic properties, the model solves the non-ideal magnetohydrodynamic
(MHD) equations for the mass conservation, magnetic field induction, energy
(coronal heating), and momentum (acceleration of the stellar wind). These last
two aspects are controlled by Alfvén waves propagating along and against the
magnetic field lines (depending on the polarity of the field). In the momentum
equation, the heat and acceleration contributions are coupled by an additional
term for the total pressure and a source term in the energy equation. The
numerical implementation is described in detail in van der Holst et al.
(2014). Once these conditions are provided, the simulation evolves all
equations locally until a global steady-state solution is reached.

###  2.1 Simulation parameters and setup

In our work, we apply the SWMF/AWSoM model to main-sequence F, G, K, and
M-type stars by assuming that their stellar winds are driven by the same
process as the solar wind. We analyze the properties of the stellar wind by a
coupled simulation covering the region of the stellar corona (SC, spherical)
and the resulting structure within the inner astrosphere (IA, cartesian).
Figure 1 illustrates the coupling procedure in one of our models. This
coupling was necessary only in the case of F, G, and K stars, in order to
completely cover the habitable zones (HZ)222The range of orbits around a star
in which an Earth-like planet can sustain liquid water on its surface., which
are larger and farther away from the star. Parameters such as stellar radius
($R_{\bigstar}$), mass ($M_{\bigstar}$), and rotation period ($P_{\rm rot}$),
are also taken into account in the simulations. We followed the approach in
Kopparapu et al. (2014) in order to determine the optimistic HZs boundaries of
each star in our sample.

####  2.1.1 Simulation domain

The star is positioned in the center of the SC spherical domain. The radial
coordinate in SC ranges from 1.05 $R_{\bigstar}$ to 67 $R_{\bigstar}$, except
for M-dwarfs, where it extends to 250 $R_{\bigstar}$. The choice of the outer
edge value of the SC domain was chosen in a way to obtain both edges of the HZ
in one domain. The habitable zones limits were calculated using Kopparapu et
al. (2014) approach and the reported measured $L_{\bigstar}$ and $T_{\rm eff}$
for each star in our sample (see Table 1). As will be discussed in Sect. 3, in
the case of M-dwarfs, the extension had to be performed in order to cover the
entire Alfvén surface (AS)333This structure sets the boundary between the
escaping wind and the magnetically coupled outflows that do not carry angular
momentum away from the star., while keeping the default parameters for AWSoM
fixed (see Sect. 2.1.3). The domain uses a radially stretched grid with the
cartesian z-axis aligned with the rotation axis. The cell sizes in the
meridional ($\phi$) and azimuthal ($\theta$) directions are fixed at $\sim$
$2.8^{\circ}$. The total number of cells in the SC domain is $\sim$ $8\times
10^{5}$.

The steady-state solutions obtained within the SC module are then used as
inner boundary conditions for the IA component. An overlap of 5 $R_{\bigstar}$
(from 62 $R_{\bigstar}$ to 67 $R_{\bigstar}$) is used in the coupling
procedure between the two domains for F, G, and K stars (more details on the
necessity of the overlap when coupling between domains can be found in Tóth et
al. 2005). The IA is a cube that extends from $62~{}R_{\bigstar}$ to
$600~{}R_{\bigstar}$ in each cartesian component. Adaptive Mesh refinement
(AMR) is performed within IA, with the smallest grid cell size of $\sim$ 1.17
$R_{\bigstar}$ increasing up to 9.37 $R_{\bigstar}$ with a total of 3.9
million cells. As the simulation evolves, the stellar wind solution is
advected from SC into the larger IA domain where the local conditions are
calculated in the ideal MHD regime.

####  2.1.2 Magnetic boundary conditions

In the initial condition of the simulation, observations are used to set the
radial component of the magnetic field $B_{\rm r}$ [G] anchored at the base of
the wind (at the inner boundary). As mentioned earlier, a finite potential
field extrapolation procedure is carried out to obtain the initial
configuration of the magnetic field throughout SC (Tóth et al., 2011). This
procedure requires setting an outer boundary (source surface, $r_{\rm s}$),
beyond which the magnetic field can be considered to be purely radial and
force-free. The magnetic field can therefore be described as a gradient of a
scalar potential and determined by solving Laplace’s equation in the domain.
For the simulations discussed here, we set $r_{\rm s}$ at $45\%$ of the SC
domain size for F, G, and K stars, and $70\%$ for M-dwarfs. While the choice
of this parameter does not alter significantly the converged solutions, it can
modify the required run time of each model to achieve convergence. Therefore,
our selection was done to guarantee convergence to the steady-state in a
comparable number of iterations between all spectral types.

The stellar magnetic field as reconstructed from Zeeman Doppler Imaging
(ZDI)444A tomographic imaging technique that allows the reconstruction of the
large-scale magnetic field (strength and polarity at the star’s surface from a
series of polarized spectra (see e.g., Donati et al. 2006; Morin et al. 2008;
Fares et al. 2009; Alvarado-Gómez et al. 2015; Hussain et al. 2016; Kochukhov
2020)., is used as the inner boundary condition of SC (Fig. 2). Therefore, the
resulting wind solutions are more realistic than models based on
simplified/idealized field geometries (Chebly et al., 2021). Although the
reconstructed maps provide the distribution of vector magnetic fields, we use
only the radial component of the observed surface field. The magnetogram is
then converted into a series of spherical harmonic coefficients with a
resolution similar to that of the original map. The order of the spherical
harmonics should be chosen so that artifacts such as the "ringing" effect do
not appear in the solution (Tóth et al., 2011). In our models, we performed
the spherical harmonics expansion up to $l_{\rm max}$ = 5.

####  2.1.3 Input parameters

After we set the initial conditions, we define several parameters for the
inner boundary. In order to reduce the degree of freedom of the parameter set,
we only modify the parameters related to the properties of the stars, such as
mass, rotation period, and radius. As for the other parameters, we implement
the same values that are commonly used in the solar case (van der Holst et al.
2014; Sachdeva et al. 2019). The Poynting flux ($S/B_{\rm\bigstar}=1.1\times
10^{6}$ J m-2 s-1 T) is a parameter that determines the amount of wave energy
provided at the base of a coronal magnetic field line. The other parameter is
the proportionality constant that controls the dissipation of Alfvén wave
energy into the coronal plasma and is also known as the correlation length of
Alfvén waves ($L\rm{{}_{\rm\perp}=1.5\times 10^{5}}~{}m~{}\sqrt{T}$). We use
the values given in Sokolov et al. (2013) to define the base temperature
($T_{\rm o}=2\times 10^{6}K$) and the base density ($n\rm{{}_{\rm o}=2\times
10^{11}}cm^{-3}$).

We note that the choice of these parameters will affect the simulation
results, as reported in several studies that followed different approaches
(e.g., Boro Saikia et al. 2020a; Jivani et al. 2023). Recently, Jivani et al.
(2023) performed a global sensitivity analysis to quantify the contributions
of model parameter uncertainty to the variance of solar wind speed and density
at 1 au. They found that the most important parameters were the phostospheric
magnetic field strength, $S/B_{\rm\bigstar}$, and $L\rm{{}_{\rm\perp}}$.
Furthermore, in Boro Saikia et al. (2020a), an increase in the mass loss rate
($\dot{M}_{\bigstar}$), and angular momentum loss rate ($\dot{J}_{\bigstar}$)
was reported when $S/B_{\rm\bigstar}$ is increased from the solar value to
$2.0\times 10^{6}$ J m-2 s-1 T), which is expected because $S/B_{\rm\bigstar}$
drives the energy of the Alfvén wave, resulting in higher $\dot{M}_{\bigstar}$
and $\dot{J}_{\bigstar}$.

In this work, however, we are interested in isolating the expected
dependencies with the relevant stellar properties (e.g., mass, radius,
rotation period, photospheric magnetic field) which can only be analyzed
consistently if the AWSoM related parameters are kept fixed between spectral
types. Moreover, as will be discussed in detail in Sect. 3.3, the results
obtained using the standard AWSoM settings are either consistent with current
stellar wind observational constraints for different types of stars or the
apparent differences can be understood in terms of other physical factors or
assumptions made in the observations. For these reasons, we have chosen not to
alter these parameters in this study, which also reduces the degrees of
freedom in our models.

###  2.2 The sample of stars

Our investigation is focused on main sequence stars, with effective
temperatures ranging from 6500 K down to 3030 K, and masses
$M_{\rm\bigstar}<1.34~{}M_{\rm\odot}$ (spectral types F to M). All of these
stars are either fully or partially convective. We use a sample of 21 stars
whose large-scale photospheric magnetic fields were reconstructed with ZDI
(See et al. 2019b and references therein). Some of these stars were observed
at different epochs. In this case, the ZDI map with the best phase coverage,
signal-to-noise ratio, and most spectra used in the reconstruction was chosen.
The sample includes radial magnetic field strengths in the ZDI reconstruction
between $5$ G and $1.5$ kG corresponding to HD 130322 (K0V) and EV Lac
(M3.5V), respectively. Spectral types range from F7 ($\tau$ Boo,
$M_{\bigstar}=1.34~{}M_{\odot}$, $R_{\bigstar}=1.46~{}R_{\odot}$) to M6 (GJ
1245 B, $M_{\bigstar}=0.12~{}M_{\odot}$, $R_{\bigstar}=0.14~{}R_{\odot}$). The
rotation periods vary between fractions of a day to tens of days, with GJ 1245
B (M6V) having the shortest rotation period ($P_{\scriptscriptstyle\rm rot}$ =
0.71 d) and HD 219134 (K3V) the longest one ($P_{\scriptscriptstyle\rm rot}$ =
42.2 d). Table 1 contains the complete list of the sample stars and a summary
of the stellar properties incorporated in our models.

Figure 3: Simulation results in the SC domain for 3 K stars: HD 12545 (panel
A), HD 6569 (panel B), 61 Cyg A (panel C) driven by ZDI magnetic field maps.
All panels contain the projection onto the equatorial plane ($z$ = 0) of the
radial wind velocity ($U_{\rm r}$).The translucent gray shade denotes the
Alfvén surface calculated from the steady-state solution. The corresponding
color scale $U_{\rm r}$ is preserved among the different panels. Selected 3D
magnetic field lines are shown in black. The absolute size of the SC domain is
indicated in each case.  Figure 4: Simulated average Alfvén surface
($\overline{AS}$, top), mass loss rate per unit surface area
($\dot{M}_{\bigstar}/R^{\rm 2}_{\rm\bigstar}$, middle), and angular momentum
loss rate ($\dot{J}$, bottom) as a function of the spectral type (left) and
the surface-average radial magnetic field ($B_{\rm R}^{\rm avg}$, right). The
mass loss rate is expressed in units of the average solar values
($\dot{M}_{\rm\odot}=2\times 10^{\rm-14}~{}M_{\rm\odot}~{}{\rm
yr^{\rm-1}}=1.265\times 10^{\rm 12}~{}{\rm g~{}s^{-1}}$), normalized to the
surface area of each star ($A_{\bigstar}=4\pi R_{\bigstar}^{2}$). Individual
points denote the results of each simulation presented in Sect. 3, Table 2.
The different symbols and colors represent the spectral types (F,
cyan/squares; G, yellow/diamonds; K, orange/circles; M, red/star). The purple
line and shaded area represent the fitted power-law with its uncertainties.

##  3 Results & Discussion

###  3.1 The effect of star properties on the wind structure

The Alfvén surface (AS) is defined by the collection of points in the 3D space
that fulfils the Alfvén radius criterion555The Alfvén radius ($R_{\rm A}$) is
defined as the distance around a star at which the kinetic energy density of
the stellar wind equals the energy density of the astrospheric magnetic
field.. Numerically, it is determined by finding the surface for which the
wind velocity reaches the local Alfvén velocity, $v_{A}=B/\sqrt{4{\pi\rho}}$,
where $B$ and $\rho$ are the local magnetic field and plasma density,
respectively. The Alfvén surface can be interpreted as the lever arm of the
wind torque –the "position" at which the torque acts to change the angular
rotation of the star666In other words, the angular momentum per unit mass
within the stellar wind can be computed as if there were solid body rotation,
at an angular velocity $\Omega_{\bigstar}=2\pi/P_{\rm rot}$, out as far as the
Alfvén surface.. The Alfvén Surface is used in numerical models to
characterize ($\dot{M}_{\bigstar}$) and ($\dot{J}_{\bigstar}$) (e.g., Vidotto
et al. 2015; Boro Saikia et al. 2020a; Garraffo et al. 2015b; Alvarado-Gómez
et al. 2016). We compute $\dot{J}_{\bigstar}$ by performing a scalar flow rate
integration over the AS and another one over a closed spherical surface (S)
beyond the AS to determine $\dot{M}_{\bigstar}$:

$\dot{M}_{\bigstar}=\int_{S}\rho(\textbf{u}\cdot\rm\textbf{dA})$ (1)
$\dot{J}_{\bigstar}=\int_{AS}\rm\Omega\rho\rm
R^{2}\sin^{2}\theta(\textbf{u}\cdot\rm\textbf{dA})$ (2)

Here $\dot{J}_{\bigstar}$ is the component of the change in angular momentum
in the direction of the axis of rotation. The distance to the Alfvén surface
is represented by $R$. The angle between the lever arm and the rotation axis
is denoted by $\theta$, which depends on the shape/orientation of the AS with
respect to the rotation axis (and accounted for in the surface integral). The
stellar angular velocity is represented by
$\Omega=2\rm\pi/P_{\scriptscriptstyle\rm rot}$. The surface element is denoted
by dA.

Figure 3 shows the AS of the stellar wind, with plasma streamers along with
the equatorial section flooded with the wind velocity ($U_{\rm r}$) for three
K stars in our sample (HIP 12545, HD 6569, 61 Cyg A). If we compare two stars
with similar $P_{\rm\scriptscriptstyle rot}$ but different $B_{\rm R}^{\rm
max}$, we can clearly see that the size of AS increases with increasing
magnetic field strength. This is a direct consequence of the dependence of the
Alfvén velocity on these quantities (Eq. 1) and the distance from the star at
which the Alfvén velocity is exceeded by the wind. For instance, for very
active stars with stronger magnetic fields, the expected coronal Alfvén
velocity is greater than for less active stars, increasing the radial distance
that the wind velocity must travel to reach the Alfvén velocity. The
associated Alfvén surface has a characteristic two-lobe configuration (Fig. 3,
gray translucent area), with average sizes of 27 $R_{\bigstar}$, 18
$R_{\rm\bigstar}$ and 13 $\rm R_{\bigstar}$ for HIP 12545, HD 6569, and 61 Cyg
A, respectively (see Table 2).

When we compare two stars with similar magnetic field strengths but different
$P_{\rm\scriptscriptstyle rot}$ (see Fig. 3, panels B and C), the change in AS
size is not as dramatic. The rotation period has primarily a geometric effect
on the resulting AS. The Alfvén surface assumes a different tilt angle in all
three cases. This tilt is mainly connected to the open magnetic field flux
distribution on the star’s surface (Garraffo et al., 2015a). We also notice in
Fig. 3 that the stellar wind distribution is mainly bipolar with a relatively
fast component reaching up to $\sim$ $891$ km s-1 for HIP 12545 $\sim$, $702$
km s-1 for HD 6569, and $\sim$ $593$ km s-1 for 61 Cyg A. In section 3.3.1 we
will discuss further the relation between the wind velocity with regard to
$P_{\rm rot}$ and $B_{\rm R}$.

Figure 4 shows the $\dot{M}_{\bigstar}$, $\dot{J}_{\bigstar}$, $\overline{AS}$
as estimated by the previously described method, against the sub-spectral type
of our star sample (left column) and the average radial magnetic field
strength ($B_{\rm R}^{\rm avg}$, right column). Similar relations have been
obtained for the maximum radial magnetic field strength and are presented in
Appendix A. The average Alfvén surface size was calculated by performing a
mean integral over the radius at each point of the 3D AS. The extracted
quantities are represented by different colors and symbols for each spectral
type (F, G, K, and M).

As expected, the AS increases as we move toward more magnetically active stars
(Fig. 4, top-right panel). From our simulations, we were able to establish a
relation between AS and $B_{\rm R}^{\rm avg}$ using the bootstrap technique
(1000 realizations) to find the mean of the slope and the intercept along with
their uncertainties. We use this approach to determine all relations from our
simulations. The relation is as follows:

$\log\overline{AS_{\rm R}}=(0.42\pm 0.06)\log B_{\rm R}^{\rm avg}+(0.71\pm
0.07)$ (3)

Our simulated steady-state $\dot{M}_{\bigstar}$ show a scatter within the
range [0.5 $\dot{M}_{\odot}$/$R_{\scriptscriptstyle\odot}^{2}$, 30
$\dot{M}_{\odot}$/$R_{\scriptscriptstyle\odot}^{2}$], which is comparable to
that estimated from the observed Ly$\rm\alpha$ absorption method of G, K, and
M-dwarfs in Wood et al. (2021). The variations in $\dot{M}_{\bigstar}$ are
related to differences in the strength and topology of the magnetic field
driving the simulations (see Alvarado-Gómez et al. 2016; Vidotto et al. 2016;
Evensberget et al. 2022), as well as to the Alfvén wave energy transfer to the
corona and wind implemented in the model (Boro Saikia et al. 2020b; Jivani et
al. 2023). For this reason, we tried to isolate the effects introduced by the
star (e.g., $M_{\bigstar}$, $R_{\bigstar}$, ${P_{\rm rot}}$, magnetic field
strength) over the ones from the Alfvén wave heating (i.e., $n_{\rm o}$,
$T_{\rm o}$, $S/B_{\bigstar}$, $L_{\perp}$).

In terms of mass loss rate, stronger winds are expected to be generated by
stronger magnetic fields (see Fig. 4) implying that the winds are either
faster or denser. This interplay determines $\dot{M}_{\bigstar}$ (Eq. 1),
which increases with increasing magnetic field strength regardless of spectral
type. We see a common increase for F, G, K, and M-dwarfs (excluding EV Lac) in
the saturated and unsaturated regime that can be defined from the simulations
as follows:

$\log\dot{M}_{\rm\bigstar}/R^{\rm 2}_{\bigstar}=(0.48\pm 0.09)\log B_{\rm
R}^{\rm avg}+(0.11\pm 0.10)$ (4)

On the other hand, we observe a slightly different behavior for M-dwarfs,
whose $\dot{M}_{\bigstar}$ and $\dot{J}_{\bigstar}$ values tend to be lower.
As discussed by Garraffo et al. (2018), the magnetic field complexity could
also affect $\dot{M}_{\bigstar}$ for a given field strength. We consider this
possibility in the following section. Note that, as has been shown in previous
stellar wind studies of M-dwarfs (e.g., Garraffo et al. 2017; Kavanagh et al.
2021; Alvarado-Gómez et al. 2022b), modifications to the base AWSoM parameters
(either in terms of the Poynting flux or the Alfvén wave correlation length)
would lead to strong variations in $\dot{M}_{\bigstar}$. This would permit
placing the M-dwarfs along the general trend of the other spectral types in
particular, the $\dot{M}_{\bigstar}$ value obtained for the star with the
strongest $B_{\rm R}$ in our sample (EV Lac). While these modifications have
physical motivations behind them (i.e. increased chromospheric activity,
stronger surface magnetic fields), in most regards, they remain unconstrained
observationally. Furthermore, the values we obtain in our fiducial AWSoM
models are still within the range of observational estimates available for
this spectral type (see Sect. 3.3.1), with the added benefit of minimizing the
degrees of freedom and isolating the effects of the stellar parameters on the
results.

Similarly, we see a large scatter of $\dot{J}_{\bigstar}$ with respect to the
spectral type (Fig. 4, bottom left column), ranging from $\rm{10^{26}~{}{\rm
g~{}cm^{2}~{}s^{-2}}}$ to $\rm{10^{31}~{}{\rm g~{}cm^{2}~{}s^{-2}}}$. This
range is within the expected $\dot{J}_{\bigstar}$ values estimated for cool
stars with the lowest value corresponding to M-dwarfs (See et al. 2019c and
references therein). The maximum $\dot{J}_{\bigstar}$ values reached in our
simulations are comparable to $\dot{J}_{\odot}$ reached at solar minimum and
maximum ($7\times 10^{\rm 30}$ and $10\times 10^{\rm 30}$ $\rm
g.cm^{2}s^{-2}$, Finley et al. 2018; Boro Saikia et al. 2020a). We note that
this is the only parameter for which we have retained units in absolute values
(as is commonly done in solar/stellar wind studies; see Cohen et al. 2010;
Garraffo et al. 2015a; Finley et al. 2018). Using absolute units, we expect a
decrease in $\dot{J}_{\bigstar}$ as we move from F to M-dwarfs, since
$\dot{J}_{\bigstar}$ is a function of $R^{\rm 2}$ (Eq. 2). The scatter around
this trend is dominated by the relatively small $\dot{M}_{\bigstar}$ values,
the distribution of $\Omega_{\bigstar}$ in our sample (variations up to a
factor of 5), and the equatorial AS size where the maximum torque is applied
($\sin{\theta}$ in Eq. 2).

We also note that the sample is biased toward weaker magnetic field strengths.
To better estimate how the magnetic field affects the properties of the
stellar winds, we need a larger sample, not only in terms of stellar
properties but also with stellar wind constraints such as
$\dot{M}_{\bigstar}$. The latter is so far the only stellar wind observable
parameter for which comparisons can be made. For this reason, we will focus on
the behavior of the $\dot{M}_{\bigstar}$ as a function of different stellar
properties in the following sections of the analysis.

###  3.2 Stellar mass-loss rate and complexity

Figure 5: Simulated mass loss rate ($\dot{M}_{\rm\bigstar}/R^{\rm
2}_{\rm\bigstar}$) versus the expected complexity value ($n$, Garraffo et al.
2018) given the Rossby number ($R_{\rm o}$) of each star. The mass loss rate
is expressed in solar units, normalized by the unit surface area of each star.
Individual points denote the results of each simulation presented in Sect. 3
(Table 2). The different colors (cyan, orange, yellow, red) represent the
different spectral types respectively (F, G, K, and M). The filled circles and
the empty ones represent stars in the saturated and unsaturated regimes,
respectively. The symbol size denotes the maximum magnetic field strength
($B_{\scriptscriptstyle\rm R}^{\rm max}$) in each case. The numbers refer to
the ID of the star in our sample. The complexity value n = 1 corresponds to
the dipole, and n increases as we move to the right.

Coronal X-ray luminosity is a good indicator of the level of magnetic activity
of a star and the amount of material heated to $10^{6}$ K temperatures. The
dependence of magnetic activity on dynamo action (i.e., dynamo number
$D=R_{\rm o}^{-2}$, Charbonneau 2020) has led a number of authors to use the
Rossby number to characterize stellar activity, for a wide range of stellar
types (Wright et al., 2011a). The Rossby number is defined as $R_{\rm
o}=P_{\rm rot}$/$\tau_{\rm c}$, where $P_{\rm rot}$ is the stellar rotation
period and $\tau_{\rm c}$ is the convective turnover time (Noyes et al. 1984;
Jordan & Montesinos 1991; Wright et al. 2011b). We adopted the approach of
Wright et al. (2018a) to calculate $\tau_{\rm c}$. In this case, the latter is
only a function of the stellar mass ($M_{\rm\bigstar}$):

$\log\tau_{\rm
c}=2.33-1.50(M_{\bigstar}/M_{\odot})+0.31(M_{\bigstar}/M_{\odot})^{2}$ (5)

As it was mentioned in Sect. 1, the study of Wood et al. (2021) suggests that
coronal activity increases with $\dot{M}_{\bigstar}$. The overall increase in
$\dot{M}_{\bigstar}$ with X-ray flux $F_{\rm X}$ ($\dot{M}_{\bigstar}\propto
F_{\rm X}^{0.77\pm 0.04}$), is most likely due to their dependence on magnetic
field strength (see Sect. 3.1). However, they report a scatter of about two
orders of magnitude of $\dot{M}_{\bigstar}$ around the trend line. This
suggests that coronal activity and spectral type alone do not determine wind
properties. The geometry of the magnetic field may also play a role.

The correlation between $\dot{M}_{\bigstar}$ and magnetic complexity has
already been suggested by Garraffo et al. (2015b), which could in principle
contribute to the scatter in (Wood et al. 2021, Fig. 10). The large-scale
distribution of the magnetic field on the stellar surface is mainly determined
by the rotation period and the mass of the star, namely $R_{\rm o}$ (Garraffo
et al. 2018; Morin et al. 2010; Gastine et al. 2013). The Rossby number was
used to determine the complexity function in Garraffo et al. (2018), which was
able to reproduce the bimodal rotational morphology observed in young open
clusters (OCs). The complexity function of Garraffo et al. (2018) is defined
as

$n=\frac{a}{R_{\rm o}}+1+bR_{\rm o}\mbox{.}$ (6)

The constant 1 reflects a pure dipole. The coefficients a = 0.02 and b = 2 are
determined from observations of OCs. The first term is derived from the ZDI
map observation of stars with different spectral types and rotation periods.
The third term is motivated by Kepler’s observations of old stars (van Saders
et al. 2016; van Saders et al. 2019; van Saders et al. 2016; van Saders et al.
2019).

We emphasize that the complexity number ($n$), estimated from Eq. 6, differs
from the complexity derived from the ZDI maps themselves (e.g., Garraffo et
al. 2022). The complexity number from $R_{\rm o}$ is expected to be higher.
This is due to the fact that many of the small-scale details of the magnetic
field are not captured by ZDI.

We expect to lose even more information about the complexity of the field
given that the ZDI maps are not really available to the community (apart from
the published images). Image-to-data transformation techniques (which we
applied to extract the relevant magnetic field information from the published
maps) can lead to some losses of information, both spatially and in magnetic
field resolution. These vary depending on the grid and the projection used to
present the ZDI reconstructions (i.e., Mercator, flattened-polar, Mollweide).
Using the star’s raw ZDI map would prevent these issues and would aid with the
reproducibility of the simulation results.

Finally, note that the expected complexity is also independent of the
spherical harmonic expansion order used to parse the ZDI information to the
simulations. The obtained $R_{\rm o}$ and $n$ values for each star in our
sample are listed in Table 2.

Figure 5 shows the behaviour of coronal activity and $\dot{M}_{\bigstar}$ with
respect to the expected magnetic field complexity (n). The coronal activity is
denoted by full and empty symbols corresponding to saturated and unsaturated
stars, respectively. We consider stars with $R_{\rm o}\leq 0.1$ in the
saturated regime and stars with $R_{\rm o}>0.1$ in the unsaturated regime
based on X-ray observations (Pizzolato et al. 2003; Wright et al. 2011b,
2018b; Stelzer et al. 2016; Wright & Drake 2016). The colors correspond to the
different spectral types, whereas the numbers indicate the ID of each star in
our sample. The symbol size represents the maximum radial magnetic field
strength of each star extracted from the ZDI observations.

We anticipate seeing a trend in which the $\dot{M}_{\bigstar}$ decreases as
the magnetic field complexity increases (leading to an increment of closed
loops on the stellar corona), for stars in saturated and unsaturated regimes.
For instance, $\epsilon$ Eri ($\\#10$, $B_{\scriptscriptstyle\rm
R}^{\scriptscriptstyle\rm max}=25$ G, n = 2.21724) has an
$\dot{M}_{\rm\bigstar}$ = 4.53 $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$ lower
than HD 6569 ($\\#9$, $B_{\scriptscriptstyle\rm R}^{\scriptscriptstyle\rm
max}=~{}\rm 29$ G, n = 1.80346 ) with $\dot{M}_{\bigstar}$ = 6.70
$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$. This is also true for $\tau$ Boo and
HD 179949 where $\tau$ Boo ($\\#1$, $B_{\scriptscriptstyle\rm
R}^{\scriptscriptstyle\rm max}=14$ G, n = 1.84728, $\dot{M}_{\bigstar}$ = 2.30
$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$) has a higher $\dot{M}_{\bigstar}$
compared to HD 179949 ($\\#2$, $B_{\scriptscriptstyle\rm
R}^{\scriptscriptstyle\rm max}=12$ G, n = 2.65746, $\dot{M}_{\bigstar}$ = 1.90
$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$).

We also noticed that as we go to more active stars, like in the case of
M-dwarfs, the field strength starts to dominate over the complexity in terms
of contribution to the $\dot{M}_{\bigstar}$. For example, GJ 1245 B ($\\#21$,
$B_{\scriptscriptstyle\rm R}^{\scriptscriptstyle\rm max}=~{}404$ G , n =
5.02602 , $\dot{M}_{\bigstar}$ = 9.27 $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm
2}$) has an $\dot{M}_{\bigstar}$ higher than DT Vir even though the complexity
of the former is almost 5 times higher (DT Vir, $\\#17$,
$B_{\scriptscriptstyle\rm R}^{\scriptscriptstyle\rm max}=~{}327$ G, n =
1.41024, $\dot{M}_{\bigstar}$ = 3.81 $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm
2}$). However, in order to better understand the contribution of the
complexity in $\dot{M}_{\bigstar}$, we will need to run simulations for a
wider range of stars with sufficiently high resolution of the driving magnetic
field to capture directly the complexity of the field (and not estimate it
from a scaling relation as it was performed here).

Moreover, our results show that whenever we have a case in which the star
properties ($M_{\bigstar}$, $R_{\bigstar}$, and $P_{\rm rot}$), magnetic field
strength and complexity are comparable, we end up with similar
$\dot{M}_{\bigstar}$. This will be the case of TYC 6878-0195-1 ($\\#13$,
$B_{\scriptscriptstyle\rm R}^{\scriptscriptstyle\rm max}=~{}162$ G, n =
1.48069, $\dot{M}_{\bigstar}$= 17.42 $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$)
and HIP 12545 ($\\#~{}15$, $B_{\scriptscriptstyle\rm R}^{\scriptscriptstyle\rm
max}=~{}184$ G, n = 1.41505, $\dot{M}_{\bigstar}$ = 20.11
$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$).

Furthermore, two stars with similar coronal activity with respect to X-ray
flux, i.e., EV Lac and YZ CMi ($F_{\scriptscriptstyle\rm X}\approx 10^{7}$
$\rm ergs$ $\rm cm^{-2}\rm s^{-1}$), but with slightly different magnetic
field complexity, result in different wind properties: respectively
$\dot{M}_{\rm\bigstar}=\rm 0.62$ $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$, and
$\dot{M}_{\rm\bigstar}=\rm 20.57$ $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$. A
similar situation occurs when two stars have a comparable field complexity but
different coronal activity i.e., YZ CMi and GJ 205 ($\\#18$,
$\dot{M}_{\rm\bigstar}=\rm 2.32$ $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$,
$F_{\scriptscriptstyle\rm X}\approx 10^{5}$ $\rm ergs$ cm${}^{-2}\rm s^{-1}$).

The lowest $\dot{M}_{\bigstar}$ corresponds to the saturated M-dwarf EV Lac
($\\#$19), which has the strongest $B_{\scriptscriptstyle\rm R}$ (1517 G) and
one of the simplest complexities in our sample (n = 1.46331). The low
complexity of the field means that the wind is dominated by open field lines,
leading to very high wind velocities in the standard AWSoM model, but with a
very low density, which in turn leads to small $\dot{M}_{\bigstar}$ values. We
remind the reader that the base density of the stellar wind is fixed at the
stellar surface and is the same for all the stars in the sample (Sect. 2.1.3).

Figure 6: Numerical results of the stellar mass loss rate
($\dot{M}_{\bigstar}$/$R^{\rm 2}_{\bigstar}$, circles), astrosphere stellar
mass loss rate (squares, Wood et al. 2021), slingshot prominences mass loss
rate (diamonds, Jardine & Collier Cameron 2019), and absorption during an
exoplanetary transit (plus within a square, Vidotto & Bourrier 2017) against
the Rossby number ($R_{\rm o}$). Colors illustrate the different spectral
types: cyan (F), yellow (G), orange (K), and red (M). The Sun is represented
by a star symbol. Dashed lines connect the common stars between our sample and
the ones with estimated $\dot{M}_{\bigstar}/R^{\rm 2}_{\bigstar}$ values by
the astrospheric Ly-$\alpha$ absorption technique. Black arrows pointing
downward correspond to the upper limits given by observations. Figure 7:
Simulated stellar wind environment for HD 179949, HD 73256, HD 189733, and DT
Vir. Multi-domain models for HD 179949, HD 73256, HD 189733 are shown. The
steady-state solution of the multi-domain cases is propagated from the
coupling region (62 $-$ 67 $R_{\bigstar}$) to the entire IA domain (1200 $\rm
R_{\bigstar}$ in each cartesian direction). The steady-state solution of DT
Vir is shown in the spherical domain (SC). The green iso-surface represents
the averaged wind velocity at 1 au for F, G, and K, as for M-dwarfs it
represents the highest averaged velocity. Color-coded is the wind dynamic
pressure ($P_{\rm dyn}=\rm{\rho U^{2}}$ ) normalized to the nominal Sun-Earth
value ($\sim$ 1.5 nPa), visualized on the equatorial plane of both domains and
on a translucent sphere (at R = 0.5 au). Selected magnetic field lines are
shown in black. The 2 white circles represent the optimistic habitable zone
boundaries. The white translucent sphere represents the coupling region (67
$R_{\bigstar}$) between the SC and IH domains in the case of F, G, and K
stars.

###  3.3 Stellar wind mass-loss rate and Rossby number

Table 2: A summary of the resulting stellar wind properties with their
corresponding driving parameters. Columns 1–11, list the star number, name,
mass-loss rate per unit surface area†
($\dot{{M}}_{\bigstar}/\rm{R}^{2}_{\bigstar}$), angular momentum loss rate
($\dot{{J}}_{\bigstar}$), average terminal velocity
($\overline{U_{\rm\scriptscriptstyle R}^{\rm\scriptscriptstyle 1au}}$),
average terminal velocity ($\overline{U_{\rm\scriptscriptstyle
R}^{\rm\scriptscriptstyle T}}$), average Alfvén surface radius
($\overline{AS}$), absolute maximum surface radial magnetic field ($|B_{\rm
R}|^{\rm max}$), average surface radial magnetic field ($B_{\rm R}^{\rm
avg}$), Rossby number‡ ($R_{\rm o}$), and the expected complexity number§
($n$).

ID number |  Name |  $\dot{M}_{\bigstar}/R^{2}_{\bigstar}$ |  $\dot{J}_{\bigstar}$ [$g.cm^{2}.s^{-1}$] |  $\overline{U_{\rm\scriptscriptstyle R}^{\rm\scriptscriptstyle 1au}}$ [km $s^{-1}$] |  $\overline{U_{\rm\scriptscriptstyle R}^{\rm\scriptscriptstyle T}}$ [km $s^{-1}$] |  $\overline{AS}$ [R★] |  $|B_{\rm R}|^{\rm max}$ [G] |  $B_{\rm R}^{\rm avg}$ [G] |  $R_{\rm o}$ |  n
---|---|---|---|---|---|---|---|---|---|---
$1$ |  $\rm\tau$ Boo |  2.30 |  1.35E+30 |  320 |  332 |  13 |  14 |  4.01 |  0.39855 |  1.84728
$2$ |  HD 179949 |  1.90 |  2.07E+29 |  345 |  356 |  11 |  12 |  3.91 |  0.81648 |  2.65746
$3$ |  HD 35296 |  7.33 |  1.20E+30 |  270 |  289 |  20 |  27 |  13.94 |  0.28396 |  1.63835
$4$ |  HN Peg |  8.42 |  1.07E+30 |  545 |  549 |  16 |  50 |  20.48 |  0.40079 |  1.85148
$5$ |  HD 190771 |  5.45 |  2.58E+29 |  431 |  432 |  14 |  24 |  11.75 |  0.71806 |  2.46397
$6$ |  TYC 1987-509-1 |  9.82 |  2.11E+29 |  547 |  530 |  17 |  54 |  25.63 |  0.55388 |  2.14387
$7$ |  HD 73256 |  3.04 |  5.18E+28 |  457 |  458 |  11 |  22 |  7.91 |  1.12036 |  3.25857
$8$ |  HD 130322 |  1.17 |  8.01E+27 |  436 |  440 |  9 |  5 |  2.08 |  1.19721 |  3.41113
$9$ |  HD 6569 |  6.70 |  7.88E+28 |  704 |  702 |  18 |  29 |  18.41 |  0.37507 |  1.80346
$10$ |  $\epsilon$ Eri |  4.53 |  6.89E+28 |  554 |  553 |  14 |  25 |  12.26 |  0.59172 |  2.21724
$11$ |  HD 189733 |  11.67 |  1.34E+29 |  535 |  534 |  16 |  51 |  28.43 |  0.61443 |  2.26141
$12$ |  HD 219134 |  1.50 |  4.97E+27 |  425 |  429 |  9 |  6 |  3.44 |  2.02732 |  5.06451
$13$ |  TYC 6878-0195-1 |  17.42 |  2.08E+29 |  734 |  692 |  24 |  162 |  60.14 |  0.18682 |  1.48069
$14$ |  61 Cyg A |  3.98 |  7.98E+27 |  609 |  593 |  13 |  18 |  10.07 |  1.14548 |  3.30842
$15$ |  HIP 12545 |  20.11 |  2.13E+29 |  906 |  891 |  28 |  184 |  106.62 |  0.13145 |  1.41505
$16$ |  TYC 6349-0200-1 |  16.55 |  1.92E+29 |  657 |  642 |  23 |  93 |  51.36 |  0.08314 |  1.40684
$17$ |  DT Vir |  3.81 |  2.67E+28 |  —- |  1102 |  37 |  327 |  125.08 |  0.07979 |  1.41024
$18$ |  GJ 205 |  2.32 |  2.28E+27 |  —- |  690 |  17 |  22 |  14.41 |  1.04304 |  3.10525
$19$ |  EV Lac |  0.62 |  4.10E+26 |  —- |  3675 |  122 |  1517 |  620.21 |  0.05738 |  1.46331
$20$ |  YZ CMi |  20.57 |  1.81E+28 |  —- |  1709 |  132 |  822 |  655.66 |  0.03637 |  1.62264
$21$ |  GJ 1245 B |  9.27 |  1.90E+27 |  —- |  1164 |  44 |  404 |  200.43 |  0.00498 |  5.02602

† Normalized to solar units ($\dot{M}_{\odot}/{R}^{2}_{\odot}=1.0$).
‡ Predicted by the empirical model of Wright et al. (2011a).

Using the results of our stellar winds models, we can study how the
$\dot{M}_{\bigstar}$ changes as a function of the Rossby number ($R_{\rm o}$).
The Rossby number is a useful quantity because it not only removes the
dependence on spectral type, but also relates the rotation period to magnetic
field strength, complexity, and even stellar coronal activity. The latter is
also important because cool stars exhibit a well-defined behavior between
$L_{\rm X}$ (or $F_{\rm X}$) and $R_{\rm o}$ (saturated and unsaturated
regimes). Thus, if we analyze $\dot{M}_{\bigstar}$ using this parameter, we
can see (to some extent) all dependencies simultaneously.

Figure 6 shows the stellar mass-loss rate per unit surface area
($\dot{M}_{\bigstar}/R^{\rm 2}_{\bigstar}$) as a function of the Rossby number
($R_{\rm o}$). The circles show our 3D MHD numerical results, while the empty,
filled, and the plus sign within a square corresponds to observational
estimates of astrospheres (Wood et al., 2021), slingshot prominences (Jardine
& Collier Cameron, 2019), and absorption during an exoplanetary transit
(Vidotto & Bourrier, 2017), respectively. We use the same method as for the
simulated stars (Eq. 5) to calculate the $R_{\rm o}$ of stars with constraints
on their mass loss rate. Spectral types are indicated by different colors:
cyan (F), yellow (G), orange (K), and red (M). The Sun is represented by a
yellow star symbol. Dashed lines connect the common stars in our models and
the observations. In this section, we will focus only on the resulting
$\dot{M}_{\bigstar}$ from the numerical results.

As was mentioned earlier, our 3D MHD simulated $\dot{M}_{\bigstar}$ values are
in the same range as the $\dot{M}_{\bigstar}$ estimates from the Ly-$\alpha$
astrospheric absorption method. Note that since we are only simulating steady-
state stellar winds, our comparison is mostly focused on the steady mass loss
$\dot{M}_{\bigstar}$ (filled squares and squares with a plus sign). As such,
it is not surprising that our $\dot{M}_{\bigstar}$ values appear 1 - 2 orders
of magnitude below the estimates associated with sporadic mass loss events
such as slingshot prominences in very active stars in the saturated regime
(filled squares, Jardine & Collier Cameron 2019).

Based on the relation between $F_{\rm X}$ and $R_{\rm o}$ (Wright et al.,
2018a), and the broad correlation observed between $\dot{M}_{\bigstar}$ and
$F_{\rm X}$ (Wood et al., 2021), we expect to see traces of a two-part trend
(albeit with significant scatter) between $\dot{M}_{\bigstar}$ and $R_{\rm
o}$: a flat or saturated part that is independent of stellar rotation ($R_{\rm
o}\lesssim 0.1$, rapidly rotating stars), and a power law showing that the
stellar wind mass loss rate decreases with increasing $R_{\rm o}$ ($R_{\rm
o}>0.1$, slowly rotating stars).

For stars in the unsaturated regime, we do see a trend in which
$\dot{M}_{\bigstar}$ increases with decreasing $R_{\rm o}$. The relationship
between $\dot{M}_{\bigstar}$ and $R_{\rm o}$ retrieved from our simulations is

$\log\dot{M}_{\bigstar}/R_{\bigstar}^{2}=(-1.13\pm 0.23)\log R_{\rm
o}+(0.50\pm 0.07).$ (7)

The majority of the $\dot{M}_{\bigstar}$ derived from observation appears to
follow the established relationship $\dot{M}_{\bigstar}$–$R_{\rm o}$, with
some scatter within the error range. We do, however, notice four outliers,
including three K stars and one G star. The K stars with the high
$\dot{M}_{\bigstar}$ correspond to the binary 70 Oph A (K0V) and 70 Oph B
(K5V). As for the $3^{\rm\scriptscriptstyle rd}$ K star and the G star, they
correspond to evolved stars: $\delta$ Eri (K0IV,
$\dot{M}_{\bigstar}~{}=~{}0.6~{}\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$,
$R{\rm o}\,\sim~{}21$) and DK UMA (G4III-I,
$\dot{M}_{\bigstar}~{}=~{}0.0077~{}\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$,
$R_{\rm o}$ $\sim$ 2.51).

We do not expect evolved stars to follow the same trend as unsaturated main
sequence stars because their winds might be generated from a different
mechanism (such as pulsations, see Vidotto 2021). As for 70 Oph A and B, we do
not have much insight into their eruptive activity levels in order to rule out
whether or not the $\dot{M}$ inferred from the astropsheric technique was
influenced by slingshot prominences or CME activity.

As can be seen in Fig. 6, our numerical results in this region are essentially
bracketed by the observations for which the $R_{\rm o}$ reaches larger values.
The largest Rossby number from our star sample corresponds to HD 219134 (K3V,
$R_{\rm o}=2.02732$), which is comparable to the accepted solar value. Since
our models use ZDI maps as inner boundary conditions to simulate stellar
winds, this implies that extending our numerical models to even larger $R_{\rm
o}$ would be very challenging as those ZDI reconstructions would require
prohibitively long observing campaigns.

While we have limited data points, we see that for objects with $R_{\rm
o}\lesssim$ 0.15, we do not obtain larger numerical values
$\dot{M}_{\bigstar}$ even when the magnetic field strengths increase
dramatically. For example, in the case of YZ CMi ($B^{\rm max}_{\rm R}$ = 822
G, $\dot{M}_{\bigstar}/R_{\bigstar}^{\rm 2}$ = 20.57
$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$) and GJ 1245 B ($B^{\rm max}_{\rm R}$
= 404 G, $\dot{M}_{\bigstar}$ = 9.27 $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm
2}$). All stars on the left-hand side of Fig. 6 lie beneath the maximum
$\dot{M}_{\bigstar}$ value obtained for YZ CMi ($B^{\rm max}_{\rm R}$ = 822 G,
$\dot{M}_{\bigstar}/R_{\bigstar}^{\rm 2}$ = 20.57
$\dot{M}_{\rm\odot}/R_{\odot}^{\rm 2}$). This is true even when $R_{\rm o}$
varies by more than one dex, magnetic field strength by factors of 100, and
the expected complexity number by $\sim$ 4.

These results indicate that the contribution from the steady wind will only
account for a small fraction of the $\dot{M}_{\bigstar}$ budget in the case of
very active stars. Furthermore, the obtained behaviour hints of a possible
saturation of the steady-state stellar wind contribution to
$\dot{M}_{\bigstar}$, while the star could still lose significant mass through
other mechanisms such as slingshot prominences or CME activity due to flares
among others.

According to Villarreal D’Angelo et al. (2017) and references therein, cool
stars can support prominences if their magnetospheres are within the
centrifugal regime (i.e. $R_{\rm K}<R_{\rm A}$, where $R_{\rm
K}=\sqrt[3]{GM_{\bigstar}/\Omega_{\bigstar}^{2}}$ is the co-rotation radius).
They provide estimates for the prominence masses ($m_{\rm p}$) and the
ejection time-scales ($t_{\rm p}$) for a sample of cool stars. According to
their analysis, DT Vir would have $m_{\rm p}=1.5\times 10^{\rm 15}$ g and
$t_{\rm p}$ = 0.1 d, while the values for GJ 1245 B would be $m_{\rm
p}=4.4\times 10^{\rm 14}$ g, $t_{\rm p}$ = 0.3 d. Using these values, they
also reported the expected mass loss rate from prominences for these two stars
in absolute units. In order to compare with the steady state wind, we convert
their results to units of $\dot{M}_{\odot}$/$R_{\odot}^{\rm 2}$. For DT Vir we
have $\dot{M}_{\bigstar}^{\rm p}$/$R_{\bigstar}^{\rm 2}$ = 0.49
$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$ and for GJ 1245 B the resulting value
is $\dot{M}_{\bigstar}^{\rm p}$/$R_{\bigstar}^{\rm 2}$ = 0.68
$\dot{M}_{\odot}$/$R_{\odot}^{\rm 2}$.

For the CMEs contribution, we can obtain an order of magnitude estimate by
following the approach in Odert et al. (2017). They estimate the mass-loss
rate from the CME ($\dot{M}_{\bigstar}^{\rm CME}$) as a function of $L_{\rm
X}$ and the power law index ($\alpha$) of the flare frequency distribution.
For the X-ray luminosity, we used the Liefke & Schmitt (2007) database, and
for the flare frequency distribution exponent we took $\alpha$ = 2 (Hawley et
al., 2014).

For DT Vir, with log(LX) = 29.75, we obtain $\dot{M}_{\bigstar}^{\rm
CME}$/$R_{\bigstar}^{\rm 2}$ $\sim$ 160 $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm
2}$. For GJ 1245 B, with log(LX) = 27.47, the estimated CME-mass loss rate is
$\dot{M}_{\bigstar}^{\rm CME}$/$R_{\bigstar}^{\rm 2}$ $\sim$
12.8$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$.

We emphasize here that this approach assumes that the solar flare-CME
association rate holds for very active stars (see the discussion in Drake et
al. 2013). As such, it does not consider the expected influence due to CME
magnetic confinement (e.g. Alvarado-Gómez et al. 2018, 2019) which currently
provides the most suitable framework to understand the observed properties of
stellar CME events and candidates (Moschou et al., 2019; Alvarado-Gómez et
al., 2022a; Leitzinger & Odert, 2022).

Still, we can clearly see that the input from CMEs to the total
$\dot{M}_{\bigstar}$ could be higher than the steady wind and prominences for
these two stars (with the latter contributing less in these cases). For
instance, the estimated contribution of CMEs to the total $\dot{M}_{\bigstar}$
of DT Vir is almost 40 times higher than the value obtained for the steady
stellar wind.

We will discuss the cases of EV Lac and YZ CMi in Section 3.3.2

####  3.3.1 Comparison between simulations and observations

In addition to analyzing the general trends, we can compare the models for
common stars between our sample and the observations in Wood et al. (2021) and
references therein. The stars in Wood et al. (2021) contain a total number of
37 stars with a mix of main-sequence and evolved stars. The sample includes 15
single K-G stars among them 4 evolved stars, and 4 binaries. Wood et al.
(2021) reports individual $\dot{M}$ values for the G-K binary pairs (this
means that it was possible to model their individual contribution to the
astrosphere of the system or they were separated enough not to share a common
astrosphere). This is important as, in principle, one could treat the binary
pairs as individual stars. The rest of the star sample includes 22 M-dwarfs
with 18 single M-dwarfs, 3 binaries, and 1 triple system. Unlike the G-K
stars, $\dot{M}_{\bigstar}$ values for the M-dwarf binaries/triple system are
listed as a single value (therefore, it means that it has to be taken as the
aggregate of all the stars in the system). For the binary system GJ 338 AB we
were unable to include it in the plot of Fig. 6 due to a lack of needed
information to estimate its $R_{\rm o}$.

Following on the results from Sect. 3.3, our simulated mass loss rates for
stars in the unsaturated regime agree well with those estimated from
astrospheric detections (see Fig. 6). Specifically, for GJ 205 (M1.5V), 61 Cyg
A (K5V), and HD 219134 (K3V) we obtain
$\dot{M}_{\rm\bigstar}/R_{\bigstar}^{\rm 2}$ of 2.32, 3.98, and 1.50,
respectively. These values are all consistent with their respective
observational estimates, taking into account the typical uncertainties of the
astrospheric absorption method777Astrospheric estimates on
$\dot{M}_{\bigstar}$ should have an accuracy of about a factor of 2 with
substantial systematic uncertainties (Wood et al., 2005b).. While further
observations could help to confirm this, the agreement between our
asynchronous models and the observations indicates that, within this $R_{\rm
o}$ range, the temporal variability of $\dot{M}_{\bigstar}$ is minimal. This
is certainly the case for the Sun ($R_{\rm o}\sim 2.0$) in which long-term
monitoring has revealed only minor variability of the solar wind mass loss
rate over the course of the magnetic cycle (Cohen 2011, Finley et al. 2018,
2019).

On the other hand, $\dot{M}_{\bigstar}$ from the 3D MHD simulations appear to
fall short by an order of magnitude or more from the available estimates for
$\epsilon$ Eri (K2V), EV Lac (M3.5V), and YZ CMi (M4.5V) with
$\dot{M}_{\bigstar}/R_{\bigstar}^{\rm 2}$ of 4.53, 0.62 and 20.57,
respectively. We will discuss different possibilities for these discrepancies
on each star in Sect. 3.3.2. However, it is important to remember that the
$\dot{M}_{\bigstar}$ estimates from the Ly-$\alpha$ absorption technique
contain systematic errors that are not easily quantified. One example is that
they depend on the assumed properties and topology of the ISM (Linsky & Wood,
2014), which have not been fully agreed upon in the literature (e.g.,
Koutroumpa et al. 2009; Gry & Jenkins 2014; Redfield & Linsky 2015). While
studies have provided a detailed characterization of the local ISM (see
Redfield & Linsky 2008a, 2015; Gry & Jenkins 2014), intrinsic uncertainties
and additional observational limitations can greatly alter the estimated mass-
loss rate values. These include column densities, kinematics, and metal
depletion rates (Redfield & Linsky 2004a, 2008b), as well as local
temperatures and turbulent velocities (Redfield & Linsky, 2004b).

Furthermore, we would also like to emphasize the variation of the
$\dot{M}_{\bigstar}$ in the astrospheric estimates with the assumed stellar
wind velocity, as we believe that this factor is one of the largest potential
source of uncertainty and discrepancy with our models. As discussed by Wood et
al. (2021), this parameter is used as input in 2.5D hydrodynamic models to
quantify the stellar wind mass loss rate. The Ly-$\alpha$ absorption
signature, leading to $\dot{M}_{\bigstar}$, is determined to first order by
the size of the astrosphere. The latter depends on the stellar wind dynamic
pressure ($P_{\rm dyn}\propto\dot{M}_{\bigstar}\,U_{\rm sw}$), which implies
an inverse relation between $\dot{M}_{\bigstar}$ and $U_{\rm sw}$ (Wood et
al., 2002).

The astrospheric analysis of Wood et al. (2021) assumed a stellar wind
velocity of 450 km s-1 at $1$ au (matching models of the heliosphere) for all
main-sequence stars. However, we find that stellar wind velocities can vary
significantly between different types of stars and even among the same
spectral type for different magnetic field strengths and rotation periods. To
quantify this, we compute the average terminal velocity of the wind,
($\overline{{U_{\rm R}^{\rm T}}}$), by averaging $U_{\rm R}$ over a sphere
extracted at 99$\%$ of the maximum extent of each simulation domain (594
$R_{\rm\bigstar}$ for F, G, and K stars and 248 $R_{\rm\bigstar}$ for
M-dwarfs; see Sect. 2). In the cases in which the spatial extension of our
numerical domain allowed, we also computed the average wind velocity at $1$
au. The resulting values, listed in Table 2, indicate variations in the wind
velocity by factors of 5 or more when moving from F-type stars
($\overline{U_{\rm R}^{\rm T}}\sim 325$ km s-1) to M-dwarf ($\overline{U_{\rm
R}^{\rm T}}\sim 1500$ km s-1). This is also illustrated in Fig. 7, which
portrays the simulated stellar wind environment for HD 179949 (F8V), HD 73256
(G8V), HD 189733 (K2V), and DT Vir (M0V). We include a green iso-surface that
corresponds to the wind velocity at $1$ au for F, G, and K stars as for
M-dwarfs it represents the average terminal wind velocity in the domain. The
visualizations also include the equatorial projection of the wind dynamic
pressure ($P_{\rm dyn}$ = $\rho U^{\rm 2}$), normalized to the nominal Sun-
Earth value, as well as on a sphere highlighting the wind 3D structure at
$0.5$ au.

What is clear from this analysis is that is not ideal to use the same wind
velocity for all spectral types. Even within the same spectral type, we can
observe a wide range of terminal velocities (e.g., the velocity in K stars
ranges from $400$ km s-1 to $700$ km s-1). As such, for models that require
wind velocity as an input parameter, we recommend using the average radial
wind velocity among a given spectral type.

For G-K stars, we obtain wind velocities at $1$ au in the range of 400 to
$700$ km s-1 which is not too different from the wind velocity assumption of
Wood et al. 2021. This is also consistent with the fact that for these
spectral types, we have a better agreement between $\dot{M}$ estimated from
our simulations and those from the astropsheric technique (Wood et al., 2021).
For lower mass stars with relatively small $R_{\rm o}$ we obtain velocities
higher than $450$ km s-1 up to $3675$ km s-1.

Note that due to computational limitations, the extent of our M-dwarf
simulations does not reach up to $1$ au (varying from $0.6$ au for DT Vir to
$0.16$ au for GJ 1245 B). Nevertheless, as indicated by the calculated
terminal velocities, even at closer distances the wind velocity is already
$>450$ km s-1, a situation that should still hold when propagated out to $1$
au. Wind velocities on the order of $1000-1500$ km s-1 at distances of $1$ au
and beyond had been reported in high-resolution AWSoM simulations of the
environment around the M5.5V star Proxima Centauri (Alvarado-Gómez et al.,
2020). This helps to explain why our simulated mass-loss rates for EV Lac, YZ
CMi, and $\epsilon$ Eri were lower than the observed ones (differences larger
than a factor of 2). We discuss these cases in more detail in the following
section.

####  3.3.2 Exploring the cases of EV Lac, YZ CMi, $\epsilon$ Eri

  1. 1.

YZ CMi & EV Lac

Frequent stellar flares have been observed at YZ CMi in several wavelength
ranges (Lacy et al. 1976; Mitra-Kraev et al. 2005; Kowalski et al. 2013; Bicz
et al. 2022). The flaring energy distribution of this star ranges from
$10^{30.6}$ to $10^{34.09}$ erg (Bicz et al., 2022) with a total flaring time
that varies from 21 to 306 minutes. Likewise, there is also significant flare
activity on EV Lac (Leto et al. 1997; Muheki et al. 2020). From spectroscopic
and photometric studies of EV Lac, Muheki et al. (2020) reports to have found
27 flares ($\sim$ 5.0 flares per day) in H $\alpha$ with energies between
$1.61\times 10^{31}$ erg $-1.37\times 10^{32}$ erg and 49 flares ($\sim$ 2.6
flares per day) from the TESS lightcurve with energies of $6.32\times 10^{31}$
erg $-1.11\times 10^{33}$ erg. With such high flare activity, it is possible
that a large fraction of the $\dot{M}_{\bigstar}$ estimated in Wood et al.
(2021) for these stars could arise from transient phenomena (e.g.,
prominences, CMEs).

Following the same approach described at the end of Section 3.3.1, we can
obtain a rough estimate of $\dot{M}$ from CMEs for EV Lac and YZ CMi. For EV
Lac we find $\dot{M}_{\bigstar}^{\rm CME}$/$R_{\bigstar}^{\rm 2}$ = 55.5
$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$ assuming log(LX) = 28.69. In the case
of YZ CMi, an log($L_{\rm X}$) = 28.53 yields $\dot{M}_{\bigstar}^{\rm
CME}$/$R_{\bigstar}^{\rm 2}$ = 47.6 $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$.
However, given the magnetic field strength observed in EV Lac and YZ CMi (a
few kG, Shulyak et al. 2010; Reiners et al. 2022; Cristofari et al. 2023), we
expect that the magnetic confinement of CMEs would play an important role in
these objects (see Alvarado-Gómez et al. 2019, Muheki et al. 2020, Maehara et
al. 2021). Therefore, it is not straightforward to estimate exactly how large
the contribution of CMEs to $\dot{M}_{\bigstar}$ is for these stars.

In addition, as discussed by Villarreal D’Angelo et al. (2017), EV Lac and YZ
CMi are considered in the slingshot prominence regime. For EV Lac they
estimate $m_{\rm p}=2.0\times 10^{\rm 16}$ g and $t_{\rm p}$ = 0.6 d, while
for YZ CMi values of $m_{\rm p}=4.5\times 10^{\rm 16}$ g and $t_{\rm p}$ = 0.6
d are given. Using the associated mass loss rate values reported in Villarreal
D’Angelo et al. (2017), we obtain $\dot{M}_{\bigstar}^{\rm
p}$/$R_{\bigstar}^{\rm 2}$ = 3.16 $\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$ for
EV Lac and $\dot{M}_{\bigstar}^{\rm p}$/$R_{\bigstar}^{\rm 2}$ = 8.32
$\dot{M}_{\rm\odot}$/$R_{\odot}^{\rm 2}$ for YZ CMi.

This suggests another possible explanation for the discrepancies between our
models and the astrospheric estimates is that some of the stellar wind
detected for EV Lac and YZ CMi contains material from the slingshot
prominences. Indeed, the location of the latter in the $\dot{M}_{\bigstar}$ –
$R_{\rm o}$ diagram (Fig. 6) appears more consistent with the mass loss rate
estimates from slingshot prominences by Jardine & Collier Cameron (2019).

Moreover, Wood et al. (2021) noted that the YZ CMi astrospheric absorption
comes primarily from neutrals near and inside the astropause, rather than from
the hydrogen wall where neutral H density is highest. Therefore, using Ly
alpha absorption to calculate $\dot{M}_{\bigstar}$ from YZ CMi will result in
substantial uncertainty.

Finally, as mentioned in Sect. 3.3.1, there is a significant difference
between the wind velocity assumed by Wood et al. (2021) and our results. Our
average terminal wind velocity for YZ CMi ($1709$ km s-1) and EV Lac ($3675$
km s-1) is significantly higher than the wind velocity of $450$ km s-1 assumed
in Wood et al. (2021) at 1 au. While the wind velocity in EV Lac might be
overestimated in our models (due to the usage of fiducial AWSoM parameters),
we still expect relatively large wind velocities for this star ($\sim$
$1000-1500$ km s-1) given its magnetic field strength and Rossby number (see
e.g., Kavanagh et al. 2021; Alvarado-Gómez et al. 2022b). As was discussed in
Sect. 3.3.1, while our terminal wind velocity for M-dwarfs is calculated
closer to the star ($0.33$ au for YZ CMi and $0.16$ au for EV Lac), we do not
expect a large reduction in the average velocity between these distances and 1
au. As such, the fast wind velocity resulting in our simulations of YZ CMi and
EV Lac would imply lower $\dot{M}_{\bigstar}$ values when analyzed following
the astrospheric technique of Wood et al. (2021).

  2. 2.

$\epsilon$ Eri

Figure 8:  Numerical results of the average Alfvén surface size (diamonds),
the inner and outer edge of the HZ (square, Kopparapu et al. 2014) derived
using scaling laws that connect $M_{\bigstar}$, $R_{\bigstar}$,
$L_{\bigstar}$, and $T_{\rm eff}$, against $M_{\bigstar}$. The filled diamonds
correspond to the $\overline{AS}$ of the 21 stars in our sample. The empty
diamonds are the $\overline{AS}$ of other stars derived from our
$\overline{AS}$ \- $B_{\rm R}$ relation: log $\overline{AS_{\rm R}}=(0.42\pm
0.06)\log B_{\rm R}^{\rm avg}+(0.71\pm 0.07)$. The habitable zone boundaries
are color-coded by the corresponding average dynamic pressure
($\overline{P_{\rm dyn}}$) in logarithmic scale. The green shaded area
represents the optimistic habitable zone. Earth is represented by $\oplus$.
The purple circles correspond to some confirmed exoplanets (taken from the
NASA exoplanet archive). The dynamic pressure at the outer edge of the HZ of
DT Vir and GJ 205 (derived from scaling laws) is missing in the plot since it
goes beyond the simulation domain that was initially established using the
measured $L_{\bigstar}$ and $T_{\rm eff}$ (Table 1). The dashed lines separate
the different spectral types.

With a relatively slow rotation period ($11$ d), and weak large-scale magnetic
field ($<$ 50 G), $\epsilon$ Eri cannot be considered within the slingshot
prominence regime (like in the cases of YZ CMi and EV Lac). Because of this,
we do not expect a significant presence of slingshot prominences in the
$\dot{M}$ value of this star. On the other hand, the analysis of Loyd et al.
(2022), estimated the contribution of flare-associated CMEs to the mass loss
rate. They reported an upper limit of $1.09$ $\dot{M}_{\rm\odot}$
$/R_{\odot}^{\rm 2}$, which is insignificant when compared to the star’s
overall estimated $\dot{M}_{\bigstar}$ value by Wood et al. (2021) and the
astrospheric technique (56 $\dot{M}_{\rm\odot}$ $/R_{\odot}^{\rm 2}$).
Therefore, the contribution from CMEs is also most likely not responsible for
the elevated astrospheric $\dot{M}_{\bigstar}$ value on this star and its
discrepancy with our steady-state models.

On the other hand, multiple observations of the large-scale magnetic field
geometry of $\epsilon$ Eri reveal that it evolves over a time-scale of months
(Jeffers et al., 2014; Jeffers et al., 2017). According to Jeffers et al.
(2014), the maximum field strength can reach up to $42$ G. As shown in Fig. 3,
a global increase in the magnetic field strength causes an increase in
$\dot{M}_{\bigstar}$. The Zeeman Doppler Imaging map of $\epsilon$ Eri used to
drive the 3D MHD model has a $B_{\rm R}^{\rm max}$ = 25 G leading to
$\dot{M}_{\bigstar}/R_{\bigstar}^{2}=4.53~{}\dot{M}_{\rm\odot}/R_{\rm\odot}^{2}$.
This value is comparable to the numerical result obtained by Alvarado-Gómez et
al. (2016) for this star ($\dot{M}_{\bigstar}/R_{\bigstar}^{2}\sim
5.3\dot{M}_{\rm\odot}/R_{\rm\odot}^{\rm 2}$). Increasing the surface magnetic
field strength of $\epsilon$ Eri to the maximum value reported in observations
will raise the mass loss rate to $\sim
10~{}\dot{M}_{\rm\odot}/R_{\rm\odot}^{\rm 2}$. As such, the variability of the
stellar magnetic field and its expected modulation of the stellar wind
properties could account for some of the differences between the simulated and
the observed mass loss rates. However, corroborating this would require
contemporaneous ZDI and astrospheric measurements which, to our knowledge,
have not been performed on any star so far. As $\epsilon$ Eri goes through a
magnetic/activity cycle (Metcalfe et al. 2013; Jeffers et al. 2017), we can
expect relatively large variations in $\dot{M}_{\bigstar}$ values in our
Alfvén-wave driven stellar wind models.

Finally and following the discussion for YZ CMi and EV Lac, the average wind
velocity for $\epsilon$ Eri at $1$ au ($554$ km s-1) resulting from our models
exceeds the one assumed in Wood et al. (2021). This will result in a smaller
estimated $\dot{M}_{\bigstar}$ value from the pressure-balance astrospheric
technique. In this way, the deviation between our models and the astrospheric
detection of $\epsilon$ Eri could be due to the combined contribution of all
the preceding elements (i.e., CMEs, cycle-related variability of the magnetic
field, higher stellar wind velocity), and therefore we do not consider this
discrepancy critical to our analysis.

###  3.4 Stellar wind and Circumstellar region

This section focuses on using the stellar wind results obtained from the 3D
MHD simulations to assess the conditions an exoplanet would experience. This
includes the characterization of the Alfvén surface for the various stellar
wind solutions, the properties of the stellar wind in the habitable zone of
these stars (in terms of the dynamical pressure of the wind), and the
resulting magnetosphere size for these stellar wind conditions (assuming that
a planet with the same properties/magnetization as Earth is in the HZ of these
stars). The obtained quantities are listed in Table 2 and 3.

####  3.4.1 Stellar wind properties and orbital distances

  1. 1.

Alfvén surface size

Figure 8 summarizes our results showing the stellar wind environment around
cool main sequence stars. We include the average size of the AS, resulting
from our 21 3D MHD models, indicated in filled diamonds. To complement this
information, empty diamonds correspond to the expected average AS size
employing the scaling relation provided in Sect. 3.1, and using the ZDI
information from 29 additional stars (See et al. 2019c and reference therein).
The green region corresponds to the optimistic HZ, calculated using the
approach provided by Kopparapu et al. (2014) and the expected behaviour of the
luminosity, temperature as a function of stellar mass on the main sequence
(Kasting et al. 1993; Kopparapu et al. 2014; Ramirez 2018). Each square
indicates the limits of the optimistic HZ for each star in our sample. These
have been color-coded by the stellar wind dynamic pressure, normalized to the
average Sun-Earth value. The position of the Earth is indicated by the
$\oplus$ symbol. In the background, a sample of the semi-major axis of some
exoplanets is included.

There are a few noteworthy aspects of Fig. 8. First of all, the 3D MHD
simulated $\overline{AS}$ values (filled diamonds) do not show a clear trend
with stellar mass. Instead, we see more or less similar $\overline{AS}$
regardless of the spectral type of the star (Table 2). We see a similar
behavior for stars whose $\overline{AS}$ were extracted from the scaling
relationship presented in Eq. 3 (empty diamonds). There is a significant
scatter in the obtained distribution of $\overline{AS}$ against
$M_{\bigstar}$, indicating that the intrinsic dependency with the surface
magnetic field properties can in principle be replicated among multiple
spectral types. However, we remind the reader that this result is also partly
a consequence of our fixed choice for the base parameters of the corona and
stellar wind solution (Sect. 2.1.3), which could in principle vary among
different spectral types and activity stages (i.e. ages). As such, the
generalization of the results presented here requires further investigation
from both, observational constraints and numerical simulations.

We can also see that for late K and M-dwarfs, $\overline{AS}$ reaches orbital
distances comparable to their HZ limits. Examples of this from our sample are
GJ 1245 B ($\overline{AS}$ = 0.028 au, $HZ_{\rm inner}$ = 0.033 au) and YZ CMi
($\overline{AS}$= 0.178 au, $HZ_{\rm inner}$= 0.09 au). This situation has
been also identified in previous case studies of stellar winds and exoplanets
(e.g, Vidotto et al. 2013; Cohen et al. 2014; Garraffo et al. 2017).

The location of the HZ relative to the stellar Alfvén surface must be
considered when studying the interactions between a star and a planet. A
planet orbiting periodically or continuously within the $AS$ region could be
directly magnetically connected to the stellar corona, which could have
catastrophic effects on atmospheric conservation (Cohen et al., 2014; Garraffo
et al., 2017; Strugarek, 2021). On the other hand, a planet with an orbit far
outside this limit will be decoupled from the coronal magnetic field and
interact with the stellar wind in a manner similar to the Earth (e.g. Cohen et
al. 2020; Alvarado-Gómez et al. 2020). In the case of a planet orbiting in and
out of the AS, the planet will experience strongly varying wind conditions,
whose magnetospheric/atmospheric influence will be greatly mediated by the
typical time-scale of the transition (Harbach et al., 2021).

Table 3: Numerical results of different parameters in our sample at the
habitable zone. Columns 1–8, respectively, list the star name, average dynamic
pressure at the middle of the HZ ($\overline{P}_{\rm dyn,HZ}$), average
dynamic pressure at the inner boundary of the HZ ($\overline{P_{\rm dyn}^{\rm
Inn,HZ}}$), average dynamic pressure at the outer boundary of the HZ
($\overline{P_{\rm dyn}^{\rm Out,HZ}}$), the average magnetopause standoff
radius ($\overline{R}_{\rm M}$), the inner habitable zone ($HZ_{\rm inner}$),
the outer habitable zone ($HZ_{\rm outer}$), the average equatorial Alfvén
surface ($\overline{AS}_{\rm eq}$). The habitable zones listed in this table
were inferred using the measured $L_{\bigstar}$ and $T_{\rm eff}$ of each star
in our sample (Table 1).

Name |  ${\overline{P}_{\rm{dyn,HZ}}}$ [$P_{\rm\oplus}$] |  $\overline{P_{\rm dyn}^{\rm Inn,HZ}}$ [$P_{\rm\oplus}$] |  $\overline{P_{\rm dyn}^{\rm Out,HZ}}$ [$P_{\rm\oplus}$] |  $\overline{R}_{\rm M,HZ}$ [$R_{\rm E}$] |  $HZ_{\rm inner}$ [au] |  $HZ_{\rm outer}$ [au] |  $\overline{AS}_{\rm eq}$ [$R_{\rm\bigstar}$]
---|---|---|---|---|---|---|---
$\tau$ Boo |  1.02 |  2.76 |  0.40 |  7.85 |  1.253 |  2.901 |  9.99
HD 179949 |  0.92 |  2.51 |  0.48 |  7.99 |  0.983 |  2.294 |  10.24
HD 35296 |  2.93 |  7.77 |  1.48 |  6.58 |  0.925 |  2.157 |  12.78
HN Peg |  6.74 |  18.58 |  3.47 |  5.73 |  0.812 |  1.904 |  13.65
HD 190771 |  3.96 |  11.01 |  2.02 |  6.26 |  0.744 |  1.748 |  12.35
TYC 1987-509-1 |  10.21 |  28.49 |  5.20 |  5.30 |  5.35 |  1.299 |  15.47
HD 73256 |  2.32 |  6.45 |  1.18 |  6.85 |  0.648 |  1.539 |  10.73
HD 130322 |  0.91 |  2.56 |  0.46 |  8.00 |  0.542 |  1.291 |  8.12
HD 6569 |  10.77 |  30.21 |  5.46 |  5.30 |  0.466 |  1.11 |  16.25
$\epsilon$ Eri |  6.45 |  18.16 |  3.23 |  5.77 |  0.426 |  1.024 |  13.64
HD 189733 |  13.56 |  38.82 |  6.88 |  5.10 |  0.458 |  1.108 |  13.4
HD 219134 |  2.01 |  5.81 |  1.01 |  7.01 |  0.410 |  0.996 |  8.33
TYC 6878-0195-1 |  13.42 |  39.39 |  6.71 |  5.11 |  0.713 |  1.751 |  23.16
61 Cyg A |  10.06 |  29.47 |  5.01 |  5.36 |  0.308 |  0.754 |  11.93
HIP 12545 |  20.28 |  60.51 |  9.98 |  4.77 |  0.509 |  1.267 |  24.45
TYC 6349-0200-1 |  13 |  40 |  6 |  5.14 |  0.443 |  1.116 |  21
DT Vir |  31.48 |  88.97 |  16.20 |  4.43 |  0.191 |  0.486 |  29.24
GJ 205 |  9.73 |  29.65 |  4.73 |  5.39 |  0.201 |  0.517 |  15.91
EV Lac |  12.15 |  33.18 |  6.27 |  —- |  0.094 |  0.245 |  119.16
YZ CMi |  193.85 |  566.34 |  97.08 |  —- |  0.09 |  0.237 |  75.23
GJ 1245 B |  140.89 |  447.39 |  68.15 |  —- |  0.033 |  0.087 |  38.12

Figure 9: Two-dimensional Mercator projections of the normalized stellar wind
dynamic pressure ($P_{\rm dyn}$) extracted from the 3D MHD models of four
stars in our sample covering F to M spectral types (HD 179949, TYC 198-509-1,
61 Cyg A, and GJ 205). Each $P_{\rm dyn}$ distribution was extracted from a
spherical surface located at the midpoint habitable zone of each star.

  2. 2.

Dynamic pressure

We also see a general trend in Fig. 8 in which the dynamic pressure at the HZ
boundaries increases as we move from earlier to later spectral types. For
example, $\overline{P_{\rm dyn}^{\rm Inn,HZ}}$ for the lowest-mass star GJ
1245 B is 447.39 $P_{\rm\oplus}$ nearly 200 times stronger than for the
highest-mass star $\tau$ Boo with $\overline{P_{\rm dyn}^{\rm Inn,HZ}}$ = 2.76
$P_{\rm\oplus}$. Our results also show a large variability in
$\overline{P_{\rm dyn}}$ as we move from the inner to the outer edge of the HZ
of G, K, and M-dwarfs (Table 3). For these stars, the $\overline{P{\rm dyn}}$
at the inner HZ is almost 6 times stronger than that at the outer edge of the
HZ (i.e., EV Lac $\overline{P_{\rm dyn}^{\rm Inn,HZ}}$ = 33.18
$P_{\rm\oplus}$, $\overline{P_{\rm dyn}^{\rm Out,HZ}}$ = 6.27
$P_{\rm\oplus}$). For F stars, the difference is smaller, around a factor of 2
like in the case of HD 179949, where $\overline{P_{\rm dyn}^{\rm Inn,HZ}}$ =
2.51 $P_{\rm\oplus}$ and $\overline{P_{\rm dyn}^{\rm Out,HZ}}$ = 0.48
$P_{\rm\oplus}$. The reason is that the HZs of these stars are farther from
the star, where the wind density starts to become less variable.

Moreover, in some cases, we have $\overline{P}_{\rm dyn}$ at the inner and
outer edge of the star HZ comparable to the typical range experienced by the
Earth (0.75 and 7 nPa, Ramstad & Barabash 2021). For example, HD 73256 (G8V,
6.45 - 1.18 $P_{\rm\oplus}$ $\sim$ 9.675 - 1.77 nPa), HD 130322 (K0V, 2.56 -
0.46 $P_{\rm\oplus}$ $\sim$ 3.84 - 0.69 nPa), $\tau$Boo (F7V, 2.76 - 0.40
$P_{\rm\oplus}$ $\sim$ 4.14 - 0.6 nPa). For the case of M-dwarfs, we have
dynamic pressures higher than those experienced by Earth, as in the case of DT
Vir (M0V, 88.97 - 16.20 $P_{\rm\oplus}$ $\sim$ 133.455 - 24.3 nPa). This is
because the HZ is located near the star where the density is highest. This
indicates that planets orbiting at very close distance to the star ($\sim$
0.03 - 0.05 au) would experience extreme space weather conditions with
${P}_{\rm dyn}$ up to $\rm 10^{3}$ and $\rm 10^{4}$ $P_{\rm\oplus}$. These
values are comparable to the ones estimated in Alvarado-Gómez et al. (2020)
for Proxima d and for Proxima b in Garraffo et al. (2016). However, the reader
is reminded here that any point from our simulations should be interpreted as
an indication of the average conditions, but should not be treated as a
specific absolute value (since it will change depending on the instantaneous
local density and velocity of the wind (both a function of the evolving
stellar magnetic field).

In addition, we notice a scatter in $\overline{P}_{\rm dyn}$ estimates at the
HZ when comparing stars of the same spectral type. This is not surprising
since the $P_{\rm dyn}$ depends on the wind velocity and density at a given
place. This also translates into having a range of dynamic pressure that a
planet will experience within the HZ. This will defer from one orbital
distance to the other as we can see in Fig. 7 where we show the equatorial
plane color-coded by the dynamic pressure.

We can use our 3D models to investigate also the influence due to the orbital
inclination. To illustrate this, Fig. 9 shows a 2D projection of the
normalized dynamic pressure $P_{\rm dyn}$ extracted from spherical surfaces
matching the midpoint of the HZ of HD 179949 (F8V), TYC 198-509-1 (G7V), 61
Cyg A (K6V), and GJ 205 (M1.5V). We notice that in the case of F and G stars
(i.e., HD 179949, and TYC-198-509-1) we have a large $P_{\rm dyn}$ variation
with inclination around a factor 7. However, $P_{\rm dyn}$ values, are still
relatively small in terms of absolute units (i.e., 0.01 - 10 $P_{\rm\oplus}$
$\sim$ 0.015 - 15 nPa). For K and M-dwarfs, we see less variability in the
$P_{\rm dyn}$ for the different inclinations, a more homogeneous $P_{\rm
dyn}$, especially in the case of the K star. However, in these cases, the
$P_{\rm dyn}$ can reach values > 100 $P_{\rm\oplus}$ (> 150 nPa). Our results
also show that even with an extreme orbit around the G-type star (TYC
198-509-1) with an inclination matching the current sheet, we would most
likely not reach the very high $P_{\rm dyn}$ values as in the case of the K
and M-dwarfs as we move closer to the star. As such, the inclination of the
orbit plays a secondary role compared to the distance. This is clearly seen in
the color gradient that gets redder and redder as we move toward lower masses
(so the HZ is closer).

On the other hand, the variability of $P_{\rm dyn}$, which we can see in Fig.
8 while represented in the same ’spatial scale’, it does not coincide in terms
of ’temporal scales’. In other words, the $x$-axis in Fig. 8 do not correspond
to the same timescale units for each star, where the 360 degrees of longitude
correspond to “1 orbital period”. However, the orbital period is very
different for a planet in the HZ of an F-type star (within a few au) compared
to a planet orbiting an M-dwarf (within a fraction of an au). A planet
orbiting an M-dwarf star experiences the variations in $P_{\rm dyn}$ on a much
faster timescale ($\sim$ 1 day for each current sheet crossing), while these
variations are much longer for more massive stars. This means that even if the
$P_{\rm dyn}$ values were the same, the faster variability over the orbital
period for low-mass stars would result in planets and their
magnetospheres/atmospheres having less time to recover from passing through
regions of high $P_{\rm dyn}$ than planets around more massive stars.

Finally, following the results compiled by Ramstad & Barabash (2021), if we
consider the presence of a rocky exoplanet with an atmosphere similar to those
of Venus and Mars at those mid-HZ locations, we would expect atmospheric ion
losses between $2\times 10^{24}$ ions s-1 and $5\times 10^{24}$ ions s-1. This
of course assumes that all processes occur in the same way as in the solar
system (which might not be necessarily true for some regions of the vast
parameter space of this problem). The ion losses will depend heavily on the
type of stars that the exoplanet orbits, both in terms of the high-energy
spectra and the properties of the stellar wind (see e.g. Egan et al. 2019a;
France et al. 2020). If the rocky exoplanet is found around the HZ of an
M-dwarf, the planet might suffer from unstable stellar wind conditions as
previously stated that might increase the ion losses in the exoplanetary
atmosphere. We will consider the case of Earth with its magnetosphere in the
following section.

####  3.4.2 Magnetopause Standoff Distances

Using the dynamic pressure, we can define a first-order approximation to
determine the magnetosphere standoff distance ($R_{\rm M}$) of a hypothetical
Earth-like planet orbiting at the HZ around each star in our sample. This is
done by considering the balance between the stellar wind dynamic pressure and
the planetary magnetic pressure (Eq. 8, Gombosi 2004; Shields et al. 2016):

$R_{\rm M}=R_{\rm E}[\frac{{B_{\rm p}}^{2}}{8\pi P_{\rm dyn}}]^{\frac{1}{6}}$
(8)

The Earth’s equatorial dipole field and radius are represented by $B_{\rm p}$
and $R_{\rm E}$ respectively. Normally the total wind pressure should be
considered (i.e., thermal, dynamic, and magnetic), but in all the cases here
considered, we can neglect the contributions of the magnetic and thermal
pressures. For this calculation, we assume an equatorial dipole magnetic field
of 0.3 G, similar to that of the Earth (Pulkkinen, 2007). The magnetospheric
standoff distance is expressed in Earth’s radii (Eq. 8). The different
$\overline{R}_{\rm\scriptscriptstyle M,HZ}$ values for the different stars in
our sample are listed in table 3. Note that we only estimate the $R_{\rm M}$
in the cases where the HZ is in the super-Alfvénic regime (Cohen et al. 2014;
Strugarek 2018).

Our estimated $R_{\rm M,HZ}$ for F, G, and early K stars have values closer to
the standard size of Earth’s dayside magnetosphere ($\sim$ 10 $R_{\rm\oplus}$,
see Pulkkinen 2007; Lugaz et al. 2015). This is comparable to the value
obtained by Alvarado-Gómez et al. (2020) for Proxima c ($\sim$ 6 - 8
$R_{\rm\oplus}$ in both activity levels), assuming an Earth-like dipole field
on the planet surface. For the late K and M-dwarfs in our star sample, $R_{\rm
M}$ starts to reach lower values < 50$\%$ from that of Earth. This suggests
that a planet orbiting these stars must have a stronger dipole magnetic field
than that of the Earth to withstand the wind conditions since $R_{\rm
M}\propto B_{\rm p}^{1/3}$. However, in Ramstad & Barabash (2021) they show
that contrary to what we have seen so far, the magnetosphere might actually
not act as a shield for the stellar wind-driven escape of planetary
atmospheres. In fact, they reported an ion loss for Earth that ranges from 6
$\times 10^{24}$ ions s-1 \- 6 $\times 10^{26}$ ions s-1 which is higher than
what Venus and Mars lose. Further modeling studies are needed in order to
characterize the stellar wind influence on the atmospheric loss of rocky
exoplanets (e.g., Egan et al. 2019b; France et al. 2020), whose input stellar
wind parameters can be extracted from this investigation.

##  4 Summary & Conclusions

In this study we employed a state-of-the-art 3D MHD model (SWMF/AWSoM) to
investigate the dependencies between different star properties
($R_{\bigstar}$, $M_{\bigstar}$, $B_{\rm R}$, and $P_{\rm rot}$) and a number
of stellar wind parameters (AS, $\dot{M}_{\bigstar}$, $\dot{J}_{\bigstar}$,
$P_{\rm dyn}$) of cool main sequence stars. We present numerical results of 21
stars going from F to M stars with magnetic field strengths between 5 and 1.5
kG and rotation periods between 0.71 d and 42.2 d. The large-scale magnetic
field distribution of these stars, obtained by previous ZDI studies, were used
to drive the solutions in the Stellar Corona domain, which are then self-
consistently coupled for a combined solution in the Inner Astrosphere domain
in the case of F, G, and K stars.

Our results showed a correlation between the average AS size and $B_{\rm
R}^{\rm avg}$, regardless of the spectral type of the star (Eq. 3). We also
obtained a strong correlation between $\dot{M}_{\bigstar}$ and $B_{\rm R}^{\rm
avg}$ for the different spectral types (excluding EV Lac, Eq. 4). The
correlation between $\dot{J}_{\bigstar}$ and $B_{\rm R}$, on the other hand,
was dominated by the absolute dependence on the stellar size, with significant
scatter resulting mainly from the variability in $\dot{M}_{\bigstar}$, the
distribution of $\Omega_{\bigstar}$ in our sample and the equatorial AS size
where the maximum torque is applied.

Having established these star-wind relations, we looked in detail at
$\dot{M}_{\bigstar}$, since it is the only observable parameter of the stellar
wind for which comparisons can be made. Using the complexity number as a
function of the Rossby number $R_{\rm o}$–defined previously in the
literature– we were able to investigate the dependence of magnetic complexity
on $\dot{M}_{\bigstar}$. Our results showed that for more active stars, as in
the case of M-dwarfs, the field strength starts to dominate over the
complexity in the contribution on shaping $\dot{M}_{\bigstar}$. Also, for
cases in which the magnetic field strength and complexity were comparable, we
obtained similar $\dot{M}_{\bigstar}$. This indicates that in these cases the
stellar properties ($R_{\bigstar}$, $M_{\bigstar}$, and $P_{\rm rot}$) play a
secondary role in changing $\dot{M}_{\bigstar}$.

We then used our stellar wind results to investigate its behaviour with
respect to the well-known stellar activity relationship ($F_{\rm X}$ vs
$R_{\rm o}$ with the saturated and unsaturated regimes). For stars in the
unsaturated regime, we see a trend where $\dot{M}_{\bigstar}$ increases with
decreasing $R_{\rm o}$ (Eq. 7). For stars in the saturated regime, we find
that the contribution of the steady wind is only a small part of the
$\dot{M}_{\bigstar}$ budget. This suggests that there could be saturation in
$\dot{M}_{\bigstar}$ due to the steady stellar wind, while the star could lose
even more mass through other mechanisms, such as transient events (i.e.
prominences, coronal mass ejections).

In addition to analyzing the general trends, we compared the model results of
stars in our sample and objects with astrospheric $\dot{M}_{\bigstar}$
constraints. Our simulated $\dot{M}_{\bigstar}$ for stars in the unsaturated
regime agree well with those estimated from astrospheric detections (namely
for GJ 205, 61 Cyg A, and HD 219134). On the other hand, $\dot{M}_{\bigstar}$
from the 3D MHD simulations appear to differ by an order of magnitude or more
from available estimates for $\epsilon$ Eri, EV Lac, and YZ CMi. We discussed
how these results might be connected with the underlying assumption made by
the observational analysis with respect to the stellar wind speed. Indeed, for
all the stars in which our models differed largely from the literature
estimates, we obtained much larger stellar wind speeds than the ones used in
the astrospheric method. As such, we emphasized the importance of using the
appropriate wind velocity when estimating $\dot{M}_{\bigstar}$ from
observations.

We further discussed various possibilities for the discrepancies in EV Lac, YZ
Cmi, $\epsilon$ Eri. For the two flaring stars, EV Lac and YZ CMi, we suspect
that the high $\dot{M}_{\bigstar}$ estimates from the Ly-$\alpha$ absorption
technique could be dominated by material from slingshot prominences and
possibly CMEs (uncertain due to the expected magnetic confinement of CMEs in
these stars). Note that this possibility was also considered by Wood et al.
(2021) in the original astrospheric analysis. In the case of $\epsilon$ Eri,
we do not expect a large contribution from prominences or CMEs to the observed
$\dot{M}_{\bigstar}$. However, as $\epsilon$ Eri undergoes a magnetic cycle,
the stellar magnetic field and its expected modulation of stellar wind
properties could explain some of the differences between the simulated and
observed $\dot{M}_{\bigstar}$.

Moreover, we used the stellar wind results from the 3D MHD simulations to
assess the conditions that an exoplanet would experience, and provide the
stellar wind conditions in the entire classical Habitable Zones of our target
stars. Our results show a scatter in the obtained distribution of AS versus
$M_{\bigstar}$, suggesting that the intrinsic dependence with the surface
magnetic field properties can be reproduced for several spectral types. With
respect to the stellar wind dynamic pressure, our results show that the
orbital inclination plays a secondary role compared to the orbital distance.
We have also found that a planet orbiting K and M stars must have a stronger
dipole magnetic field than that of Earth to withstand the wind conditions, if
the planetary magnetic field is indeed acting as a shield (this paradigm,
however, is starting to be challenged by solar system observations).

Finally, the properties of the stellar wind in the HZ of different spectral
types obtained here can be used in future studies to, for instance, estimate
the expected radio emission due to wind-magnetosphere interactions or the
planetary atmospheric mass loss due to erosion of the stellar wind from ion
escape processes.

## Acknowledgements

The authors would like to thank the referee for valuable comments that
improved the quality of the paper. The authors gratefully acknowledge the
Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this
project by providing computing time on the GCS Supercomputer SuperMUC-NG at
Leibniz Supercomputing Centre (www.lrz.de) under application ID 21761 (PI:
Alvarado-Gómez). JJC and KP acknowledge funding from the German Leibniz
Community under project number P67/2018. CG was supported by NASA contract
NAS8-03060 to the Chandra X-ray Center. This research has made use of NASA’s
Astrophysics Data System Bibliographic Services.

##  5 Data availability

The data would be made available to the community on reasonable request due to
the volume of the 3D simulations. Extractions of specific quantities discussed
in the paper could be requested from the corresponding author.

## References

  * Ahuir et al. (2020) Ahuir J., Brun A. S., Strugarek A., 2020, A&A, 635, A170
  * Airapetian et al. (2017) Airapetian V. S., Glocer A., Khazanov G. V., Loyd R. O. P., France K., Sojka J., Danchi W. C., Liemohn M. W., 2017, ApJ, 836, L3
  * Alvarado-Gómez et al. (2015) Alvarado-Gómez J. D., et al., 2015, A&A, 582, A38
  * Alvarado-Gómez et al. (2016) Alvarado-Gómez J. D., Hussain G. A. J., Cohen O., Drake J. J., Garraffo C., Grunhut J., Gombosi T. I., 2016, A&A, 594, A95
  * Alvarado-Gómez et al. (2018) Alvarado-Gómez J. D., Drake J. J., Cohen O., Moschou S. P., Garraffo C., 2018, ApJ, 862, 93
  * Alvarado-Gómez et al. (2019) Alvarado-Gómez J. D., Drake J. J., Moschou S. P., Garraffo C., Cohen O., NASA LWS Focus Science Team: Solar-Stellar Connection Yadav R. K., Fraschetti F., 2019, ApJ, 884, L13
  * Alvarado-Gómez et al. (2020) Alvarado-Gómez J. D., Drake J. J., Garraffo C., Cohen O., Poppenhaeger K., Yadav R. K., Moschou S. P., 2020, ApJ, 902, L9
  * Alvarado-Gómez et al. (2022a) Alvarado-Gómez J. D., Drake J. J., Cohen O., Fraschetti F., Garraffo C., Poppenhäger K., 2022a, Astronomische Nachrichten, 343, e10100
  * Alvarado-Gómez et al. (2022b) Alvarado-Gómez J. D., et al., 2022b, ApJ, 928, 147
  * Aschwanden (2006) Aschwanden M., 2006, Physics of the Solar Corona: An Introduction with Problems and Solutions. Springer Praxis Books, Springer Berlin Heidelberg, https://books.google.de/books?id=W7FE5_aowEQC
  * Bicz et al. (2022) Bicz K., Falewicz R., Pietras M., Siarkowski M., Preś P., 2022, ApJ, 935, 102
  * Biermann (1957) Biermann L., 1957, The Observatory, 77, 109
  * Boro Saikia et al. (2020a) Boro Saikia S., Jin M., Johnstone C., Lüftinger T., Güdel M., Airapetian V., Kislyakova K., Folsom C., 2020a, Astronomy & Astrophysics, 635
  * Boro Saikia et al. (2020b) Boro Saikia S., Jin M., Johnstone C. P., Lüftinger T., Güdel M., Airapetian V. S., Kislyakova K. G., Folsom C. P., 2020b, A&A, 635, A178
  * Charbonneau (2020) Charbonneau P., 2020, Living Reviews in Solar Physics, 17, 4
  * Chebly et al. (2021) Chebly J. J., Alvarado-Gómez J. D., Poppenhaeger K., 2021, arXiv e-prints, p. arXiv:2111.09707
  * Cohen (2011) Cohen O., 2011, MNRAS, 417, 2592
  * Cohen et al. (2010) Cohen O., Drake J. J., Kashyap V. L., Sokolov I. V., Gombosi T. I., 2010, ApJ, 723, L64
  * Cohen et al. (2014) Cohen O., Drake J. J., Glocer A., Garraffo C., Poppenhaeger K., Bell J. M., Ridley A. J., Gombosi T. I., 2014, ApJ, 790, 57
  * Cohen et al. (2020) Cohen O., Garraffo C., Moschou S.-P., Drake J. J., Alvarado-Gómez J. D., Glocer A., Fraschetti F., 2020, ApJ, 897, 101
  * Cranmer (2012) Cranmer S. R., 2012, Space Sci. Rev., 172, 145
  * Cranmer et al. (2015) Cranmer S. R., Asgari-Targhi M., Miralles M. P., Raymond J. C., Strachan L., Tian H., Woolsey L. N., 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140148
  * Cristofari et al. (2023) Cristofari P. I., et al., 2023, MNRAS, 522, 1342
  * Donati et al. (2006) Donati J.-F., et al., 2006, Monthly Notices of the Royal Astronomical Society, 370, 629
  * Drake et al. (1993) Drake S. A., Simon T., Brown A., 1993, ApJ, 406, 247
  * Drake et al. (2013) Drake J. J., Cohen O., Yashiro S., Gopalswamy N., 2013, ApJ, 764, 170
  * Egan et al. (2019a) Egan H., Jarvinen R., Brain D., 2019a, MNRAS, 486, 1283
  * Egan et al. (2019b) Egan H., Jarvinen R., Ma Y., Brain D., 2019b, Monthly Notices of the Royal Astronomical Society, 488, 2108
  * Evensberget et al. (2022) Evensberget D., Carter B. D., Marsden S. C., Brookshaw L., Folsom C. P., Salmeron R., 2022, MNRAS, 510, 5226
  * Fares et al. (2009) Fares R., et al., 2009, MNRAS, 398, 1383
  * Feldman (1977) Feldman W. C., 1977, J. Geophys. Res., 82, 667
  * Fichtinger et al. (2017) Fichtinger B., Gúdel M., Mutel R. L., Hallinan G., Gaidos E., Skinner S. L., Lynch C., Gayley K. G., 2017, A&A, 599, A127
  * Finley et al. (2018) Finley A. J., Matt S. P., See V., 2018, ApJ, 864, 125
  * Finley et al. (2019) Finley A. J., See V., Matt S. P., 2019, ApJ, 876, 44
  * France et al. (2020) France K., et al., 2020, AJ, 160, 237
  * Gaidos et al. (2000) Gaidos E. J., Güdel M., Blake G. A., 2000, Geophys. Res. Lett., 27, 501
  * Gallet & Bouvier (2013) Gallet F., Bouvier J., 2013, A&A, 556, A36
  * Gallet & Bouvier (2015) Gallet F., Bouvier J., 2015, A&A, 577, A98
  * Garraffo et al. (2015a) Garraffo C., Drake J. J., Cohen O., 2015a, ApJ, 807, L6
  * Garraffo et al. (2015b) Garraffo C., Drake J. J., Cohen O., 2015b, ApJ, 813, 40
  * Garraffo et al. (2016) Garraffo C., Drake J. J., Cohen O., 2016, ApJ, 833, L4
  * Garraffo et al. (2017) Garraffo C., Drake J. J., Cohen O., Alvarado-Gómez J. D., Moschou S. P., 2017, ApJ, 843, L33
  * Garraffo et al. (2018) Garraffo C., et al., 2018, ApJ, 862, 90
  * Garraffo et al. (2022) Garraffo C., Alvarado-Gómez J. D., Cohen O., Drake J. J., 2022, ApJ, 941, L8
  * Gastine et al. (2013) Gastine T., Morin J., Duarte L., Reiners A., Christensen U. R., Wicht J., 2013, A&A, 549, L5
  * Gombosi (2004) Gombosi T. I., 2004, Physics of the Space Environment. Cambridge University Press
  * Gombosi et al. (2018) Gombosi T. I., van der Holst B., Manchester W. B., Sokolov I. V., 2018, Living Reviews in Solar Physics, 15, 4
  * Gombosi et al. (2021) Gombosi T. I., et al., 2021, Journal of Space Weather and Space Climate, 11, 42
  * Gry & Jenkins (2014) Gry C., Jenkins E. B., 2014, A&A, 567, A58
  * Güdel et al. (1997) Güdel M., Guinan E. F., Skinner S. L., 1997, ApJ, 483, 947
  * Gurnett & Kurth (2019) Gurnett D., Kurth W., 2019, Nature Astronomy, 3
  * Hansteen & Velli (2012) Hansteen V. H., Velli M., 2012, Space Sci. Rev., 172, 89
  * Harbach et al. (2021) Harbach L. M., Moschou S. P., Garraffo C., Drake J. J., Alvarado-Gómez J. D., Cohen O., Fraschetti F., 2021, ApJ, 913, 130
  * Hawley et al. (2014) Hawley S. L., Davenport J. R. A., Kowalski A. F., Wisniewski J. P., Hebb L., Deitrick R., Hilton E. J., 2014, ApJ, 797, 121
  * Hussain et al. (2016) Hussain G. A. J., et al., 2016, A&A, 585, A77
  * Jardine & Collier Cameron (2019) Jardine M., Collier Cameron A., 2019, MNRAS, 482, 2853
  * Jeffers et al. (2014) Jeffers S. V., Petit P., Marsden S. C., Morin J., Donati J. F., Folsom C. P., 2014, A&A, 569, A79
  * Jeffers et al. (2017) Jeffers S. V., Boro Saikia S., Barnes J. R., Petit P., Marsden S. C., Jardine M. M., Vidotto A. A., the BCool collaboration 2017, Monthly Notices of the Royal Astronomical Society: Letters, 471, L96
  * Jivani et al. (2023) Jivani A., et al., 2023, Space Weather, 21, e2022SW003262
  * Johnstone et al. (2015a) Johnstone C. P., Güdel M., Lüftinger T., Toth G., Brott I., 2015a, A&A, 577, A27
  * Johnstone et al. (2015b) Johnstone C. P., Güdel M., Brott I., Lüftinger T., 2015b, A&A, 577, A28
  * Jordan & Montesinos (1991) Jordan C., Montesinos B., 1991, Monthly Notices of the Royal Astronomical Society, 252, 21P
  * Kasper et al. (2021) Kasper J. C., et al., 2021, Phys. Rev. Lett., 127, 255101
  * Kasting et al. (1993) Kasting J. F., Whitmire D. P., Reynolds R. T., 1993, Icarus, 101, 108
  * Kavanagh et al. (2021) Kavanagh R. D., Vidotto A. A., Klein B., Jardine M. M., Donati J.-F., Ó Fionnagáin D., 2021, MNRAS, 504, 1511
  * Kislyakova et al. (2014) Kislyakova K. G., et al., 2014, A&A, 562, A116
  * Klimchuk (2015) Klimchuk J. A., 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140256
  * Kochukhov (2020) Kochukhov O., 2020, in Wade G., Alecian E., Bohlender D., Sigut A., eds, Vol. 11, Stellar Magnetism: A Workshop in Honour of the Career and Contributions of John D. Landstreet. pp 89–101 (arXiv:1912.07241), doi:10.48550/arXiv.1912.07241
  * Kopparapu et al. (2014) Kopparapu R. K., Ramirez R. M., SchottelKotte J., Kasting J. F., Domagal-Goldman S., Eymet V., 2014, The Astrophysical Journal, 787, L29
  * Koutroumpa et al. (2009) Koutroumpa D., Lallement R., Raymond J. C., Kharchenko V., 2009, ApJ, 696, 1517
  * Kowalski et al. (2013) Kowalski A. F., Hawley S. L., Wisniewski J. P., Osten R. A., Hilton E. J., Holtzman J. A., Schmidt S. J., Davenport J. R. A., 2013, ApJS, 207, 15
  * Lacy et al. (1976) Lacy C. H., Moffett T. J., Evans D. S., 1976, The Astrophysical Journal Supplement Series, 30, 85
  * Leitzinger & Odert (2022) Leitzinger M., Odert P., 2022, Serbian Astronomical Journal, 205, 1
  * Leto et al. (1997) Leto G., Pagano I., Buemi C., Rodono M., 1997, VizieR Online Data Catalog, pp J–A+
  * Liefke & Schmitt (2007) Liefke C., Schmitt J. H. M. M., 2007, NEXXUS 2 - The Database for Nearby X-Ray and Extreme UV emitting Stars, https://hsweb.hs.uni-hamburg.de/projects/nexxus/index.html
  * Lim & White (1996) Lim J., White S. M., 1996, The Astrophysical Journal Letters, 462, L91
  * Linsky & Wood (2014) Linsky J. L., Wood B. E., 2014, ASTRA Proceedings, 1, 43
  * Loyd et al. (2022) Loyd R. O. P., et al., 2022, ApJ, 936, 170
  * Lugaz et al. (2015) Lugaz N., Farrugia C., Huang C.-L., Spence H., 2015, Geophysical Research Letters, 42
  * Maehara et al. (2021) Maehara H., et al., 2021, PASJ, 73, 44
  * Matt et al. (2012) Matt S. P., MacGregor K. B., Pinsonneault M. H., Greene T. P., 2012, ApJ, 754, L26
  * McComas et al. (1998) McComas D., et al., 1998, Geophysical Research Letters, 25, 1
  * McComas et al. (2003) McComas D. J., Elliott H. A., Schwadron N. A., Gosling J. T., Skoug R. M., Goldstein B. E., 2003, Geophys. Res. Lett., 30, 1517
  * Metcalfe et al. (2013) Metcalfe T. S., et al., 2013, ApJ, 763, L26
  * Mitra-Kraev et al. (2005) Mitra-Kraev U., et al., 2005, Astronomy & Astrophysics, 431, 679
  * Morin et al. (2008) Morin J., et al., 2008, MNRAS, 390, 567
  * Morin et al. (2010) Morin J., Donati J. F., Petit P., Delfosse X., Forveille T., Jardine M. M., 2010, MNRAS, 407, 2269
  * Moschou et al. (2019) Moschou S.-P., Drake J. J., Cohen O., Alvarado-Gómez J. D., Garraffo C., Fraschetti F., 2019, ApJ, 877, 105
  * Muheki et al. (2020) Muheki P., Guenther E. W., Mutabazi T., Jurua E., 2020, MNRAS, 499, 5047
  * Noyes et al. (1984) Noyes R. W., Hartmann L. W., Baliunas S. L., Duncan D. K., Vaughan A. H., 1984, ApJ, 279, 763
  * Odert et al. (2017) Odert P., Leitzinger M., Hanslmeier A., Lammer H., 2017, MNRAS, 472, 876
  * Ofman (2010) Ofman L., 2010, Living Reviews in Solar Physics, 7, 4
  * Oran et al. (2017) Oran R., Landi E., van der Holst B., Sokolov I. V., Gombosi T. I., 2017, ApJ, 845, 98
  * Parker (1958) Parker E. N., 1958, ApJ, 128, 664
  * Pecaut & Mamajek (2013) Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9
  * Pizzolato et al. (2000) Pizzolato N., Maggio A., Sciortino S., 2000, A&A, 361, 614
  * Pizzolato et al. (2003) Pizzolato N., Maggio A., Micela G., Sciortino S., Ventura P., 2003, A&A, 397, 147
  * Powell et al. (1999) Powell K. G., Roe P. L., Linde T. J., Gombosi T. I., De Zeeuw D. L., 1999, Journal of Computational Physics, 154, 284
  * Pulkkinen (2007) Pulkkinen T., 2007, Living Reviews in Solar Physics, 4, 1
  * Ramirez (2018) Ramirez R. M., 2018, Geosciences, 8, 280
  * Ramstad & Barabash (2021) Ramstad R., Barabash S., 2021, Space Sci. Rev., 217, 36
  * Redfield & Linsky (2004a) Redfield S., Linsky J. L., 2004a, ApJ, 602, 776
  * Redfield & Linsky (2004b) Redfield S., Linsky J. L., 2004b, ApJ, 613, 1004
  * Redfield & Linsky (2008a) Redfield S., Linsky J. L., 2008a, The Astrophysical Journal, 673, 283
  * Redfield & Linsky (2008b) Redfield S., Linsky J. L., 2008b, ApJ, 673, 283
  * Redfield & Linsky (2015) Redfield S., Linsky J. L., 2015, ApJ, 812, 125
  * Reiners et al. (2022) Reiners A., et al., 2022, A&A, 662, A41
  * Ribas et al. (2005) Ribas I., Guinan E. F., Güdel M., Audard M., 2005, ApJ, 622, 680
  * Riley et al. (2014) Riley P., et al., 2014, Sol. Phys., 289, 769
  * Sachdeva et al. (2019) Sachdeva N., et al., 2019, ApJ, 887, 83
  * Schmidt et al. (1995) Schmidt K. H., Bögner M., Buschhorn G., Kotthaus R., Oberhuber R., Rzepka M., 1995, Nuclear Instruments and Methods in Physics Research A, 367, 215
  * See et al. (2019a) See V., et al., 2019a, ApJ, 876, 118
  * See et al. (2019b) See V., et al., 2019b, ApJ, 876, 118
  * See et al. (2019c) See V., et al., 2019c, The Astrophysical Journal, 886, 120
  * Shields et al. (2016) Shields A. L., Ballard S., Johnson J. A., 2016, Phys. Rep., 663, 1
  * Shulyak et al. (2010) Shulyak D., Reiners A., Wende S., Kochukhov O., Piskunov N., Seifahrt A., 2010, A&A, 523, A37
  * Skumanich (1972) Skumanich A., 1972, ApJ, 171, 565
  * Sokolov et al. (2013) Sokolov I. V., et al., 2013, ApJ, 764, 23
  * Stelzer et al. (2016) Stelzer B., Damasso M., Scholz A., Matt S. P., 2016, Monthly Notices of the Royal Astronomical Society, 463, 1844
  * Stone et al. (1998) Stone E. C., Frandsen A. M., Mewaldt R. A., Christian E. R., Margolies D., Ormes J. F., Snow F., 1998, Space Sci. Rev., 86, 1
  * Strugarek (2018) Strugarek A., 2018, in Deeg H. J., Belmonte J. A., eds, , Handbook of Exoplanets. p. 25, doi:10.1007/978-3-319-55333-7_25
  * Strugarek (2021) Strugarek A., 2021, arXiv e-prints, p. arXiv:2104.05968
  * Suzuki (2006) Suzuki T. K., 2006, ApJ, 640, L75
  * Tian (2008) Tian F., 2008, in AAS/Division for Planetary Sciences Meeting Abstracts #40. p. 14.06
  * Tóth et al. (2011) Tóth G., van der Holst B., Huang Z., 2011, ApJ, 732, 102
  * Tóth et al. (2012) Tóth G., et al., 2012, Journal of Computational Physics, 231, 870
  * Tóth et al. (2005) Tóth G., et al., 2005, Journal of Geophysical Research, 110, A12226
  * Van Doorsselaere et al. (2020) Van Doorsselaere T., et al., 2020, Space Sci. Rev., 216, 140
  * Velli et al. (2015) Velli M., Pucci F., Rappazzo F., Tenerani A., 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140262
  * Vidotto (2021) Vidotto A. A., 2021, Living Reviews in Solar Physics, 18, 3
  * Vidotto & Bourrier (2017) Vidotto A. A., Bourrier V., 2017, MNRAS, 470, 4026
  * Vidotto et al. (2013) Vidotto A. A., Jardine M., Morin J., Donati J. F., Lang P., Russell A. J. B., 2013, A&A, 557, A67
  * Vidotto et al. (2014a) Vidotto A. A., Jardine M., Morin J., Donati J. F., Opher M., Gombosi T. I., 2014a, MNRAS, 438, 1162
  * Vidotto et al. (2014b) Vidotto A. A., et al., 2014b, Monthly Notices of the Royal Astronomical Society, 441, 2361
  * Vidotto et al. (2015) Vidotto A. A., Fares R., Jardine M., Moutou C., Donati J. F., 2015, MNRAS, 449, 4117
  * Vidotto et al. (2016) Vidotto A. A., et al., 2016, MNRAS, 455, L52
  * Villadsen et al. (2014) Villadsen J., Hallinan G., Bourke S., Güdel M., Rupen M., 2014, ApJ, 788, 112
  * Villarreal D’Angelo et al. (2017) Villarreal D’Angelo C., Jardine M., See V., 2017, Monthly Notices of the Royal Astronomical Society: Letters, 475, L25
  * Wargelin & Drake (2002) Wargelin B. J., Drake J. J., 2002, ApJ, 578, 503
  * Weber & Davis (1967) Weber E. J., Davis Leverett J., 1967, ApJ, 148, 217
  * Wood (2004) Wood B. E., 2004, Living Reviews in Solar Physics, 1, 2
  * Wood et al. (2002) Wood B. E., Müller H.-R., Zank G. P., Linsky J. L., 2002, ApJ, 574, 412
  * Wood et al. (2005a) Wood B. E., Redfield S., Linsky J. L., Müller H.-R., Zank G. P., 2005a, The Astrophysical Journal Supplement Series, 159, 118
  * Wood et al. (2005b) Wood B. E., Müller H. R., Zank G. P., Linsky J. L., Redfield S., 2005b, ApJ, 628, L143
  * Wood et al. (2014) Wood B. E., Müller H.-R., Redfield S., Edelman E., 2014, ApJ, 781, L33
  * Wood et al. (2021) Wood B. E., et al., 2021, ApJ, 915, 37
  * Wright & Drake (2016) Wright N. J., Drake J. J., 2016, Nature, 535, 526
  * Wright et al. (2011a) Wright N. J., Drake J. J., Mamajek E. E., Henry G. W., 2011a, ApJ, 743, 48
  * Wright et al. (2011b) Wright N. J., Drake J. J., Mamajek E. E., Henry G. W., 2011b, ApJ, 743, 48
  * Wright et al. (2018a) Wright N. J., Newton E. R., Williams P. K. G., Drake J. J., Yadav R. K., 2018a, MNRAS, 479, 2351
  * Wright et al. (2018b) Wright N. J., Newton E. R., Williams P. K. G., Drake J. J., Yadav R. K., 2018b, Monthly Notices of the Royal Astronomical Society, 479, 2351
  * van Saders et al. (2016) van Saders J. L., Ceillier T., Metcalfe T. S., Silva Aguirre V., Pinsonneault M. H., García R. A., Mathur S., Davies G. R., 2016, Nature, 529, 181
  * van Saders et al. (2019) van Saders J. L., Pinsonneault M. H., Barbieri M., 2019, ApJ, 872, 128
  * van den Oord & Doyle (1997) van den Oord G. H. J., Doyle J. G., 1997, A&A, 319, 578
  * van der Holst et al. (2014) van der Holst B., Sokolov I. V., Meng X., Jin M., Manchester W. B. I., Tóth G., Gombosi T. I., 2014, ApJ, 782, 81
  * van der Holst et al. (2019) van der Holst B., Manchester W. B. I., Klein K. G., Kasper J. C., 2019, ApJ, 872, L18

##  Appendix A Trends with Maximum Radial Magnetic Field

We have also quantified $\overline{AS}_{\rm R}$ and $\dot{M}_{\rm\bigstar}$ /
$R^{\rm 2}_{\bigstar}$ as a function of the absolute maximum radial magnetic
field strength ($|B_{\rm R}|^{\rm max}|$). It is important to also investigate
$|B_{\rm R}|^{\rm max}$, since the average radial magnetic field strength may
suffer from cancellations, especially if the star has a symmetric surface
magnetic field distribution. Figure 10 shows the simulated average Alfvén
surface area ($\overline{AS}$, top) and the mass-loss rate per unit surface
area ($\dot{M}_{\bigstar}/R^{\rm 2}_{\rm\bigstar}$, bottom) as a function of
the maximum absolute radial magnetic field on the stellar surface ($|B_{\rm
R}|^{\rm max}$). We see a trend where $\overline{AS}$ and
$\dot{M}_{\bigstar}/R^{\rm 2}_{\rm\bigstar}$ increase with increasing magnetic
field strength. We fit this trend to a power law by applying the bootstrap
method used to derive this parameter as a function of the average radial
magnetic field (similar to the procedure used in Sect. 3.1, Eqs. 3 and 4).

Figure 10: Simulated average Alfvén surface ($\overline{AS}$, top), and the
mass loss rate per unit surface area ($\dot{M}_{\bigstar}/R^{\rm
2}_{\rm\bigstar}$, bottom) as a function of the absolute surface-max radial
magnetic field ($|B_{\rm R}|^{\rm max}$). The mass loss rate is expressed in
units of the average solar values, normalized to the surface area of each
star. Individual points denote the results of each simulation presented in
Sect. 3, Table 2. The different symbols and colors represent the spectral
types (F, cyan/squares; G, yellow/diamonds; K, orange/circles; M, red/star).
The purple line with the shaded purple area represents the fitted power-law
with its uncertainties.

$\log\overline{AS_{\rm R}}=(0.44\pm 0.05)\log|B_{\rm R}|^{\rm max}+(0.54\pm
0.08)$ (9)

$\log\dot{M}_{\rm\bigstar}/R^{\rm 2}_{\bigstar}=(0.83\pm 0.07)\log|B_{\rm
R}|^{\rm max}-(0.48\pm 0.10)$ (10)