The winds of young Solar-type stars in the Pleiades, AB Doradus, Columba, and β Pictoris Open Access

ABSTRACT

Solar-type stars, which shed angular momentum via magnetized stellar winds, enter the main sequence with a wide range of rotational periods P_rot. This initially wide range of rotational periods contracts and has mostly vanished by a stellar age |$t\sim {0.6}\, {\rm Gyr}$|⁠, after which Solar-type stars spin according to the Skumanich relation |$P_\text{rot}\propto \sqrt{t}$|⁠. Magnetohydrodynamic stellar wind models can improve our understanding of this convergence of rotation periods. We present wind models of 15 young Solar-type stars aged ∼24 Myr to ∼0.13 Gyr. With our previous wind models of stars aged ∼0.26 and ∼0.6 Gyr we obtain 30 consistent three-dimensional wind models of stars mapped with Zeeman–Doppler imaging – the largest such set to date. The models provide good cover of the pre-Skumanich phase of stellar spin-down in terms of rotation, magnetic field, and age. We find the mass-loss rate |$\dot{M}\propto \Phi ^{{0.9\pm 0.1}}$| with a residual spread of ∼150 per cent and the wind angular momentum loss rate |$\dot{J}\propto {}P_\text{rot}^{-1} \Phi ^{1.3\pm 0.2}$| with a residual spread of ∼500 per cent where Φ is the unsigned surface magnetic flux. When comparing different magnetic field scalings for each single star we find a gradual reduction in the power-law exponent with increasing magnetic field strength.

1 INTRODUCTION

The pre-Skumanich spin-down phase of a cool star’s life is the era between the end of the contractive spin-up phase and the Skumanich (1972) spin-down phase where the stellar period of rotation |$P_\text{rot}\propto \sqrt{t}$| , where t is the stellar age. For Solar-type stars the pre-Skumanich spin-down phase ends around 0.6 Gyr (Gallet & Bouvier 2013, 2015). As stars enter the pre-Skumanich spin-down phase with a spread of rotational periods of at least an order of magnitude (Edwards et al. 1993), angular momentum shedding in the pre-Skumanich spin-down phase must be instrumental in permitting the convergence of rotational periods with increasing stellar age.

Angular momentum shedding occurs by means of magnetized stellar winds (Schatzman 1962). Semi-empirical models of angular momentum loss on the general form |$\dot{J} \propto \dot{M} P_\text{rot}^{-1} R_{\small {\rm A}}^n$| have been proposed in many studies (Weber & Davis 1967; Mestel 1968, 1984; Kawaler 1988). In these models |$\dot{M}$| is the stellar winds mass-loss rate and |$R_{\small {\rm A}}$| denotes the average radius of the Alfvén surface where the wind’s kinetic energy and magnetic energy have similar magnitudes. The average Alfvén radius |$R_{\small {\rm A}}$| depends not only on the large-scale coronal magnetic field strength |B| but also on |$\dot{M}$| as both the wind kinetic energy and mass-loss are determined by the wind speed and density. The wind mass-loss rate |$\dot{M}$| is poorly constrained by observations, particularly for stars younger than 0.6 Gyr (Wood et al. 2005, 2014); this uncertainty propagates to our knowledge of |$\dot{J}$| and |$R_{\small {\rm A}}$|⁠. The exponent n ≲ 2 is a geometric parameter which decreases with increasing complexity of the large-scale coronal magnetic field.

The magnitude of the large-scale photospheric magnetic field decreases with increasing P_rot (Marsden et al. 2014), as does stellar magnetic activity in general (Noyes et al. 1984). This decrease is presumably caused by a throttling of the internal stellar dynamo as the amount of available rotational energy decreases. Qualitative dynamo changes with increasing P_rot (see Barnes 2003; Jeffers et al. 2011; Morin et al. 2011; Marsden et al. 2011a, b; Brown 2014) would also affect the complexity of the photospheric magnetic field as it emerges from lower stellar layers.

Both the strength and complexity of the large-scale stellar photospheric magnetic field can be studied with spectropolarimetric stellar observations and Zeeman–Doppler imaging (ZDI, Semel 1989; Donati et al. 1997). While the average large-scale magnetic field intensity found using ZDI is subject to uncertainty (Lehmann et al. 2019), on average the ZDI-derived field decays as |$|B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}|\propto P_\text{rot}^{{-1.32\pm 0.14}}$|⁠; however, the relation is subject to significant scatter (Vidotto et al. 2014b). ZDI tends to yield more complex magnetic field geometries in younger stars (Folsom et al. 2018), but this may be due to the technique’s increased resolving power for rapidly spinning stars (Morin et al. 2010).

The early three-dimensional numerical wind models by Vidotto et al. (2009a, 2012) and Cohen et al. (2010) show that field complexity can significantly affect the angular momentum loss rate |$\dot{J}$|⁠. This highlights the need for incorporating realistic large-scale magnetic fields when studying the |$|B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}|$|–P_rot relation. Numerical models also suggest that |$\dot{M}$| (and thus also |$\dot{J}$|⁠) may be inhibited by a complex magnetic field (Garraffo, Drake & Cohen 2015).

In this work, we investigate the effect of the photospheric magnetic field and stellar rotation rate on wind mass- and angular momentum loss rates as well as other parameters of interest by driving the state-of-the-art numerical magnetohydrodynamic Alfvén wave solar model (awsom, Sokolov et al. 2013; van der Holst et al. 2014) with the photospheric radial magnetic field from the magnetic maps of young, Solar-type stars published by Folsom et al. (2016, 2018, hereafter F16, F18). We consider 11 stars with well-constrained ages in the ∼130 Myr old Pleaides and AB Doradus clusters and four very young stars in the ∼42 Myr old Columba association and ∼24 Myr old β Pictoris association.

By letting the resulting three-dimensional stellar wind models relax towards a steady state, we obtain simultaneous simulated values for the wind mass-loss rate |$\dot{M}$|⁠, wind angular momentum loss rate |$\dot{J}$| and the Alfvén radius |$R_{\small {\rm A}}$|⁠, wind pressure for an Earth-like planet,¹ and many other quantities of interest. This permits us to investigate the degree to which the coronal magnetic field complexity affects stellar spin-down.

For Solar-type stars the pre-Skumanich spin-down phase ends around 0.6 Gyr (Gallet & Bouvier 2013, 2015). By combining the results of this work with our previous wind models of ∼0.6 Gyr old Solar-type stars in the Hyades (Evensberget et al. 2021, hereafter Paper I) and ∼0.26 Gyr old stars in Coma Berenices and Hercules-Lyra (Evensberget et al. 2022, hereafter Paper II), we get good age coverage of the pre-Skumanich spin-down period and a large sample of 30 stellar models from which trends may be robustly extracted.

This paper is outlined as follows: In Section 2, we briefly describe the spectropolarimetric observations upon which the magnetic maps used in this work are based, as well as the process by which the surface magnetic field maps are obtained. In Section 3, we describe the stellar wind model in terms of physical effects included, model equations, three-dimensional numerical model, and boundary conditions. In Section 4, we present visualizations of the wind model coronal magnetic field, Alfvén surface and current sheet, and wind pressure out to 1 au, and tabulate aggregate wind parameters calculated from the three-dimensional wind maps. In Section 5, we investigate the effects of scaling the magnetic field, the statistical trends, and variation of key wind parameters when plotted against the magnetic field strength and related parameters. In Section 6, we conclude our work and put our results into context with other recent works in the field.

2 OBSERVATIONS

The stellar surface magnetic maps used in this work were derived from ZDI (Semel 1989) observations and modelling by F16, F18. As part of ‘TOwards Understanding the sPIn Evolution of Stars’ (TOUPIES)² project, spectropolarimetric observations in circularly polarized and unpolarized light of the stellar targets in Table 1 were carried out using the ESPaDOnS instrument (Donati 2003; Silvester et al. 2012) at the Canada–Hawaii–France Telescope (CHFT) and using the Narval instrument (Aurière 2003) at the Télescope Bernard Lyot (TBL) between 2009 and 2015. The ESPaDOnS observations of BD-07 2388 and HD 6569 were part of the History of the Magnetic Sun³ Large Program at the CFHT. The TOUPIES reduced spectra can be downloaded from the Polarbase (Donati et al. 1997; Petit et al. 2014) website.⁴ Each star was observed over a period of a few weeks to obtain good phase coverage of the star while minimizing intrinsic stellar variation. The resulting surface magnetic field maps are published in (F16, F18);⁵ the radial component of the maps, which are used to drive the wind models of this work, are reproduced in Fig. 1.

2.1 Magnetic mapping with ZDI

This section gives a short overview of the ZDI technique on which the stellar magnetic maps are based. ZDI (Semel 1989) is a tomographic technique that permits the reconstruction of the stellar surface magnetic field. The technique requires high-resolution spectropolarimetric observations from multiple epochs in order to reconstruct a two-dimensional image of the vector surface magnetic field. In ZDI, least-squares deconvolution (LSD, Donati & Brown 1997; Kochukhov, Makaganiuk & Piskunov 2010) is used to combine magnetically sensitive spectral line profiles in velocity space into a single profile with a much higher signal-to-noise ratio than individual lines.

By comparing the resulting circularly polarized and unpolarized LSD line profiles, and by combining multiple LSD profiles from different time and stellar phases, it is possible to reconstruct a two-dimensional surface map of the large-scale magnetic features on the stellar surface in a robust way (Hussain et al. 2000). Linearly polarized light may also be included and has the potential to break degeneracies, but as the linearly polarized signal is ∼10 times weaker than the circularly polarized signal, this is only feasible for a subset of the stars to which ZDI has been applied, see Rosén, Kochukhov & Wade (2015). The review by Donati & Landstreet (2009) gives an overview of ZDI applied to cool stars.

The amount of detail in the magnetic map depends on unpolarized line width, projected rotational velocity vsin i, signal-to-noise ratio of the observations, and (implementation dependent) choices of fitting parameters such as the χ² target values used in F16, F18, see the discussion in Donati & Brown (1997) and Morin et al. (2010).

Modern ZDI gives a parametric representation of the vector surface magnetic field |$\boldsymbol{B}(\theta ,\varphi)$| in terms of spherical harmonics coefficients (Jardine et al. 1999; Donati et al. 2006). Here, θ and φ denote the polar and azimuth position on the stellar surface in a co-rotating spherical coordinate system (similar to colatitude and longitude on Earth). For the magnetograms in this study, the surface radial magnetic field B_r(θ, φ) is parametrized by a set of complex-valued harmonic coefficients α_ℓm and calculated via a spherical harmonics expansion as

where |$\operatorname{Re}$| denotes taking the real part, P_ℓm(μ) is the associated Legendre polynomial of order m and degree ℓ, and μ = cos θ. Since we discard the imaginary component of the right-hand side in equation (1) the non-negative orders m ≥ 0 provide sufficient degrees of freedom to represent any magnetic field configuration. The maximum degree ℓ_max is set to 15, giving a minimum angular resolution of 180°/ℓ_max = 12° in the magnetograms.

When comparing magnetograms with similar ℓ_max values, it is important to note that ℓ_max is only an upper bound on the representable magnetogram complexity. For slow rotators there may not be sufficient signal in the ZDI profile to reconstruct features past the first few degrees, as the signal scales with vsin i (Morin et al. 2010) and thus decreases with increasing P_rot. To quantify the physical complexity of the magnetograms we calculate an ‘effective degree’ ℓ_.90 and ℓ_.99 (see Paper II) for which, respectively, 90 per cent and 99 per cent of the wind model magnetic energy is contained in degrees smaller than or equal to ℓ. The effective degree of the magnetograms in Fig. 1 is discussed in Section 4.1.

While there is general agreement that ZDI is able to reconstruct the large-scale polarity structure of the field (e.g. Hussain et al. 2000), uncertainty remains surrounding the field strength (Lehmann et al. 2019). When comparing with the complementary Zeeman broadening technique, which is sensitive to both the large- and small-scale field, Yadav et al. (2015) found only about 20 per cent of field to be reconstructed with ZDI. In this work, we control for this possible underreporting of the large-scale magnetic field strength by creating two series of wind models, one of which has its surface magnetic field increased by a factor of 5; we return to this in Section 3 and in the analysis in Section 5.1.

3 SIMULATIONS

We obtain the wind models of this work by numerically solving a two-temperature extended ideal MHD model with physical heating and cooling terms, and letting the models evolve towards a steady-state solution. For this we use the Alfvén wave Solar model (awsom, Sokolov et al. 2013; van der Holst et al. 2014) driven by the surface magnetic fields in Fig. 1. In the awsom model the mechanism of coronal heating is Alfvén wave dissipation (see e.g. Barnes 1968). The awsom model was created to model the Solar wind, in which context it has been the subject of extensive validation (Meng et al. 2015; Sachdeva et al. 2019; van der Holst et al. 2019). The model is also used in stellar winds modelling (e.g. Alvarado-Gómez et al. 2016a, b; Cohen 2017; Garraffo et al. 2017; Kavanagh et al. 2021, and Paper I, Paper II). The awsom model is part of the block-adaptive tree solarwind Roe upwind scheme (bats-r-us, Powell et al. 1999; Tóth et al. 2012) code, which in turn is part of the space weather modelling framework (swmf, Tóth et al. 2005, 2012; Gombosi et al. 2021). We refer the reader to Paper I and Paper II for a detailed description of the model equations and modelling parameters used in this work; here, we provide a brief summary of the coronal heating mechanism in the awsom model and the magnetic field scaling applied to compensate for the issues considered in Section 2.1.

As the awsom module extends outwards from the stellar chromosphere, an Alfvén wave energy flux is prescribed at the inner boundary. The Alfvén wave energy flux is modelled as a boundary Poynting flux-to-field ratio |$\Pi _{\small {\rm A}}/ B$|⁠, such that the local amount of Alfvén wave energy crossing the inner model boundary is proportional to the local value of the magnetic field |$|\boldsymbol{B}(\theta , \varphi)|$|⁠. We note that the parameters used in this paper and Paper I and Paper II have been chosen to match Solar values and reproduce Solar conditions based on Solar magnetograms. As such, they are likely to be more accurate the closer the age of the modelled star is to the age of the Sun, and the values may be more questionable for the |${\le } 0.13\, {\rm Gyr}$| young stars of this study compared to the older stars in Paper I and Paper II. We also consider this issue in Section 6.

In our wind models the numerical domain extends from 1 to 400 stellar radii, split into an inner domain using a spherical grid and an outer domain (beyond 45 stellar radii) using a Cartesian grid. We first find a steady state on the inner grid, then we use this solution to initialize the outer grid. For the outer grid we apply a local artificial wind flux scheme (Sokolov et al. 2002) with an mc3 slope limiter (Koren 1993) as we did in Paper I and Paper II.

For the inner, spherical domain we now apply an HLLE scheme as modified by Linde (2002) and a minmod-type slope limiter (Roe 1986). To incorporate electron heat conduction (see Paper I) we apply a semi-implicit time-stepping method, which effectively treats the thermal conduction operator implicitly (Toth, Keppens & Botchev 1998; Keppens et al. 1999). Furthermore, we have extended the inner domain with an outer buffer region in order to prevent unphysical inflow situations from destabilizing the solution before a steady state is reached. We apply ‘outflow’ boundary conditions at the outer edge of the buffer region with an outflow pressure of 10⁻¹² Pa. This minimizes unphysical inwards-propagating disturbances caused by sub-Alfvénic (see Section 4.2) flows across the outer boundary (these can form in the early stages of a model run but are not present in the final solution). We emphasize that the final, steady-state solutions do not exhibit destabilizing sub-Alfvénic flow at outer boundary in any of our models.

In this work, all the wind models have been created with the open-source version of the SWMF.⁶ The options discussed here are all available via configuration parameters. The application of the above changes to our wind modelling approach has enabled us to model the winds of the young, rapidly rotating stars in Table 1. We have also re-processed our previous wind models (from Paper I, Paper II) with this new and stable numerical configuration, so that the same methodology and model configuration is applied across the whole age range from 24 to 650 Myr. The re-processing of our previous wind models has led to only minor differences, and thus we do not present the re-processed wind models in Section 4; we focus instead on our modelling of the stars of Table 1. In Section 5 onwards we do, however, use the aggregate wind parameters computed from the re-processed wind models and the models in this work, so that we have a fully consistent modelling methodology. The aggregate wind parameters of the re-processed wind models are also given in Appendix  A.

In an attempt to control for ZDI’s inherent uncertainty concerning the magnetic field strength (see Section 2.1), we create two wind models for each star corresponding to different scalings of the surface magnetic field maps of Fig. 1. The |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series uses the unscaled magnetic field of equation (1), while |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series applies a scaling factor of 5 to the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| values (this approach was also used in Paper I and Paper II). In this work, we refer to each wind model by its case name of the form ‘HII 296’ or ‘5|$\, \times \,$|HII 296’; these refer to the star HII 296 from Table 1 with the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| surface magnetic fields, respectively. The symbol shape indicates whether the case belongs to |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| or the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series, and its colours indicate the star’s association (see Table 1).

4 RESULTS

In this section, we present the wind model results and aggregate quantities that may be calculated from the wind model output. We present plots and key quantities describing the final, steady-state coronal magnetic field, the wind speed and Alfvén surface location, the wind mass-loss rate and angular momentum loss rate, and plots of the wind pressure in the equatorial plane extending out to planetary distances.

4.1 Coronal magnetic field

In this section, we describe the final, steady-state coronal magnetic field in our wind model solutions where the hydrodynamic, magnetic, and thermal forces are in balance inside our model domain.

Table 2 gives aggregate magnetic quantities calculated from the steady-state wind models. These quantities are calculated at the inner boundary of the wind model. The quantity |B_r| is the average value of the unsigned radial field over the stellar surface, whereas max |B_r| is similarly the maximum value of the radial field over the stellar surface. The quantities |B_r| and max |B_r| can also be calculated from equation (1) and the two methods are in agreement.

Continuing with the quantities of Table 2, the quantity |$|\boldsymbol{B}|$| is the final, relaxed average (vector) magnetic field strength over the stellar surface in the wind models. The vector field |$\boldsymbol{B} = B_r \boldsymbol{\hat{r}} +\boldsymbol{B}_\perp$| will relax towards a steady state, and the non-radial magnetic field |$\boldsymbol{B}_\perp$| at the stellar surface is free to settle since only the radial field B_r is fixed by a boundary condition. As the bats-r-us awsom code is driven with the radial magnetic field B_r only, |$\boldsymbol{B}_\perp$| is determined by the physics above the chromosphere. This is in contrast with the ZDI magnetic field, where the non-radial magnetic field is also affected by photospheric currents. As the photosphere lies below the chromosphere it is not part of the awsom model; we instead observe that the relaxed values of |$\boldsymbol{B}$| are close to the potential component of the ZDI magnetic field. The potential component of the ZDI magnetic field is thought to determine the magnetic geometry of the corona, while the non-potential field may be stirred by rapid stellar rotation. The ZDI field component originating from photospheric currents is not thought to affect the wind to a significant degree (Jardine et al. 2013).

The quantity Φ in Table 2 is the unsigned surface magnetic flux. The unsigned magnetic flux through a closed surface is

where B is the magnetic field and |$\boldsymbol{\hat{n}}$| is the normal of the surface S. The surface magnetic flux |$\Phi (S_{\mathchoice{}{}{\scriptscriptstyle }{}\bigstar })$| at the inner boundary |$S_{\mathchoice{}{}{\scriptscriptstyle }{}\bigstar }$| of our model is given in Table 2. This quantity can also be found from the equivalent |$\Phi (S_{\mathchoice{}{}{\scriptscriptstyle }{}\bigstar })=4\pi R_{\mathchoice{}{}{\scriptscriptstyle }{}\bigstar }^2 |B_r|$| with excellent agreement between the two methods.

The remaining columns of Table 2 give several measures of the surface magnetic field complexity in each wind model. The ‘Dip.’, ‘Quad.’, ‘Oct.’, and ‘16 + ’ columns give the fraction of magnetic energy contained in magnetogram terms of degree ℓ = 1, ℓ = 2, ℓ = 3, and ℓ ≥ 4, respectively, i.e. the amount of energy in the dipolar, quadrupolar, octupolar, and hexadecapolar and higher degrees. Here, we also see good agreement with the values of F16 and F18. In Section 2.1, we mentioned that the magnetogram ‘effective degree’ may lie significantly below the maximum ℓ_max of the spherical harmonics expansion in equation (1) when most of the magnetic energy is contained in spherical harmonics terms of low degree. In the last two columns of Table 2 we provide two measures of the effective degree, ℓ_.90 and ℓ_.99, which are the magnetogram degrees ℓ for which 90 per cent and 99 per cent of the surface magnetic field energy is contained in degrees less than or equal to ℓ. It may be seen that the rapid rotators BD-07 2388 and LO Peg are the only stars with ℓ_.99 = ℓ_max in our sample. The second highest effective degree is found in PW And with ℓ_.99 = 12, while the rest of the stars in our sample range from ℓ_0.99 = 4 (for HD 6569) to ℓ_.99 = 10. A similar variation is seen for the calculated ℓ_.90 values. These variations in level of reconstructed detail can be recognized in Fig. 1 when comparing the radial magnetic fields of e.g. BD-07 2388 and HD 6569.

Fig. 2 shows the steady-state coronal magnetic field structure for the wind models created in this work. Each panel of the plot shows the coronal magnetic field structure by tracing a large number of magnetic field lines from an evenly sampled set of points on the stellar surface. Due to the presence of hydrodynamic forces in the wind solution, magnetic field lines can be dragged along with the wind velocity field when the hydrodynamic pressure exceeds the magnetic pressure. In our simulations these magnetic field lines terminate at the outer boundary of the wind model; they are termed ‘open’ since they have only one footpoint on the stellar surface. For display purposes the open field lines are truncated at four stellar radii, while closed magnetic field lines, which connect two different points on the stellar surface, are not truncated. The stellar surface and the magnetic field lines are coloured according to the local value of the radial magnetic field strength B_r; we use a colour scale which is linear from −10 to 10 G and logarithmic outside of this range.

Considering a closed spherical surface S of radius R centred on the star, the ‘open magnetic flux’ is the part of the magnetic flux associated with regions of S crossed by open magnetic field lines. The closed magnetic flux is, vice versa, the amount of magnetic flux in regions of S where the magnetic field lines close back on S. Splitting the magnetic flux Φ(S) through S into an open and a closed component Φ(S) = Φ_open(S) + Φ_closed(S) it can be seen that Φ_closed(S) tends towards zero for |$R \gg R_{\mathchoice{}{}{\scriptscriptstyle }{}\bigstar }$|⁠. We use this to calculate open magnetic flux Φ_open by evaluating Φ(S) for |$R \gg R_{\mathchoice{}{}{\scriptscriptstyle }{}\bigstar }$| where Φ_open ≫ Φ_closed. The values we find are consistent with the above argument; we give the open flux values as Φ_open in Table 3.

Table 3 also gives the area of the stellar surface that is crossed by open magnetic field lines in the S_open column. This quantity is found by tracing the evenly spaced magnetic field lines shown in Fig. 2 until the field line either reaches the edge of the computational domain (open field line), or loops back to the stellar surface (closed field line). S_open is then the number of open field lines divided by the total number of field lines. From Table 2 we can see that Φ_open and S_open are well correlated with |$S_\text{open}/S \sim {0.6}\, \Phi _\text{open}/\Phi$|⁠. The complex coronal magnetic fields of HII 739 and PELS 031 (in Fig. 2) do not seem to produce lower values of Φ_open and S_open than the rest of the sample.

In most of our model cases we observe dipole-dominated coronal magnetic fields, except HII 739, PELS 031, and their scaled counterparts, for which the coronal magnetic fields appear more complex. Indeed the amount of dipolar magnetic field energy of these two stars (in Table 2) is the lowest of our sample at 25 per cent and 10 per cent, respectively. The stars BD-07 2388 and LO Peg are notable as they have a very high effective magnetogram degree measure ℓ_0.99 = 15; this does not, however, stop them from forming a dipole-like coronal magnetic field as both stars have ∼40 per cent of their magnetic energy in the dipolar terms of the magnetogram. TYC 0486 also exhibits a dipole-like coronal field with 37.9 per cent of its magnetic energy in dipolar magnetogram terms. This suggests that the lower threshold of dipolar magnetic energy for forming a dipolar coronal field lies at around ∼30 per cent of the total magnetic energy in our wind models.

We observe a range of ‘magnetic inclinations’, measured as the angle between the stellar axis of rotation |$\boldsymbol{\hat{\Omega }}\parallel \boldsymbol{\hat{z}}$| and the ‘north and south pole’ of the magnetic dipole in Fig. 2; in order to quantify this we calculate magnetic inclination in the corona by fitting a plane to the surface where B_r = 0. Letting |$\boldsymbol{\hat{n}}_{{\scriptscriptstyle B}}$| be the normal vector of this plane, we compute the magnetic inclination |$i_{{\scriptscriptstyle B}}= \cos ^{-1}(\boldsymbol{\hat{\Omega }}\cdot \boldsymbol{\hat{n}}_{{\scriptscriptstyle B}})$|⁠. The values of |$i_{{\scriptscriptstyle B}}$| are given in Table 3; we see that HII 296, BD-07 2388, HIP 76768, and TYC 5164 have |$i_{{\scriptscriptstyle B}}\lesssim {15}{^{\circ }}$|⁠; these stars also have low magnetic inclinations in Fig. 2, with the exception of HII 296, which does not exhibit a dipole-dominated coronal magnetic field. TYC 0486 and BD-16 351, which have values of |$i_{{\scriptscriptstyle B}}\gtrsim {70}{^{\circ }}$| are seen to be highly magnetically inclined in Fig. 2. High magnetic inclinations and rapid stellar rotation create more complex wind geometries (see Section 4.2).

The final coronal field quantity we calculate is the axisymmetric open magnetic flux, which has been linked to the amount of cosmic rays propagating in the inner Solar system and impacting the Earth  (Wang, Sheeley & Rouillard 2006) following Vidotto et al. (2014a):

Here, |$\boldsymbol{B}_\text{axi}$| is to be understood as the azimuthally averaged value in spherical coordinates |$\boldsymbol{B} = B_r\boldsymbol{\hat{r}} + B_\theta \boldsymbol{\hat{\theta }}+ B_\varphi \boldsymbol{\hat{\varphi }}$|⁠. We obtain the tabulated values of Φ_open by integrating over a large spherical surface S with radius R ≫ R_⋆, similarly to the calculation of Φ_open. We note that there is a close correlation between Φ_axi/Φ_open and |$\text i_{\small B}$|⁠; this is expected as both parameters measure the lack of magnetic symmetry around the star’s rotation axis.

4.2 Alfvén surface, wind velocity, and current sheet

Alfvén waves propagate with a speed |$u_{\small {\rm A}}= B / \sqrt{\mu _0 \rho }$|⁠, where ρ is the wind mass density and μ₀ is the magnetic constant. Alfvén waves can thus not propagate upstream against the flow velocity |$\boldsymbol{u}$| when the flow speed |$|\boldsymbol{u}| \gt u_{\small {\rm A}}$|⁠. As the wind travels away from the stellar surface, the local magnetic field strength and wind density drops, while the wind speed increases. This leads the local wind speed u to eventually exceed |$u_{\small {\rm A}}$|⁠; the set of points where this first occurs forms a closed surface surrounding the star. This closed surface |$S_{\small {\rm A}}$| is called the Alfvén surface; it separates the inner sub-Alfvénic region from the outer super-Alfvénic region of the wind model. Since the wind flow is nearly radial, the sub-Alfvénic region is disconnected from the super-Alfvénic region in the sense that no Alfvén wave disturbance will cross the Alfvén surface inwards (in the |$-\boldsymbol{\hat{r}}$| direction).

The inability of Alfvén waves to propagate inwards suggests that co-rotating wind cannot effectively drain the star’s supply of angular momentum once outside |$S_{\small {\rm A}}$|⁠. This is the conceptual ‘lever-arm model’ of angular momentum loss proposed by Schatzman (1962). While the conceptual lever-arm model is not accurate, it remains useful as the angular momentum loss is strongly affected by the distance to the Alfvén surface; the average Alfvén radius |$R_{\small {\rm A}}$| is an important parameter in models of stellar spin-down (Weber & Davis 1967; Mestel 1968, 1984; Kawaler 1988) which can be written in the form |$\dot{J} \propto \dot{M} P^{-1}_\text{rot}R_{\small {\rm A}}^n$|⁠, where n ≲ 2 depends on the large-scale magnetic structure of the corona, decreasing with increasing field complexity.

The circumstellar current sheet (Schatten 1971), here characterized by B_r = 0, is a thin sheet-like region of magnetic reconnection and high currents as magnetic field lines of opposite polarity are forced close together. The condition B_r = 0 was also used (in Section 4.1) to calculate the star’s magnetic inclination |$i_{{\scriptscriptstyle B}}$|⁠.

Fig. 3 shows the Alfvén surface as an opaque surface and the inner current sheet as a translucent grey surface. The Alfvén surface and the xz plane is coloured by the local radial wind speed. We observe two-lobed Alfvén surfaces for every star, including the ones without clear dipolar coronal fields. This shows that in our models only largest scale magnetic features become entrained in the wind. The wind radial velocity and the size of the Alfvén surface lobes appear correlated with the surface magnetic field strength |$|\boldsymbol{B}|$|⁠; this is one of many trends we investigate in Section 5. When the stellar rotation rate Ω and the magnetic inclination |$i_{{\scriptscriptstyle B}}$| are high the star’s rotation pulls the current sheet around it in a spiraling structure; this is most evident in HII 739 and PELS 031. When the magnetic field exhibits less symmetry around the |$\boldsymbol{\hat{\Omega }}$| axis the spiral undulations appear more pronounced (see e.g. Vidotto et al. 2010).

For the same star, e.g. HII 739, we see from Fig. 3 that the Alfvén surface tends to lie about a third closer to the star in the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series than in the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series. The average distance to the Alfvén surface is expected to increase as the surface magnetic field B_r increases in value (e.g. Kawaler 1988); we see this increase for all the models although the Alfvén surface of BD-07 2388 extends about as far as that of 5|$\, \times \,$|BD-07 2388 in the +z direction. We calculate |$R_{\small {\rm A}}$| by defining the radial distance to the Alfvén surface |$\boldsymbol{r}_{\small {\rm A}}(\theta , \varphi)$| over each point θ, φ on the stellar surface. The average Alfvén radius is then given by

From the resulting values, which are tabulated in Table 3, it is clear that |$R_{\small {\rm A}}$| increases between the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series.

We also calculate the cylindrical Alfvén radius; this is the mean length of the Alfvén surface’s torquing arm around the |$\boldsymbol{\hat{\Omega }}$| axis. The value is found by crossing |$\boldsymbol{r}_{\small {\rm A}}$| with the axis of rotation, |$\boldsymbol{r}_{\small {\rm A}}(\theta , \varphi)\times \boldsymbol{\hat{\Omega }}$|⁠, which gives the following expression:

The calculated cylindrical Alfvén radii are also given in Table 3. Comparing the growth of |$R_{\small {\rm A}}$| and |$\left|\boldsymbol{r}_{\small {\rm A}}\times \boldsymbol{\hat{\Omega }}\right|$| between BD-07 2388 and 5|$\, \times \,$|BD-07 2388 we see that both quantities have grown by a similar amount, even though the northern Alfvén surface lobe of 5|$\, \times \,$|BD-07 2388 appears to have grown mostly thicker and only a bit taller.

Note the unusually large and conical Alfvén surfaces of the very rapidly (⁠|$P_\text{rot}={0.3}\, {\rm d}$|⁠) rotating BD-07 2388, LO Peg(⁠|$P_\text{rot}={0.4}\, {\rm d}$|⁠), and their scaled counterparts 5|$\, \times \,$|BD-07 2388 and 5|$\, \times \,$|LO Peg. In the 5|$\, \times \,$|BD-07 2388 model the conical southern Alfvén lobe extends past the computational domain so that no value of |$R_{\small {\rm A}}$| and |$|\boldsymbol{r}_{\small {\rm A}}\times \boldsymbol{\hat{\Omega }}|$| may be computed for this case; this is indicated by dashes in Table 3. The conical Alfvén surface shapes arise from the presence of a significant azimuthal |$\boldsymbol{B}$|-component which increases the local |$u_{\small {\rm A}}$| value and thus shifts the local radial position of the Alfvén surface further away from the stellar surface. We note here that the Southern hemisphere of stars mapped with ZDI is usually incomplete (see e.g. F16) due to the observing geometry and that the tall southern Alfvén surfaces may not exist in reality.

4.3 Mass-loss and angular momentum loss

One of the strengths of numerical wind models as opposed to early, analytical models, which were often highly symmetrical and/or restricted to the equatorial plane, is the ease with which the fully three-dimensional wind mass-loss rate may be calculated. As we already have the density ρ and the wind velocity |$\boldsymbol{u}$| everywhere in the solution domain we find the wind mass-loss rate by integrating over a closed surface S containing the star:

The wind mass-loss of Solar-type stars is difficult to constrain by observations. Based on observations of ∼16 Solar-type stars Wood (2004) and Wood et al. (2005) found increasing mass-loss rates when going backwards in time from the Sun’s age to about ∼0.7 Gyr, before which they found lower rates of mass-loss. The existence of this ‘wind dividing line’ has, however, been called into question (Vidotto 2021; Wood et al. 2021). While |$\dot{M}$| is too small to significantly affect the mass of Solar-type stars over their lifetime, the parameter scales the angular momentum loss in semi-empirical models where |$\dot{J} \propto \dot{M} P_\text{rot}^m R_{\small {\rm A}}^n$|⁠.

The wind angular momentum loss, of course, is what slows the star’s spin with age. To compute this quantity we follow the general approach of Mestel (1999) and Vidotto et al. (2014a), but in a slight deviation from Paper I and Paper II (where this quantity was calculated at the Alfvén surface) we calculate the angular momentum loss rate over the same spherical surface S as used in (6), as this is numerically more convenient for very extended Alfvén surfaces such those of as 5|$\, \times \,$|BD-07 2388 and 5|$\, \times \,$|LO Peg. The two methods give comparable results. At a spherical surface S

The scalar angular velocity Ω = 2π/P_rot is found from Table 1, |$\boldsymbol{\varpi }= \boldsymbol{r} - (\boldsymbol{r} \cdot \boldsymbol{\hat{z}}) \boldsymbol{\hat{z}}$| is the cylindrical distance from the origin, and |$\boldsymbol{v} = \boldsymbol{u} - \boldsymbol{\Omega }\times \boldsymbol{r}$| is the wind velocity in the co-rotating frame. The notation (·)₃ refers to the z component of the enclosed vector. The values calculated for |$\dot{M}$| and |$\dot{J}$| are given in Table 3. The Ω term in equation (7) suggests that |$\dot{J}/\Omega$| may vary less than |$\dot{J}$| itself; we return to this topic in Section 5.

4.4 Wind pressure in the equatorial plane

The stellar wind pressure exerted on a body travelling with velocity |$\boldsymbol{w}$| through the wind is (Vidotto, Jardine & Helling 2011)

where the right-hand side terms are the thermal pressure, magnetic pressure, and ram pressure, respectively. Fig. 4 shows the wind pressure |$\text P_{\small W}$| exerted on a stationary (⁠|$\boldsymbol{w}=0$|⁠) object in the stellar equatorial plane extending out to the would-be orbit of an Earth-like planet (dotted white line). The intersection of the Alfvén surface (of Fig. 3) and the stellar equatorial plane is plotted with a black line. In Fig. 4, we observe as expected that |$\text P_{\small W}$| drops with increasing stellar distance, and that spiral arm-like corotating interaction regions (CIRs, see Belcher & Davis 1971; Gosling 1996) form when regions of fast stellar wind encounters regions of slow stellar winds. We mainly see two-armed structures corresponding to the dipole-dominated coronal magnetic fields in Fig. 2. We note that HD 6569, HIP 76768, BD-07 2388, TYC 5164, and HIP 12545, with their low magnetic inclinations, give less pronounced spiral structures than the highly inclined TYC 0486 and BD-16 351, with the other models occupying a middle ground. It is clear that the rapid rotators HII 739, BD-07 2388, PW And, and LO Peg produce more tightly wound spirals than the other, slower rotators.

The average wind pressure for an Earth-like planet is tabulated as |$\text P_{\small W}^{\oplus }$| in Table 3. This value is calculated by taking the average of equation (8) over both planetary and stellar phase. At such a large distance from the central star the ram pressure term |$\rho |\boldsymbol{u}- \boldsymbol{w}|^2$| dominates the thermal and magnetic pressure contributions to |$\text P_{\small W}$|⁠. Furthermore, the wind is mainly radial |$|\boldsymbol{u}| \approx u_r$|⁠, and the wind speed exceeds the planet’s orbital speed |$|\boldsymbol{u}| \approx u_r \gg |\boldsymbol{w}|$| so that |$\text P_{\small W}\approx \rho u_r^2$|⁠. For reference, monthly |$\text P^{\oplus }_{\small W}$| averages in the OMNI⁷ (see King & Papitashvili 2005) data set range from ∼0.1 to |${\sim }{10}\, {\rm nPa}$|⁠.

When the super-Alfvénic stellar wind encounters a magnetized body, such as a planet, a standing shock wave will form upstream of the body. The magnetospheric stand-off distance R_mag is a measure of the distance from the planet to the shock wave; this distance can be approximated by considering pressure balance between the wind pressure and the body’s dipolar magnetic field strength,

here, R_p is the radius of the planet (Chapman & Ferraro 1931; Vidotto et al. 2009b) and B_p is the strength of the planet’s magnetic field. For our Earth-like planet we use |$B_\text{p}={0.7}\, {\rm G}$|⁠; this value accounts for magnetospheric currents (Mead 1964). The calculated values of R_mag/R_p are given in Table 3. We see that the magnetospheric stand-off distance ranges from |$3.6 \, R_\text{p}$| to |$8.2 \, R_\text{p}$| which is past the suggested threshold of |${\sim }2\, R_\text{p}$| (of Lammer et al. 2007) where an Earth-like planet would be protected from atmospheric erosion. As R_p increases with stellar age, this suggests that an Earth-like planet would be protected from atmospheric erosion throughout its star’s main-sequence lifetime.

5 DISCUSSION

In this section, we examine trends in our results in the magnetic field strength as well as in other parameters. In Section 5.1, we study the effect of scaling the magnetic field on the wind mass-loss rate and other aggregate quantities. In Section 5.2, we quantify the correlation between magnetic field and wind parameters using ordinary least-squares analysis. In Section 5.3, we compare our computed values of wind mass-loss and angular momentum loss with literature values.

5.1 Magnetic field scaling

As the ZDI reported magnetic field strength is subject to uncertainty, it is of interest to study the effect of pure magnetic field scaling on relevant aggregate quantities of Table 2 and Table 3. We carry out such a study by comparing wind models of the same star using the differing magnetic scalings of 1 and 5 times the radial surface magnetic field, i.e. by comparing e.g. the wind mass-loss rate of the PELS 031 and 5|$\, \times \,$|PELS 031 model cases of the star PELS 031 in Table 1. Fig. 5 shows the effect of scaling the magnetic field by comparing the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series with the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series in this way. In addition to the Pleiades, AB Doradus, Columba, and β Pictoris stars modelled in this work, we have included models of the Hyades stars from Paper I and the Coma Berenices and Hercules-Lyra stars from Paper II (see Appendix  A). The symbols representing each model case of the same star are linked together by dashed line segments with equations

where the index i ranges over the stars in Table 1. The symbol shapes are the same as is used throughout the paper; a circle for the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series and a five-pointed star for the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series; the colours are the same as in Table 1. In keeping with Paper I and Paper II we refer to the dashed line segments and symbols of each star as its ‘barbell’.

While we found tight bounds on the slopes of the fitted lines so that the slopes α_i had similar values for each star in Paper II, we observe significantly more spread in the younger stars modelled in this work. The |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series models HD 6569, HIP 10272, TYC 0486, and their scaled counterparts essentially fall inside the min/max range spanned by the older stars of Paper I and Paper II; the remainder of the new stars modelled in this work fall outside Paper II limits and they do not display a clear linear trend. Instead, we observe a trend where the power-law indices (exponents) α_i decrease with increasing unsigned surface flux Φ. To highlight this trend we have plotted the slopes α_i of each fitted barbell with equation |$\hat{y}_i(x)\propto x^{\alpha _i}$| in the right-hand side of Fig. 5. The dashed lines indicate the Φ values at the two ends of the barbell, and the star symbols is placed at the geometric mean of the endpoints.

From the right-hand side of Fig. 5 we see that trend of decreasing α_i with increasing unsigned surface flux Φ is pronounced for the mass-loss rate |$\dot{M}$|⁠, the angular momentum loss rate |$\dot{J}$|⁠, and the scaled angular momentum loss rate |$\dot{J}/\Omega$|⁠, where α_i falls between 0.6 and 1.4. The range of variation is smaller for Φ_open although HII 739 and PELS 031 appear as outliers. Finally, for the wind pressure at an Earth-like planet |$\text P_{\small W}^{\oplus }$| there is no clear trend of α_i in Φ except that BD-07 2388, LO Peg, and PW And lie below the rest of the model values.

5.2 Statistical trends and correlations

In this section, we carry out a complementary analysis to that of Section 5.1; instead of considering pairs of models of the same star in the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series, we study the effect of the magnetic field strength, surface magnetic flux, and open magnetic flux on the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series separately.

By log-transforming our data the assumptions of ordinary least-squares analysis are satisfied (see e.g. Draper 1998). Fig. 6 shows trend lines of the form

fitted to the variation of mass-loss rate, angular momentum loss rate, rotation-scaled angular momentum loss rate, wind pressure for an Earth-like planet and magnetospheric stand-off distance for an Earth-like planet as a function of the average surface radial magnetic field strength, the unsigned surface magnetic flux, and the open magnetic flux (all these quantities are discussed in Section 4). In each of the panels of Fig. 6, there are two dashed line segments of the form of equation (11); the leftmost one corresponds to the fitted trend line of the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series, while the rightmost one corresponds to the fitted trend line of the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series. Surrounding each fitted trend line are two grey shaded regions; the innermost (dark grey) region is a 95 per cent confidence band of the a and b parameters of the fitted trend line, while the outermost (light grey) region corresponds to a 95 per cent prediction band. The prediction band variance is the sum of the confidence band variance and the regression mean square error. Given further stars and magnetograms drawn from similar populations as the ones in this study, there is a 95 per cent chance of them falling inside these light grey shaded regions. The confidence interval of a for the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series are also printed in the plots, for example we have |$\dot{M} \propto |B_r|^{0.81\pm 0.18}$| for the models of the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series.

The confidence intervals of a and b, and several other statistical parameters describing the fitted lines of Fig. 6 are given in Table 4. Using the numerical confidence intervals of a and b we can recover the full equation (11) of the fitted lines and confidence bands, e.g. |$\dot{M}=(10^{9.01\pm 0.24}\, {\rm kg\, s}^{-1})(|B_r|/{1}\, {\rm G})^{0.81\pm 0.18}$| for the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series |B_r| fitted line. Note that when reconstructing the full fitted equations of the form y = bx^a as in equation (11), the choice of physical units for x and y will affect the numerical value of b and its confidence interval. For the numerical values in Table 4 all quantities are in SI base units except that the magnetic field is expressed in Gauss |$({1}\, {\rm G} = 10^{-4}\, {\rm T})$|⁠, |$R_{\small {\rm A}}$| and |$|\boldsymbol{r_{\small {\rm A}}} \times \boldsymbol{\hat{\Omega }}|$| are expressed in stellar radii (see Table 1), and R_mag is expressed in terms of planetary radii R_p, as is the case throughout this paper. The magnetic fluxes Φ and Φ_open are expressed in Weber (⁠|${1}\, {\rm Wb}={1}\, {\rm T\, m}^{2} = 10^{8}\, {\rm G\, cm}^{2} = 10^{8}\, {\rm Mx}$|⁠).

Continuing with the columns of Table 4, the coefficients of determination (r² values) quantify the amount of variation that is explained by the fitted curve, the y_0.975/y_0.025 values are a physical measure of the amount of variation not explained by the fitted trend line, and the probability values (p values) give the likelihood that the observed trend is a spurious one that would not exist in a larger sample of stellar models. The confidence intervals, prediction bands, r² values, and p values are calculated in the standard way using the statsmodels (Seabold & Perktold 2010) Python package.

The most notable difference between Fig. 6 and the similar analysis of Paper II is that the Φ and Φ_open ranges are now overlapping between the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series due to the inclusion of the more powerful and complex magnetic fields of the stars in the AB Doradus moving group, the Pleiades cluster, Columba, and β Pictoris, as well as these stars’ rapid rotation. In general, we also observe stronger trends with smaller uncertainties (i.e. larger power-law indices with reduced error estimates) for the quantities plotted except for |$\dot{J}/\Omega$| against Φ_open where the fitted power-law indices are slightly smaller and accompanying error estimates are slightly larger. More specifically, in the leftmost column of Fig. 6, where the quantities are plotted against the mean surface radial field strength B_r, the trend lines from the current full set of models are generally steeper than the trend lines of Paper II with the exception of |$\dot{J}/\Omega$| against |B_r| which is similar to the previous results. This also occurs in the middle column. The fits for |$\dot{J}/\Omega$| are, however, very similar to those of Paper II. In the rightmost column, where the model output Φ_open is on the x-axis, we see that |$\dot{M}$| and |$\dot{J}$| exhibit stronger trends than in Paper II, similarly to what was seen in the middle column. The fits to |$\dot{J}/\Omega$| against Φ_open have slightly weaker trends, while the trends in |$\text P_{\small W}^{\oplus }$| and R_mag are stronger.

The model HII 739 and its scaled counterpart 5|$\, \times \,$|HII 739 are clear outliers in the top three rows of Fig. 6; in the left column of panels this can be attributable to this star’s large radius of |$1.5\, R_{\odot }$|⁠. HII 739 and 5|$\, \times \,$|HII 739 remain high in the centre and right columns but their x position is shifted rightwards because of Φ’s dependence on |$R_{\mathchoice{}{}{\scriptscriptstyle }{}\bigstar }$|⁠. The very young star models TYC 6349 and 5|$\, \times \,$|TYC 6349 also exhibit large values of |$\dot{J}/\Omega$|⁠. The other notable outliers, possibly due to their star’s rapid |$P_\text{rot}={0.3}\, {\rm d}$| rotation, are BD-07 2388 and 5|$\, \times \,$|BD-07 2388. We do not, however, see a similar effect for the other rapid rotator in the LO Peg and 5|$\, \times \,$|LO Peg models.

It is interesting to compare the mass-loss rate |$\dot{M}$| panels in Fig. 6 (top row) where |$\dot{M}$| is plotted against the average unsigned radial magnetic field strength |B_r| (left), the unsigned magnetic flux Φ (middle) and the open unsigned magnetic flux Φ_open (right). The four stars BD-07 2388, LO Peg, PW And, and TYC 6878, and their counterparts in the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| give the appearance of a weak saturation effect limiting the increase of |$\dot{M}$| as a function of |B_r|, but the effect is not seen in the panel showing |$\dot{M}$| as a function of Φ (although the four stars still lay below the main trend line). This apparent saturation limiting |$\dot{M}$| with increasing |$|\boldsymbol{B}|$| can be seen more clearly in Fig. 7 where we compare the results of our wind modelling with literature values. We emphasize to the reader that the apparent saturation limiting |$\dot{M}$| as a function of increased surface magnetism in our data set is a spurious effect arising from differences in stellar radius (see Table 1), and that it is not present when plotting |$\dot{M}$| against |$\Phi =4 \pi R^2 \, |B_r|$| and thus accounting for where differences in stellar radii.

In Appendix  B, Fig. B1 and Table B1 show the results of a similar analysis where the model cases in the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| have been pooled into a single series of 60 stellar wind models.

5.3 Comparison with literature values

In Fig. 7, we compare the full set of wind mass-loss rates and wind angular momentum loss rates described in this work and Paper I and Paper II with literature values obtained from three-dimensional wind modelling and from scaling laws. As the choices and assumptions behind this plot are laid out in Paper I, we provide this plot here to emphasize that the trends that we observed in Paper I and Paper II still apply. When taken together with the blue symbols representing awsom type three-dimensional wind models, we see an increasing trend in mass-loss rate up to about |$3 \times 10^{10}\, {\rm kg\, s^{-1}}$| or ∼24 times the fiducial Solar mass-loss rate of |$2 \times 10^{-14}\, {\rm M_{\odot } \, yr^{-1}}$| (e.g. Vidotto 2021) in the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series of models, and a similar increasing trend up to about 1 × 10¹¹ kg s⁻¹ or ∼80 times the Solar value for the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series. As we noted near the end of Section 5.2 the apparent threshold of |$\dot{M}$| values at |${\sim } 80 \dot{M}_{\odot }$| is mostly due to stellar radius variation (see Table 1) in our sample.

Our calculated values of |$\dot{M}$| are within the spread of the values calculated by See et al. (2019) (see Fig. 7 and see Paper I for the differences between the CS11 and mod M15 methods). We note, however, that in our estimate of mean |$|\boldsymbol{B}|$| from their work (see Paper I) does not yield mean |$|\boldsymbol{B}|$| values past ∼200 G such that many of the strongest magnetic fields in the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series have no comparison in their work. The frequently cited work of Wood et al. (2002, 2005) also predicts mass-loss rates up to |${\sim } 100\, \dot{M}_{\odot }$| for 0.7 Gyr old stars, but lower mass-loss rates for even younger stars. We do not find strong evidence of a drop in |$\dot{M}$| values for the youngest stars in our data set but our results are otherwise well matched to those of their work.

In the bottom panel of Fig. 7 we compare our calculated angular momentum loss rates to those of the literature. The right-hand side axis shows the angular momentum loss rate in terms of the average Solar angular momentum loss rate, with the caveat that the average Solar wind angular momentum loss rate is not as well constrained as the average Solar wind mass-loss rate. We use the observation-based angular momentum loss rate of 2.2 × 10²³ N m from Finley et al. (2019) as our Solar baseline.

As can be seen from the red-outlined symbols in Fig 7, numerical wind models give different |$\dot{M}$| and |$\dot{J}$| values depending on the adopted Solar magnetic field strength (i.e. phase of the Solar cycle). As in Paper I and Paper II we continue to find excellent agreement between our models and the awsom based wind models of Alvarado-Gómez et al. (2016a) and Pognan et al. (2018), and good agreement with the models of Réville et al. (2016).

6 CONCLUSIONS

In this work, we have modelled the winds of 15 young, Solar-type stars with well-constrained ages in the Pleiades cluster, the AB Doradus moving group, the Columba association, and the β Pictoris association. These stars aged 24 to 125 Myr have a wide range of rotation periods from 0.3 to 7 d; this is expected for stars in this age range, and matches the observations of e.g. Gallet & Bouvier (2013). The models are driven using observationally based magnetic maps derived with ZDI.

We have studied the coronal magnetic field, Alfvén surface, and wind pressure out to 1 au of our models and we find mostly dipole-dominated coronal magnetic fields; in our models these appear if ∼30 per cent or more of the surface magnetic energy is dipolar. Our ZDI maps have a wide range of inclinations of the magnetic dipole, and produce a wide range of Alfvén surface lobe and current sheet orientations. The effective magnetogram degree measures in ℓ_.99 (see Table 2) does not seem to greatly affect the shape of the Alfvén surface in Fig. 3 with the possible exception of the conical Alfvén surfaces of BD-07 2388, LO Peg, 5|$\, \times \,$|BD-07 2388, and 5|$\, \times \,$|LO Peg. We do not consider these shapes to result from the ℓ_.99 values however, but rather from the two stars’ very rapid rotational periods which also leads to more detailed ZDI magnetograms as Section 2.1 briefly addressed.

In order to study the isolated effect of the surface magnetic field strength independently of the other model parameters, we created the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series where all magnetic field maps are scaled by a factor of 5. Such an increase of the magnetic field strength is generally accompanied by an increase in the mass- and angular momentum loss rate and wind pressure for an Earth-like planet as was seen in Section 5.1. The comparison of the wind models in the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series, including the models originally presented in Evensberget et al. (2021, 2022) and re-processed as part of this work (see Section 3), showed that for this wide range of surface magnetic field strengths there is no single scaling law of the form y ∝ |B_r|^α or y ∝ Φ^α that is a good fit to the data in the left panel of Fig. 5. Instead, as can be seen in the right panel of Fig. 5, the power-law index α appears to decrease with increasing magnetic field strength.

In Section 5.2, we fitted power laws and prediction bands to the |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| and |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series separately. This kind of fitted curve differs from the barbell curves in Section 5.1 as it includes the effect of other stellar parameter that are known to vary with the surface magnetic field, most notably the stellar rotation rate Ω. We find somewhat wider prediction bands than in Evensberget et al. (2022), suggesting that the more rapid rotation of the stars modelled in this work yield slightly different characteristics than the older stars of Evensberget et al. (2021, 2022); this can be due to either Ω itself or the influence of Ω on the effective magnetogram degree. In general, we find tighter correlations between the unsigned surface flux Φ = 4πR²|B_r| and our parameters of interest than between |B_r| and our parameters of interest. This is expected as Φ accommodates stellar radius variations (from stellar type, see Table 1) on the prediction band widths. The |$\dot{J} / \Omega$| parameter has a significantly tighter correlation with the magnetic field strength than the |$\dot{J}$| parameter itself, which is expected from the dependence of |$\dot{J}$| on Ω in eq. (7).

In Section 5.3, we observe that our results agree well with other models using the awsom model, but predict lower values of |$\dot{M}$| and |$\dot{J}$| than models using polytropic MHD models. Only for the strongest magnetic fields in our |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| model series do we obtain comparable values of |$\dot{M}$| and |$\dot{J}$| to polytropic MHD models with significantly lower surface magnetic fields. As was previously noted in Evensberget et al. (2022) this gap may possibly be closed by applying an age- or rotation-based scaling to awsom model parameters such as the Poynting flux-to-field ratio |$\Pi _{\small {\rm A}}/B$|⁠, which was found by Boro Saikia et al. (2020) to affect the mass-loss rate such that |$\dot{M} \propto \Pi _{\small {\rm A}}/ B$| (see also Kavanagh et al. 2021).

In this work and Evensberget et al. (2021, 2022), we have created what we believe to be the largest set of ZDI-driven three-dimensional stellar wind models to date. Between the models all physical and numerical parameters are the same except the surface magnetic maps, and three scalar stellar parameters mass, radius, and rotation rate. This consistency has enabled us to formulate robust scaling relations that includes uncertainty estimates and give predictive ranges for values such as wind mass- and angular momentum loss rates for stars in the 0.02 to 0.7 Gyr range. We hope that this data set and our derived scaling relations will be useful in comparing numerical wind models with observational constraints on |$\dot{M}$|⁠, for creating more detailed models of wind–planet interactions, and to give better constraints on how Solar parameters values are applicable to young, Solar-type stars in magnetohydrodynamic modelling.

ACKNOWLEDGEMENTS

For part of this work, DE has been funded by a University of Southern Queensland (UniSQ) International Stipend Research Scholarship (ISRS) and a UniSQ International Fees Research Scholarship (IFRS). Parts of this project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 817540, ASTROFLOW). Parts of this material is based upon work supported by the National Science Foundation under Grant No. (AST-2108985). Parts of this research was undertaken using the University of Southern Queensland Fawkes HPC which is co-sponsored by the Queensland Cyber Infrastructure Foundation (QCIF), see www.unisq.edu.au/hpc. We thank SURF (www.surf.nl) for the support in using the National Supercomputer Snellius. This work utilizes data obtained by the Global Oscillation Network Group (GONG) Program, managed by the National Solar Observatory, which is operated by AURA, Inc. under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Interamerican Observatory. This work has used the Vienna Atomic Line Database (VALD), operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. We acknowledge the use of the SIMBAD data base. This research has used the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in A&AS 143, 23. This research has used NASA’s Astrophysics Data System. This work was carried out using the swmf tools developed at The University of Michigan Center for Space Environment Modelling (CSEM) and made available through the NASA Community Coordinated Modelling Center (CCMC). The open-source version of swmf is available at https://github.com/MSTEM-QUDA/SWMF.git. Version 2.40.a5beb110f (2022-10-08) has been used for all processing. This work has used the following additional numerical software, statistics software and visualization software: numpy (van der Walt, Colbert & Varoquaux 2011) version 1.23.5, scipy (Virtanen et al. 2020) version 1.9.3, matplotlib (Hunter 2007) version 3.7.1, statsmodels (Seabold & Perktold 2010) version 0.13.2, tecplot version 2022.1.0.14449, and pytecplot version 1.5.0. We would like to express our gratitude to the anonymous reviewers for their valuable contributions, which have enhanced the quality and rigour of our article.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: REPROCESSED AGGREGATE PARAMETERS

In Section 3, we noted that the Hyades (Paper I) and Coma Berenices and Hercules-Lyra (Evensberget et al. 2022) models have been re-processed with the same numerical configuration as the models presented in this work. In Table A1, we provide the aggregate parameters (in the same format as Table 3) resulting from the re-processing. In the re-processing we found only minor differences to the results in Evensberget et al. (2021, 2022) so that the overall conclusions of this work and Evensberget et al. (2021, 2022) were not affected by the re-processing.

Note that the Solar maximum/minimum models were based on magnetograms from the Global Oscillation Network Group (GONG) spherical harmonic transform coefficients⁸ at Carrington Rotation 2157 (maximum) and 2211 (minimum) and are thus referred to as Sun-G2157 and Sun-G2211 in Paper I. 

APPENDIX B: POOLED SERIES

In Evensberget et al. (2021), which considered five stars in the Hyades, the same analysis as in Section 5.2 was given for a single, pooled data series that comprises |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series models and the |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series models. Thus, in order to maintain consistency with Evensberget et al. (2021, 2022), we provide here the full pooled series analysis which includes the 30 |$B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series models and the 30 |$5B_{\mathchoice{}{}{\scriptscriptstyle }{}\mathrm{ZDI}}$| series models, for a total of 60 wind models. The fitted trend lines and prediction bands are shown in Fig. B1 with the corresponding data in Table B1.

For a large range of B_r values, such as what is presented here, it should be remembered that the radial magnetic field strength is correlated with the rotation period.