The active lives of stars: A complete description of the rotation and XUV evolution of F, G, K, and M dwarfs

C. P. Johnstone; M. Bartel; M. Güdel

doi:10.1051/0004-6361/202038407

Home

All issues

Volume 649 (May 2021)

A&A, 649 (2021) A96

Full HTML

Free Access

Issue		A&A Volume 649, May 2021


Article Number		A96
Number of page(s)		26
Section		Stellar structure and evolution
DOI		https://doi.org/10.1051/0004-6361/202038407
Published online		20 May 2021

A&A 649, A96 (2021)

The active lives of stars: A complete description of the rotation and XUV evolution of F, G, K, and M dwarfs^⋆

C. P. Johnstone¹^,2, M. Bartel² and M. Güdel²

¹ Natural History Museum, Burgring 7, 1010 Vienna, Austria
e-mail: colinjohnstone@gmail.com
² University of Vienna, Department of Astrophysics, Türkenschanzstrasse 17, 1180 Vienna, Austria

Received: 12 May 2020
Accepted: 28 August 2020

Abstract

Aims. We study the evolution of the rotation and the high energy X-ray, extreme ultraviolet (EUV), and Ly-α emission for F, G, K, and M dwarfs, with masses between 0.1 and 1.2 M_⊙, and provide a freely available set of evolutionary tracks for use in planetary atmosphere studies.

Methods. We develop a physical rotational evolution model constrained by observed rotation distributions in young stellar clusters. Using rotation, X-ray, EUV, and Ly-α measurements, we derive empirical relations for the dependences of high energy emission on stellar parameters. Our description of X-ray evolution is validated using measurements of X-ray distributions in young clusters.

Results. A star’s X-ray, EUV, and Ly-α evolution is determined by its mass and initial rotation rate, with initial rotation being less important for lower mass stars. At all ages, solar mass stars are significantly more X-ray luminous than lower mass stars and stars that are born as rapid rotators remain highly active longer than those born as slow rotators. At all evolutionary stages, habitable zone planets receive higher X-ray and EUV fluxes when orbiting lower mass stars due to their longer evolutionary timescales. The rates of flares follow similar evolutionary trends with higher mass stars flaring more often than lower mass stars at all ages, though habitable zone planets are likely influenced by flares more when orbiting lower mass stars.

Conclusions. Our results show that single decay laws are insufficient to describe stellar activity evolution and highlight the need for a more comprehensive description based on the evolution of rotation that also includes the effects of short-term variability. Planets at similar orbital distances from their host stars receive significantly more X-ray and EUV energy over their lifetimes when orbiting higher mass stars. The common belief that M dwarfs are more X-ray and EUV active than G dwarfs is justified only when considering the fluxes received by planets with similar effective temperatures, such as those in the habitable zone.

Key words: stars: activity / stars: chromospheres / stars: coronae / stars: late-type / stars: magnetic field / stars: rotation

^⋆

Our code for calculating stellar rotation and XUV evolution is available at https://github.com/ColinPhilipJohnstone/Mors and an extensive grid of evolutionary tracks is available at https://zenodo.org/record/4266670#.X6rMuq4o9H5.

© ESO 2021

1. Introduction

Stellar magnetic activity is crucially important for the formation and evolution of planetary atmospheres and surface conditions. Magnetic fields cause the outer atmospheres of stars to be heated to million degree temperatures, leading to the emission of X-ray and ultraviolet (XUV) radiation as well as the acceleration of magnetised stellar winds that have a diverse range of effects on planets and their atmospheres. This high energy radiation is absorbed in the upper atmospheres of planets leading to dissociation, ionisation, and heating (Roble et al. 1988; Murray-Clay et al. 2009), which enhances and drives atmospheric loss processes, including rapid hydrodynamic losses if the star’s activity is high enough (Tian et al. 2005). The details of how the activities of stars evolve are important for determining the eventual state of a planet’s atmosphere (Luger et al. 2015; Johnstone et al. 2015b; Kubyshkina et al. 2019).

When stars are first born, they have very strong magnetic fields and high XUV luminosities (Yang & Johns-Krull 2011), which both decay rapidly during the first few million years (Gregory et al. 2012, 2016). At later ages, rotational spin-down causes a significant decay in magnetic fields and high energy emission (Güdel et al. 1997; Vidotto et al. 2014). For the XUV emission, this decay is more rapid at shorter wavelengths, causing the shape of the XUV spectrum to evolve (Ribas et al. 2005; Claire et al. 2012). The timescale for activity evolution is a strong function of mass, with lower mass stars remaining active for much longer (West et al. 2008; Reiners & Basri 2008). It is known empirically that a star’s X-ray emission is determined by its mass, age, and rotation rate (Pizzolato et al. 2003), meaning that its activity evolution is linked to its rotational evolution. Stars born as fast rotators remain active longer than those born as slow rotators (Tu et al. 2015). Stars with different initial rotation rates can have very different influences on how the atmospheres of planets evolve (Johnstone et al. 2015b; Johnstone 2020).

A star’s long-term rotational evolution depends primarily on its mass and initial rotation rate (Gallet & Bouvier 2015; Matt et al. 2015), though other factors play a role, such as the lifetime of the early circumstellar gas disk during the classical T Tauri phase (Herbst et al. 2002) and the star’s metallicity (Amard & Matt 2020). At ages of ∼1 Myr, the rotation rates of stars are distributed from approximately a few to a few ten times the rate of the current Sun (Affer et al. 2013). In the first few million years (Myr), this rotation distribution appears to remain approximately constant with no clear evolution (Rebull et al. 2004), except possibly for stars with masses below 0.4 M_⊙ (Henderson & Stassun 2012). After this phase, pre-main-sequence contraction causes them to spin up until the zero-age main-sequence (ZAMS), and during this time the initial wide distribution of rotation rates gets even wider. After the ZAMS, angular momentum removal by stellar winds causes stars to spin down and the distribution converges until most stars follow a single mass- and age-dependent value. How long this convergence takes depends on stellar mass, with solar mass stars converging in the first billion years and lower mass stars taking longer; for fully convective M dwarfs, it is possible that no convergence takes place (Irwin et al. 2011). At ages of a few billion years (Gyr), lower mass stars tend to be slower rotators than higher mass stars (Nielsen et al. 2013). A detailed review is available in Bouvier et al. (2014).

In this paper, we study the evolution of stellar rotation and XUV emission on the pre-main-sequence and the main-sequence for stars with masses between 0.1 and 1.2 M_⊙ and provide a comprehensive empirical description of these phenomena. We explore how long stars with different masses and initial rotation rates remain active and discuss the influences on the evolution of planetary atmospheres. In Sect. 2, we develop a description of rotational evolution, in Sect. 3, we study the evolution of X-ray emission, in Sect. 4, we study extreme ultraviolet (EUV) and Ly-α emission, in Sect. 5, we discuss XUV fluxes in the habitable zone, in Sect. 6, we discuss the contributions of flares, and in Sect. 7, we summarise our results. In Appendix A, we describe our freely available grid of rotation and XUV evolution models.

2. Rotational evolution

2.1. Rotational evolution model

In this paper, we extend our rotation model developed in previous studies (Johnstone et al. 2015a, 2019b; Tu et al. 2015) to describe the full rotational evolution between 1 Myr and the end of the main-sequence for stars with masses between 0.1 and 1.2 M_⊙. Our model describes the star as being composed of an envelope and a core, with the envelope corresponding to the outer convective zone and the core corresponding to everything else, and these two components are assumed to rotate as solid bodies and at separate rates. The rotation rates of the two components are influenced by core-envelope angular momentum exchanges and changes in the internal structure of the star, and the rotation of the envelope is additionally influenced by angular momentum loss by a magnetised stellar wind. We use the stellar evolution models of Spada et al. (2013) to get all parameters related to the star’s internal structure, including convective turnover times, and use their set of models with initial metallicities, compositions, and mixing length parameters that best match the case of the Sun. Throughout this paper, we define the Sun’s rotation rate, Ω_⊙, as 2.67 × 10⁻⁶ rad s⁻¹.

To understand how a star’s rotation evolves, we must consider the both the changes in the angular momentum values and the moments of interia of the core and the envelope. In our model, this is given by

$\begin{matrix} \frac{d Ω_{core}}{d t} = \frac{1}{I_{core}} (- τ_{ce} - τ_{cg} - Ω_{core} \frac{d I_{core}}{d t}), \end{matrix}$ $\begin{aligned} \frac{\mathrm{d}\Omega _\mathrm{core} }{\mathrm{d}t} = \frac{1}{I_\mathrm{core} } \left( -\tau _\mathrm{ce} - \tau _\mathrm{cg} - \Omega _\mathrm{core} \frac{\mathrm{d} I_\mathrm{core} }{\mathrm{d}t} \right), \end{aligned}$ (1)

$\begin{matrix} \frac{d Ω_{env}}{d t} = \frac{1}{I_{env}} (τ_{w} + τ_{ce} + τ_{cg} + τ_{d} l - Ω_{env} \frac{d I_{env}}{d t}), \end{matrix}$ $\begin{aligned} \frac{{{\rm{d}}{\Omega _{{\rm{env}}}}}}{{{\rm{d}}t}} = \frac{1}{{{I_{{\rm{env}}}}}}\left( {{\tau _{\rm{w}}} + {\tau _{{\rm{ce}}}} + {\tau _{{\rm{cg}}}} + {\tau _{\rm{d}}}l - {\Omega _{{\rm{env}}}}\frac{{{\rm{d}}{I_{{\rm{env}}}}}}{{{\rm{d}}t}}} \right), \end{aligned}$ (2)

where Ω_core and Ω_env are the core and envelope angular velocities, I_core and I_env are the core and envelope moments of inertia, τ_w is the stellar wind spin-down torque, τ_ce is the core-envelope coupling torque, τ_cg the core-growth torque, and τ_dl the disk-locking torque. For fully convective stars (M_⋆ ≲ 0.35M_⊙), the distinction between core and envelope is not considered and the stars are assumed to rotate entirely as solid bodies. For these stars, we replace the above equations with

$\begin{matrix} \frac{d Ω_{⋆}}{d t} = \frac{1}{I_{⋆}} (τ_{w} + τ_{dl} - Ω_{⋆} \frac{d I_{⋆}}{d t}), \end{matrix}$ $\begin{aligned} \frac{\mathrm{d}\Omega _\star }{\mathrm{d}t} = \frac{1}{I_\star } \left( \tau _\mathrm{w} + \tau _\mathrm{dl} - \Omega _\star \frac{\mathrm{d} I_\star }{\mathrm{d}t} \right), \end{aligned}$ (3)

where Ω_⋆ and I_⋆ are the star’s rotation rate and moment of inertia.

As in previous work, we calculate the stellar wind torque using

$\begin{matrix} τ_{w} = - K_{τ} τ^{'}, \end{matrix}$ $\begin{aligned} \tau _\mathrm{w} = - K_\tau \tau ^{\prime }, \end{aligned}$ (4)

where K_τ is a parameter that we use to better reproduce the Skumanich spin-down of the modern Sun and we use K_τ = 11 as derived in Sect. 4.1 of Johnstone et al. (2015a). This value depends on the specific wind torque formula used and our assumption of the solar wind mass loss rate and the Sun’s dipole field strength. We calculate τ′ using

$\begin{matrix} τ^{'} = K_{1}^{2} B_{dip}^{4 m} {\dot{M}}_{⋆}^{1 - 2 m} R_{⋆}^{4 m + 2} \frac{Ω_{env}}{{(K_{2}^{2} v_{esc}^{2} + Ω_{env}^{2} R_{⋆}^{2})}^{m}}, \end{matrix}$ $\begin{aligned} \tau ^{\prime } = K_1^2 B_\mathrm{dip} ^{4m} \dot{M}_\star ^{1-2m} R_\star ^{4m+2} \frac{\Omega _\mathrm{env} }{(K_2^2 v_\mathrm{esc} ^2 + \Omega _\mathrm{env} ^2 R_\star ^2)^m}, \end{aligned}$ (5)

where B_dip is the strength of the dipole component of the star’s magnetic field, Ṁ_⋆ is the wind mass loss rate, R_⋆ and M_⋆ are the star’s radius and mass, and v_esc is the surface escape velocity ( $= \sqrt{2 G M_{⋆} / R_{⋆}}$ $=\sqrt{2 G M_\star / R_\star}$ ). The other parameters in this equation are given by K₁ = 1.3, K₂ = 0.0506, and m = 0.2177. Matt et al. (2012) derived this relation from a large number of magnetohydrodynamic wind simulations assuming axis-symmetry, a dipole magnetic field, and a polytropic equation of state for the heating with approximately constant wind temperatures. We do not consider the influence of non-dipolar magnetic field geometries, although this can have an influence on the wind torque (Holzwarth & Jardine 2005; Réville et al. 2015; Garraffo et al. 2016). Stellar magnetic fields are mostly composed of a combination of dipole and higher order components and since the dipole component rapidly becomes dominant further from the star’s surface, it tends to be the dominant component for angular momentum loss even for stars with very non-dipolar fields (Finley & Matt 2018). However, since the thermal structure of the wind influenced the angular momentum loss (Cohen 2017; Pantolmos & Matt 2017), we should expect some deviation from Eq. (5) due to the importance of the magnetic field geometry on the wind heating and acceleration and due to the observed dependence of stellar coronal temperatures on activity level (Johnstone & Güdel 2015). Additionally, See et al. (2019) argued that non-dipolar geometries could be more important for stars with very high mass loss rates. For our purposes, it is not necessary to consider these details.

Both the dipole field strength and the mass loss rate in Eq. (5) are manifestations of the star’s magnetic activity, which depends on the star’s mass, age, and rotation rate. Many activity parameters are tightly correlated with the Rossby number, Ro, which is a dimensionless parameter defined by Ro = P_rot/τ_c, where P_rot is the rotation period and τ_c is the convective turnover time. Based on the observational correlation between global magnetic field and rotation given by Vidotto et al. (2014), we use

$\begin{matrix} B_{dip} = {\begin{matrix} B_{dip, ⊙} {(\frac{R o_{sat}}{R o_{⊙}})}^{- 1.32}, & if R o \leq R o_{sat}, \\ B_{dip, ⊙} {(\frac{Ro}{R o_{⊙}})}^{- 1.32}, & otherwise, \end{matrix} \end{matrix}$ $\begin{aligned} B_\mathrm{dip} = \left\{ \begin{array}{ll} B_{\mathrm{dip} ,\odot } \left( \frac{Ro_\mathrm{sat} }{Ro_\odot } \right)^{-1.32},&\mathrm{if}\ Ro \le Ro_\mathrm{sat} ,\\ B_{\mathrm{dip} ,\odot } \left( \frac{Ro}{Ro_\odot } \right)^{-1.32},&\mathrm{otherwise},\\ \end{array} \right. \end{aligned}$ (6)

where Ro_⊙ and B_dip, ⊙ are the Rossby number and dipole field strength of the modern Sun and Ro_sat is the Rossby number of the saturation threshold. For B_dip, ⊙, we take the value 1.35 G (Johnstone et al. 2015a). Very little is known about the mass loss rates of low mass stars other than the modern Sun and our best option is to assume a plausible scaling law for the mass loss rate, Ṁ_⋆, given by

$\begin{matrix} {\dot{M}}_{⋆} = {\begin{matrix} f {\dot{M}}_{⊙} {(\frac{R_{⋆}}{R_{⊙}})}^{2} {(\frac{R o_{sat}}{R o_{⊙}})}^{a_{w}} {(\frac{M_{⋆}}{M_{⊙}})}^{b_{w}}, & if R o \leq R o_{sat}, \\ f {\dot{M}}_{⊙} {(\frac{R_{⋆}}{R_{⊙}})}^{2} {(\frac{Ro}{R o_{⊙}})}^{a_{w}} {(\frac{M_{⋆}}{M_{⊙}})}^{b_{w}}, & otherwise, \end{matrix} \end{matrix}$ $\begin{aligned} \dot{M}_\star = \left\{ \begin{array}{ll} f \dot{M}_\odot \left( \frac{R_\star }{R_\odot } \right)^2 \left( \frac{Ro_\mathrm{sat} }{Ro_\odot } \right)^{a_\mathrm{w} } \left( \frac{M_\star }{M_\odot } \right)^{b_\mathrm{w} },&\mathrm{if} \ Ro \le Ro_\mathrm{sat} ,\\ f \dot{M}_\odot \left( \frac{R_\star }{R_\odot } \right)^2 \left( \frac{Ro}{Ro_\odot } \right)^{a_\mathrm{w} } \left( \frac{M_\star }{M_\odot } \right)^{b_\mathrm{w} },&\mathrm{otherwise} ,\\ \end{array} \right. \end{aligned}$ (7)

where Ṁ_⊙ = 1.4 × 10⁻¹⁴ M_⊙ yr⁻¹ is the current Sun’s mass loss rate, a_w and b_w are fit parameters, and f is the magneto-centrifugal factor discussed below. As described in Appendix B, we find a_w = −1.76 and b_w = 0.649. We use the Rossby number instead of the rotation rate in the above relations since this captures the effect of young pre-main-sequence stars being saturated even at low rotation rates observed for activity related phenomena such as X-ray emission and magnetic field strength (Briggs et al. 2007; Johnstone et al. 2014). Based on activity indicators such as X-ray luminosity, we assume that saturation occurs at a constant Rossby number for all stellar masses and ages and is the same for all activity related phenomena, allowing us to derive Ro_sat from X-ray observations. The value of Ro_sat depends on the convective turnover times used, and in Sect. 3.1 we find Ro_sat = 0.0605 for the convective turnover times given by Spada et al. (2013).

For very rapidly rotating stars, we include also magneto-centrifugal effects. As winds propagate away from their stars, azimuthal components form for both the propagation velocity and the magnetic field, causing additional radial acceleration by centrifugal and Lorentz forces (Weber & Davis 1967). While these forces are negligible for the modern solar wind, they can be significant for the winds of rapidly rotating stars with strong magnetic fields (Belcher & MacGregor 1976; Johnstone 2017). When a star’s rotation is close to the break-up rotation rate, its mass loss is enhanced by these effects, likely causing additional spin-down which should stop pre-main-sequence stas from spinning up to or past break-up. We therefore include in Eq. (7) an additional multiplicative factor given by

$\begin{matrix} f (Ω_{env}) = {\begin{matrix} 1, & if x \leq 0.1, \\ 0.93 {(1.01 - x)}^{- 0.43} e^{0.31 x^{7.5}}, & otherwise, \end{matrix} \end{matrix}$ $\begin{aligned} f \left( \Omega _\mathrm{env} \right) = \left\{ \begin{array}{ll} 1,&\mathrm{if}\ x \le 0.1,\\ 0.93 \left( 1.01 - x \right)^{-0.43} e^{0.31 x^{7.5}},&\mathrm{otherwise},\\ \end{array} \right. \end{aligned}$ (8)

where x = Ω_env/Ω_break. This equation is based on the 1D MHD simulations of Johnstone (2017) and is shown in Fig. 1. As they spin up, stars reach break-up when the Keplerian co-rotation radius becomes equal to the star’s radius at the equator, which happens at a rotation rate of $Ω_{break} = {(G M_{⋆} / R_{e}^{3})}^{1 / 2}$ $\Omega_\mathrm{break} = \left( G M_\star / R_\mathrm{e}^3 \right)^{1/2}$ , where R_e is the star’s equatorial radius at the break-up threshold. Defining R_p as the star’s polar radius at the break-up threshold, bulging at the equators of very rapidly rotating stars means that R_e = 1.5R_p (Maeder 2009), which leads to

$\begin{matrix} Ω_{break} = {(\frac{2}{3})}^{\frac{3}{2}} {(\frac{G M_{⋆}}{R_{p}^{3}})}^{\frac{1}{2}} . \end{matrix}$ $\begin{aligned} \Omega _\mathrm{break} = \left( \frac{2}{3} \right)^{\frac{3}{2}} \left( \frac{G M_\star }{R_\mathrm{p} ^3} \right)^{\frac{1}{2}}. \end{aligned}$ (9)

Fig. 1.

Upper panel: factor f in Eq. (7), as a function of rotation with the blue circles showing the results of MHD simulations for rapidly rotating stars by Johnstone (2017) and the dashed black line showing our fit formula given by Eq. (8). Lower panel: evolution of the break-up rotation rate for different stellar masses.

We assume here that R_p = R_⋆ where R_⋆ is the star’s radius in the absence of rotation and therefore do not consider the influence of rapid rotation on the polar radius, which is anyway likely to change R_p by only a few percent (Ekström et al. 2008). The evolution of Ω_break is shown in Fig. 1.

For the exchanges of angular momentum between the core and the envelope, we adopt the approach used in MacGregor & Brenner (1991), Spada et al. (2011), and Gallet & Bouvier (2015), where the torque is given by

$\begin{matrix} τ_{ce} = \frac{Δ J}{t_{ce}}, \end{matrix}$ $\begin{aligned} \tau _\mathrm{ce} = \frac{\Delta J}{ t_\mathrm{ce} }, \end{aligned}$ (10)

where t_ce is the core-envelope coupling timescale and ΔJ is the angular momentum that must be transferred between the two components in order to make them rotate with the same speed. We define this torque such that positive values imply angular momentum transfer from the core to the envelope, meaning that ΔJ is given by

$\begin{matrix} Δ J = \frac{I_{env} I_{core}}{I_{env} + I_{core}} (Ω_{core} - Ω_{env}), \end{matrix}$ $\begin{aligned} \Delta J = \frac{ I_\mathrm{env} I_\mathrm{core} }{ I_\mathrm{env} + I_\mathrm{core} } \left( \Omega _\mathrm{core} - \Omega _\mathrm{env} \right), \end{aligned}$ (11)

which implies that ΔJ = 0 when Ω_core = Ω_env. For t_ce, we assume

$\begin{matrix} t_{ce} = a_{ce} {(| Ω_{env} - Ω_{core} |)}^{b_{ce}} {(\frac{M_{⋆}}{M_{⊙}})}^{c_{ce}}, \end{matrix}$ $\begin{aligned} t_\mathrm{ce} = a_\mathrm{ce} \left( |\Omega _\mathrm{env} - \Omega _\mathrm{core} | \right)^{b_\mathrm{ce} } \left( \frac{M_\star }{M_\odot } \right)^{c_\mathrm{ce} }, \end{aligned}$ (12)

where a_ce, b_ce, and c_ce are fit parameters. This is similar to the assumption made by Spada et al. (2011) and Johnstone et al. (2019b) with the additional mass dependence. As described in Appendix B, we find a_ce = 25.6, b_ce = −3.25 × 10⁻², and c_ce = −0.448 when Ω_env and Ω_core have units of Ω_⊙ and t_ce has units of Myr.

The core-growth torque, τ_cg, in Eqs. (1) and (2) represents a different type of angular momentum exchange between the core and the envelope. The growth of the core at young ages means that material that is part of the envelope becomes part of the core, and this material possesses angular momentum that is therefore transported from the envelope to the core. Assuming a positive value corresponds to angular momentum transport from the envelope to the core, this is given by

$\begin{matrix} τ_{cg} = - \frac{2}{3} R_{core}^{2} Ω_{env} \frac{d M_{core}}{d t}, \end{matrix}$ $\begin{aligned} \tau _\mathrm{cg} = - \frac{2}{3} R_\mathrm{core} ^2 \Omega _\mathrm{env} \frac{\mathrm{d}M_\mathrm{core} }{\mathrm{d}t}, \end{aligned}$ (13)

where M_core and R_core are the core mass and radius. During the pre-main-sequence phase, both I_core and I_env change rapidly due to two effects: firstly, the growth of the core increases I_core and decreases I_env, and secondly, the radial distribution of mass changes. This torque balances the contribution of the former effect to both dI_env/dt and dI_core/dt in Eqs. (1) and (2), meaning that during the core growth phase, no unphysically rapid spin-down of the core takes place as I_core increases. The above is valid when dM_core/dt > 0, and when dM_core/dt < 0, Ω_env should be replaced by Ω_core.

The final ingredient in this model is disk-locking, which is an ad hoc assumption that is commonly used to reproduce the lack of observed spin-up during the early pre-main-sequence when stars still possess circumstellar gas disks (e.g. Allain 1998; Gallet & Bouvier 2013). This lack of spin-up is surprising since these stars are contracting and accreting high angular momentum material from their disks, meaning that something must be removing angular momentum from the star during this phase. This could be enhanced stellar winds driven by the accretion of material from the disk onto the stellar surface (Matt & Pudritz 2008; Cranmer 2009). It is normal in rotation models to simply keep the star’s surface rotation rate constant for the first few million years of the star’s life. To do this, we assume a disk-locking torque acting on the envelope that cancels all other terms in Eq. (2), given by

$\begin{matrix} τ_{dl} = {\begin{matrix} - τ_{w} - τ_{ce} - τ_{cg} + Ω_{env} \frac{d I_{env}}{d t}, & if t \leq t_{disk}, \\ 0, & otherwise . \end{matrix} \end{matrix}$ $\begin{aligned} \tau _\mathrm{dl} = \left\{ \begin{array}{ll} -\tau _\mathrm{w} - \tau _\mathrm{ce} - \tau _\mathrm{cg} + \Omega _\mathrm{env} \frac{\mathrm{d} I_\mathrm{env} }{\mathrm{d}t},&\mathrm{if}\ t \le t_\mathrm{disk} ,\\ 0,&\mathrm{otherwise}.\\ \end{array} \right. \end{aligned}$ (14)

where t_disk is the disk-locking time. As in Tu et al. (2015) and Johnstone et al. (2019b), we assume

$\begin{matrix} t_{disk} = 13.5 {(\frac{Ω_{0}}{Ω_{⊙}})}^{- 0.5}, \end{matrix}$ $\begin{aligned} t_\mathrm{disk} = 13.5 \left( \frac{ \Omega _0 }{ \Omega _\odot } \right)^{-0.5}, \end{aligned}$ (15)

where Ω₀ is the initial (1 Myr) rotation rate of the star in units of Ω_⊙ and t_disk is in Myr. The inverse dependence means that the envelopes of fast rotators start spinning up earlier, which is consistent with the observed distribution of fast rotators in the young ∼12 Myr old cluster h Per (Moraux et al. 2013). For Ω₀ below Ω_⊙, this equation can predict unreasonably large values of t_disk and to avoid this, we suggest setting a maximum value of 15 Myr, though we do not consider any cases with Ω₀ < Ω_⊙ in this paper.

The five unconstrained parameters in our models are a_w and b_w from Eq. (7), and a_ce, b_ce, and c_ce from Eq. (12). We determine these parameters using a large number of observational constraints, described in the next section, and a Markov chain Monte Carlo (MCMC) method for parameter optimisation, described in Appendix B. We obtain a_w = −1.76, b_w = 0.649, a_ce = 25.6, b_ce = −3.25 × 10⁻², and c_ce = −0.448. While these results could yield important information about stellar wind properties and angular momentum transport within stellar interiors, this is not the aim of this paper and we do not attempt physical interpretations of our parameter fitting results.

2.2. Observational constraints

To constrain the parameters in our model, we use measured rotation distributions in young clusters. We assume that the sequence of clusters provides a good representation of the rotational evolution of individual groups of stars, which is reasonable given that multiple clusters with similar ages have similar distributions (Fig. 14 of Hartman et al. 2010). A large number of the rotation rates that we use were presented in Sect. 3 of Johnstone et al. (2015a) who collected data for seven clusters with ages between 100 and 1000 Myr: these were the Pleiades (∼125 Myr), M50 (∼130 Myr), M35 (∼150 Myr), NGC 2516 (∼150 Myr), M37 (∼550 Myr), Praesepe (∼580 Myr), and NGC 6811 (∼1000 Myr). We include additional newer rotation measurements for Pleiades, Praesepe, and NGC 6811. We also include an additional five clusters, with rotation rates mostly measured by K2, with ages between 12 and 2500 Myr. These are h Per (∼12 Myr), NGC 2547 (∼40 Myr), Hyades (∼600 Myr), NGC 752 (∼1300 Myr), and NGC 6819 (∼2500 Myr). We put the twelve clusters into seven age bins, with ages of 12, 40, 150, 550, 650, 1000, and 2500 Myr, as summarised in Table 1 with references to the original sources for the rotation period measurements. For most clusters, we use stellar masses reported in the original rotation studies and rederive masses for Pleiades, M35, M37, and NGC 6819 using colour-mass conversions from the data given by Pecaut & Mamajek (2013).

Table 1.

Sources of the observational constraints on rotational evolution used in this paper.

It is necessary in our fitting procedure also to have information about the rotation rates of stars with ages later than our oldest cluster, which has an age of ∼2.5 Gyr. One possibility is to use rotation measurements from K2 observations of the ∼4.0 Gyr cluster M67. However, inconsistent results for the rotation distribution in this cluster have been found (Barnes et al. 2016; Gonzalez 2016a,b), and Esselstein et al. (2018) showed that accurately determining the rotation periods from the available data is challenging, especially given their long 20–35 day rotation periods and the limited 75 day observation window. Another problem with using M67 is that the age of the cluster is uncertain, and many age determinations are based on the gyrochonology method, which for our purpose is not useful. This is also a problem with using the large number of rotation periods measured for field stars. In this paper, we use the gyrochronological relation given by Mamajek & Hillenbrand (2008) to estimate rotation periods at 4.5 Gyr for each mass. This is given by

$\begin{matrix} P_{rot} = a {[{(B - V)}_{0} - c]}^{b} t^{n} \end{matrix}$ $\begin{aligned} P_\mathrm{rot} = a \left[ (B-V)_0 - c \right]^b t^n \end{aligned}$ (16)

where a = 0.407, b = 0.325, c = 0.495, n = 0.566, t is the age in Myr, and P_rot is in days. We convert from (B − V)₀ to M_⋆ using the data from Pecaut & Mamajek (2013).

In Fig. 2, we show as red circles the observed rotation distribution as a function of mass in each observed age bin. We additionally add for comparison the 2 Myr cluster NGC 6530 with rotation rates and masses determined by Henderson & Stassun (2012). We also show in the lower right panel the rotation distribution for approximately 40 000 field stars observed by Kepler (black circles) using the rotation rates determined by McQuillan et al. (2014) and Santos et al. (2019) and for approximately 320 mid and late M dwarfs (red stars) from the MEarth Project (Newton et al. 2016). Neither NGC 6530 nor the Kepler and MEarth field stars are used in the parameter determination for our rotational evolution model.

Fig. 2.

Stellar rotation distribution at each observed age bin listed in Table 1. Red and blue circles show observed and modelled distributions, with the modelled distributions calculated by evolving the observed 150 Myr distribution between 1 and 5000 Myr. In the 150 Myr panel, red circles show stars that we are unable to fit with our rotation model using realistic initial rotation rates. In the lower right panel, black show the distribution measured by Kepler (McQuillan et al. 2014; Santos et al. 2019), red stars show the distribution determined by the MEarth Project (Newton et al. 2016), blue and green circles show the 1 and 5 Gyr modelled distributions, and solid lines show predictions from the gyrochrological relation of Mamajek & Hillenbrand (2008).

2.3. Results: The evolution of rotation

In Fig. 2, we show the observed 150 Myr distribution evolved between 2 and 2500 Myr using our rotation model at each observed age bin. We use the 150 Myr cluster since it has the most stars and these stars are distributed over the entire mass range. The cluster evolution shown in Fig. 2 is as expected, with an initially wide distribution that becomes wider during the pre-main-sequence spin-up phase and then converges to a mass- and age-dependent value on the main-sequence as stars spin down.

The 2 Myr panel in Fig. 2 shows the starting distribution implied by the 150 Myr cluster and our best fit rotation model. This fits well the observed distribution in the ∼2 Myr cluster NGC 6530, though it appears that between 0.3 and 0.5 M_⊙, the observed distribution has more fast rotators. Our 2 Myr distribution is approximately mass-independent down to 0.3 M_⊙, with stars distributed between 1 and 40 Ω_⊙. At lower masses, the distribution becomes much tighter and shifted towards faster rotation, which has good observational support from the 40 Myr old cluster NGC 2547. This could have two intepretations: either the initial rotation distribution is mass-dependent at such low masses in the way that we show, or the initial distribution is mass-independent and the disk-locking times are shorter for such low masses. The latter is supported by the results of Henderson & Stassun (2012), who found evidence for a mass-independent initial rotation distribution at 1–2 Myr followed by mass-dependent spin-up of stars with masses below 0.5 M_⊙ in the first 10 Myr (shown in their Fig. 16). Since the activities of all of these stars are saturated, both possibilities lead to the same early activity evolution.

A feature clearly visible in Fig. 2 is the lack of observed spin-down between 650 Myr and 1000 Myr for stars with masses between 0.5 and 0.9 M_⊙ despite clear spin-down of solar mass stars, as has been discussed in the recent literature (Agüeros et al. 2018; Curtis et al. 2019). This epoch of stalled spin-down must be temporary since rotation measurements of field stars do not show a lack of slowly rotating K dwarfs. This lack of spin-down is not reproduced in our model, or any other rotation model to our knowledge, suggesting that additional physical processes are needed to explain the early main-sequence spin-down of K dwarfs. This could be a temporary reduction in angular momentum loss, due maybe to changes in the global magnetic field structures and strengths, or potentially a change in core-envelope coupling timescales causing more rapid angular momentum transfer from the core to the envelope. Interestingly, the 550 Myr age bin also shows a similar effect.

In the 150 Myr panel, the red circles show observed stars not included in our model distribution. While evolutionary tracks can be fit to all slow rotators in the 150 Myr distribution, we do not consider stars that require initial rotation rates below 1Ω_⊙ since there is no observational support for such slow rotation at 1 Myr. For the rapid rotators, there is a difficulty that we cannot fit realistic evolutionary tracks for some ∼1 M_⊙ stars, many of which are above the break-up threshold and therefore unrealistic. While some stars are below break-up, they must have undergone significant spin-down since the ZAMS and our models suggest that they could not have arrived at their current rotation rates without being above break-up in the past. The 550 and 650 Myr clusters also show stars above the upper bound of our model distribution. A similar difficulty was found in the models of Matt et al. (2015), who suggested that such stars require a modified torque, and this could suggest that at young ages, the torques for stars with the same basic parameters (mass, age, rotation) are not uniform. It is also possible that these stars have or had short-period binary or planetary companions or it could be that the measured rotation periods are unrealistically short, which can happen for specific distributions of surface inhomogeneities.

In the lower right panel of Fig. 2, we show our model distributions at 1 and 5 Gyr. At 5 Gyr, the rotation distribution fits very well the expectations from the gyrochonological relation of Mamajek & Hillenbrand (2008) at masses above 0.4 M_⊙ but shows a dip towards slow rotation for lower masses. The lower bound of the main cloud of stars in the Kepler distribution fits well these 5 Gyr expectations for masses above 0.6 M_⊙, but is above our expectations at lower masses. This could suggest that our description of the later evolution of low mass stars is unrealistic, but it could also be that the Kepler sample is missing the slowest rotators since their longer rotation periods and lower photometric variability make period determinations difficult (McQuillan et al. 2014). The latter interpretation is supported by the rotation distribution from the MEarth Project which has many M dwarfs with periods slower than 100 days and is more consistent with our model predictions. In our model, this dip is due to a peak at this mass in the convective turnover times, τ_c, that we use which causes more rapid spin-down at masses around 0.3 M_⊙. This feature, realistic or not, does not influence the results of this paper since we use the same τ_c values to calculate X-ray emission and the peak in τ_c influences those calculations in such a way the feature is not present in our model X-ray distributions.

In Fig. 3, we show rotation tracks for stellar masses of 1.0, 0.75, 0.5, and 0.25 M_⊙. In each bin, we show tracks for slow, medium, and fast rotators, defined as the tracks that pass through the 5th, 50th, and 95th percentiles of the observed 150 Myr rotation distribution. These percentiles are calculated using all stars within 0.1 M_⊙ of the specified masses and we use the distribution at 150 Myr since that is the most complete and has the largest number of stars. For the 1.0, 0.75, and 0.5 M_⊙ cases, the tracks clearly give an excellent description of the observed rotational evolution sequence, with the only issue being the slow rotator track being slower than the observed 5th percentile in the first two clusters for the 1.0 and 0.75 M_⊙ cases. For the 0.25 M_⊙ case, the solid lines are the tracks defined, as above, to go through the corresponding percentiles of the observed distribution at 150 Myr and the dashed lines are defined to go through the percentiles in NGC 6530. In neither case, do these tracks give a good description at all ages of the evolution implied by the observed percentiles. Since the saturation threshold for activity is at very low rotation rates for fully convective M dwarfs, this is only an uncertainty for the X-ray evolution after ∼2 Gyr. The rotational evolution of these low mass stars, and the possibility that some fully convective M dwarfs do not spin down significantly, will be studied in more detail in Bartel et al. (in prep).

Fig. 3.

Rotational evolution for stars with masses of 1.0, 0.75, 0.5, and 0.25 M_⊙. The red, green, and blue lines show our slow, medium, and fast rotator tracks with the solid and dotted lines representing the envelope and core rotation rates. The slow, medium, and fast rotator are defined as the observed 5th, 50th, and 95th percentiles at 150 Myr. The grey circles show observed rotation rates and the short horizontal lines show the percentiles from the measurements in each age bin. In the 0.25 M_⊙ mass bin, the dashed lines show the tracks that instead pass through the observed percentiles in the 2 Myr cluster NGC 6530.

3. X-ray evolution

3.1. X-ray relations

It is known empirically that a star’s X-ray luminosity depends on its rotation rate, mass, and age and that this dependence can be broken down into the unsaturated and saturated regimes. For main-sequence solar mass stars, the saturation threshold is at a rotation rate of approximately 15Ω_⊙ and for lower mass stars, this threshold is at slower rotation rates, meaning that K and M dwarfs need to spin down more than G dwarfs before they enter the unsaturated regime. For pre-main-sequence stars, the saturation threshold is at much slower rotation, and in the first 10 Myr, all stars are saturated regardless of their rotation rates. The saturation X-ray luminosity is approximately a constant fraction of the bolometric luminosity meaning that among saturated stars, higher mass stars are more X-ray luminous and pre-main-sequence stars tend to be more X-ray luminous than their main-sequence counterparts. A final factor is a significant spread in measured X-ray luminosities for stars with similar parameters, much or all of which is a result of random and cyclic variability.

The X-ray dependence on stellar parameters can be well described with a power-law dependence between the X-ray luminosity normalised to the bolometric luminosity, R_X = L_X/L_bol, and the Rossby number, Ro. In the unsaturated regime, the dependence can also be described by a single mass independence relation between L_X and rotation rate, and Reiners et al. (2014) argued that this description is preferable. In this paper, we use the Ro–R_X relation since it is able to describe the observed evolution of the saturation threshold for pre-main-sequence stars. This relation can be described as a broken power-law given by

$\begin{matrix} R_{X} = {\begin{matrix} C_{1} R o^{β_{1}}, & if R o \geq R o_{sat}, \\ C_{2} R o^{β_{2}}, & if R o \leq R o_{sat}, \end{matrix} \end{matrix}$ $\begin{aligned} R_\mathrm{X} = \left\{ \begin{array}{ll} C_1 Ro^{\beta _1},&\mathrm{if}\ Ro \ge Ro_\mathrm{sat} ,\\ C_2 Ro^{\beta _2},&\mathrm{if}\ Ro \le Ro_\mathrm{sat} ,\\ \end{array} \right. \end{aligned}$ (17)

where C₁, C₂, β₁, and β₂ are constants to be determined empirically and Ro_sat is the saturation Rossby number. It is common in the literature to assume that R_X has a constant value of R_X, sat in the saturated regime, meaning that β₁ = 0 and C₁ = R_X, sat, though a weak dependence of R_X on Ro has been pointed out in the literature (Reiners et al. 2014; Magaudda et al. 2020). We make no assumption for β₁, but fit it to the observed data. We do not consider the ‘supersaturation’ phenomenon, which is a possible decrease in X-ray emission for the most rapid rotators in the saturated regime (Jardine 2004; Argiroffi et al. 2016).

As described in Appendix C, we constrain the constants in the above equation using the distribution of stars with known X-ray luminosities and rotation periods presented by Wright et al. (2011). The values of these constants depend on the model or relation used to get the convective turnover times, τ_conv, and therefore to be consistent with our rotational evolution model, we use the values of τ_conv from the stellar evolution models of Spada et al. (2013). The result for our relation is β₁ = −0.135 ± 0.030, β₂ = −1.889 ± 0.079, Ro_sat = 0.0605 ± 0.00331, and R_X, sat = 5.135 × 10⁻⁴ ± 3.320 × 10⁻⁵, where R_X, sat is the value of R_X at Ro_sat. The values of C₁ and C₂ can be derived from the fact that the two power-laws have equal values of R_X at the saturation point, meaning that $R_{X, sat} = C_{1} R o_{sat}^{α_{1}} = C_{2} R o_{sat}^{α_{2}}$ $R_\mathrm{X,sat} = C_1 Ro_\mathrm{sat}^{\alpha_1} = C_2 Ro_\mathrm{sat}^{\alpha_2}$ . This relation is shown in Fig. 4 and the evolution of the saturation rotation rate and X-ray luminosity implied by our fit is shown in Fig. 5. Our relation is shallower in the unsaturated regime than many previous estimates and this shallow relation suggests that the Sun is less X-ray active than other stars with similar parameters, which is consistent with results from other activity indicators (Reinhold et al. 2020).

Fig. 4.

Stellar R_X, defined as L_X/L_bol, as a function of Rossby number for the sample of main sequence stars collected by Wright et al. (2011). The Rossby numbers were calculated using the convective turnover times given by Spada et al. (2013). The dashed red line shows our best fit relation, given by Eq. (17). The green line shows the range of values for the non-flaring Sun, with the green circles showing the 10th and 90th percentiles of the Sun’s X-ray luminosity.

Fig. 5.

Convective turnover time (upper panel), the rotation rate of the saturation threshold (middle panel), and the X-ray luminosity at the saturation threshold (lower panel) as functions of age for different stellar masses. In the middle and lower panels, the transition from solid to dashed lines shows the time in the evolution for each mass when 90% of all stars have dropped below the saturation threshold. The convective turnover times are from Spada et al. (2013) and the saturation threshold is calculated assuming a saturation Rossby number of 0.0605.

Since the stellar sample of Wright et al. (2011) contains only main-sequence stars and does not contain unsaturated fully convective M dwarfs, it is useful to check if the above scaling law applies on the pre-main-sequence and for all stellar masses, especially when combined with the convective turnover times that we use. Most very young stellar clusters are not useful since the large convective turnover times put most stars in the saturated regime. A useful cluster is h Per since its 12 Myr age means that we might expect some slowly rotating solar mass stars to be unsaturated. Argiroffi et al. (2016) studied X-ray emission in h Per and found evidence of a relation between Ro and R_X for slowly rotating solar mass stars. While their results suggest a much shallower power law relation than the relation found for main-sequence stars, this is based on a very small number of stars with X-ray detections and if we rederive the Rossby numbers using the convective turnover times from Spada et al. (2013), we find that the distribution for h Per appears consistent with the distribution for main-sequence stars, as we show in the upper panel of Fig. 6.

Fig. 6.

Upper panel: comparison of the Ro–R_X distribution shown in Fig. 4 with the stars in the ∼12 Myr cluster h Per with the black circles and green triangles showing X-ray detections and upper limits from Argiroffi et al. (2016). We recalculate the Rossby numbers using the convective turnover times given by Spada et al. (2013) assuming an age of 12 Myr. Lower panel: comparison of the Ro–R_X distribution shown in Fig. 4 with a sample of fully convective M dwarfs from Jeffries et al. (2011) and Wright et al. (2018), with detections and upper limits shown as circles and triangles.

For fully convective M dwarfs, Jeffries et al. (2011) showed that stars with masses below 0.35 M_⊙ follow approximately the same Ro–R_X relation as higher mass stars, and this result has been supported by later studies (Wright & Drake 2016; Wright et al. 2018; Magaudda et al. 2020). In the lower panel of Fig. 6, we show the Ro–R_X distribution for fully convective main-sequence M dwarfs from several studies with Ro recalculated using the convective turnover times from Spada et al. (2013). The low activity M dwarfs look mostly consistent with the entire sample, though it is notable that so many of the lowest activity M dwarfs are upper limits which could put them below the Ro–R_X distribution. This could be consistent with the results of Magaudda et al. (2020) who found a steeper slope in the unsaturated regime for fully convective stars. We conclude that it is appropriate with our current knowledge to use the Ro–R_X relation constrained for main-sequence stars on the pre-main-sequence and for all stellar masses.

The uncertainties given for the fit parameters in Eq. (17) are calculated using the bootstrap method and do not take into account uncertainties in the measured X-ray luminosities and rotation periods, uncertainties in the modelled τ_conv, biases in our sample of stars, and the method used to fit the relation. Several previous studies have estimated β₂ and a range of values have been derived. Considering the entire sample, Wright et al. (2011) found β₂ = −2.18 and from a subset of the sample they found β₂ = −2.7. The difference between our estimate and their lower estimate is mostly because we use different convective turnover times. Similarly, Reiners et al. (2014) estimated β₂ = −2 and showed that this depends on which stars are used in the fit and which method is used for the fitting, with the OLS(Y|X) method consistently estimating shallower relations than the OLS Bisector method. We use the former method for reasons discussed in Appendix C.

3.2. Variability and the nature of the spread

One reason why the parameters in Eq. (17) are difficult to determine is the large spread of approximately two orders of magnitude in the Ro–R_X distribution. In Fig. 7, we show the distribution of Δlog R_X, defined as the difference between the measured log R_X and the value predicted by Eq. (17). The distribution can be described as a normal distribution centred around zero with a standard deviation of 0.359, or as a generalised normal distribution as described in Fig. 7. There are three factors that could cause this spread and it is important to understand the contributions of each. Firstly, inaccuracies in measurements of the X-ray luminosities should cause some spread, and the sample we use is especially susceptible to this since it contains X-ray measurements from several instruments, especially Einstein IPC, ROSAT PSPC, and XMM-Newton. Each instrument is sensitive to a different photon energy range and the count-rate to flux conversions used to calculate fluxes in the entire X-ray range are often derived using different methods. Secondly, random and cyclic variability surely play a substantial role in causing the spread in the Ro–R_X relation. The magnitude of the variability for the Sun, not including flares, is shown in green in Fig. 4 and is similar to but smaller than the overall spread. Thirdly, the spread could be caused by intrinsic differences between stars, with some stars being intrinsically more active that other stars with similar masses, ages, and rotation rates.

Fig. 7.

Upper panel: histogram showing the distribution of Δ log R_X, defined as the difference between the log of the observed R_X and the log of the R_X from our fit formula. The black line shows a normal distribution with μ = 0 and σ = 0.359 and the red line shows the generalised normal distribution with a probability density function given by βexp[−(|Δlog R_X−μ|/α)^β]/[2αΓ(1/β)], where β = 1.43, μ = 0, α = 0.4, and Γ is the gamma-function. Lower panel: cumulative distribution functions for the variability amplitudes of several samples of stellar L_X measurements. The black line represents the spread in the Ro–R_X distribution, the dotted line represents variability of the Sun, and the other lines show variability distributions derived from the literature as described in Sect. 3.2.

Several studies have estimated the X-ray variability of individual stars. Stern et al. (1995) and Schmitt et al. (1995) combined measurements from Einstein and ROSAT to study variability in the Hyades and nearby K and M dwarfs on approximately decadal timescales and found that L_X typically varies by less than a factor of two. Some of this variability is cyclic in nature since activity cycles have been seen in other stars (Boro Saikia et al. 2018) and much of the variability is a result of stellar flares (Güdel et al. 2004; Stelzer et al. 2007), though for our purposes, the exact nature of the variability is not important. Marino et al. (2000, 2002) studied variability on timescales of a few months for M, K, G, and F dwarfs with multiple ROSAT measurements and found that on these timescales, M dwarfs are more variable than their higher mass counterparts. Marino et al. (2003) found similar results among a more active sample of stars in the Pleiades. By comparing with Einstein measurements, Marino et al. (2000) found larger variability for M dwarf on 8–15 year timescales. These studies did not find evidence that the variability amplitude depends on activity level or age, consistent with the fact that the spread in the Ro–R_X distribution is also independent of Ro. They defined the variability amplitude for stars with more than one X-ray luminosity measurement as |log L_X, 1 − log L_X, 2| and since they considered only variability of individual stars, this is equivalent to |log R_X, 1 − log R_X, 2|.

To compare these results to the spread in the Ro–R_X distribution, we calculate a distribution of variability amplitudes from the sample of stars used to constrain the Ro–R_X relation. For every star in this sample, we calculate |log R_X, 1 − log R_X, 2|, where log R_X, 1 is the value for that star and log R_X, 2 is the value for the star with the nearest Rossby number. The reason to use R_X here is that this corrects for different bolometric luminosities. We also calculate a variability amplitude distribution for the six stars in the Sun-in-time sample using L_X measurements by Güdel et al. (1997) and Telleschi et al. (2005). These measurements, summarised in Table 1 of Telleschi et al. (2005), were based on ROSAT and XMM-Newton observations taken 7–10 years apart. Finally, we also calculate the solar X-ray variability using daily average L_X values derived from the < 10 nm part of the quiescent (non-flaring) solar spectra provided by the Flare Irradiance Spectrum Model (Chamberlin et al. 2007) which uses activity proxies to reconstruct the full solar spectrum. For the solar distribution, we calculate the variability amplitudes for 10 000 pairs of days randomly selected from the time period 1948–2010.

In Fig. 7, we show the variability amplitude distribution for the sample of stars used to constrain the Ro–R_X relation and compare it to measured variability distributions from Marino et al. (2000, 2002). The blue and green lines show variability on monthly timescales for M dwarf and K, G, and F dwarfs, and these distributions are clearly inconsistent with the spread in the Ro–R_X distribution. The red line shows variability for M dwarfs on the longer timescale of 8–15 years and this matches the spread distribution much better, suggesting that the spread in the Ro–R_X relation is consistent with variability on longer timescales. The solar variability is consistent with the spread for small variability amplitudes but is inconsistent with the high variability tail of the spread, possibly because we use non-flaring solar spectra only. Strangely, the much smaller Sun-in-time sample shows far less variability than all other samples which could be a result of having only six stars in the sample. Using Kolmogorov-Smirnov tests, we find that the consistency between the spread in the Ro–R_X distribution and the variability in each of the other distributions can be rejected, except for the long-term M dwarf variability which cannot be rejected.

In conclusion, while the nature of the observed spread in the Ro–R_X distribution remains unclear, it appears unnecessary to assume that the spread suggests intrinsic differences between the X-ray activities of stars with similar masses, ages, and rotation rates. Despite the uncertainties, it is reasonable to assume that the spread is mostly caused by stellar variability on timescales of decades or longer. If true, stars have L_X values more than twice the average 20% of the time, more than three times the average 9% of the time, and more than five times the average 3% of the time, and they spend similar amounts of time with these factors below the average. These variations could be important when considering the effects of stellar XUV emission on the long-term evolution of planetary atmospheres.

3.3. Observations of X-ray evolution in young clusters

To verify the description of X-ray evolution presented in the next section, we use measured X-ray distributions from 10 young clusters with known ages collected from the literature. These are Taurus (∼2 Myr), h Per (∼12 Myr), NGC 2547 (∼40 Myr), α Per (∼50 Myr), Blanco I (∼50 Myr), Pleiades (∼125 Myr), NGC 2516 (∼150 Myr), NGC 6475 (∼300 Myr), M37 (∼550 Myr), and Hyades (∼650 Myr). Combining α Per and Blanco I into one bin at 50 Myr and NGC 2516 and Pleiades into one bin at 150 Myr gives us eight age bins with measured X-ray distributions. These age bins and the references for our sources of the X-ray measurements are given in Table 2 and the distributions are shown in Fig. 8. For six of these clusters, we use the X-ray luminosities determined by Núñez & Agüeros (2016) instead of the values from the original studies: these clusters and the original studies are NGC 2547 (Jeffries et al. 2006), NGC 2516 (Pillitteri et al. 2006), Pleiades (Stauffer et al. 1994; Micela et al. 1999; Stelzer et al. 2000), NGC 6475 (Prosser et al. 1995; James & Jeffries 1997), and Hyades (Stern et al. 1995; Stelzer et al. 2000). For an additional comparison, we include also X-ray measurements for nearby field stars from ROSAT collected from the NEXXUS database (Schmitt & Liefke 2004), converting between (B-V)₀ and mass using the data from Pecaut & Mamajek (2013).

Fig. 8.

Comparison between stellar X-ray luminosity distributions for several young clusters and our predicted distributions at these ages. In each panel, the blue circles shows our model distribution, the red circles show measured values, and the green triangles show upper limits. The final panel shows X-ray luminosities determined by ROSAT for nearby field stars (red circles) compared to our model distributions at ages of 1 and 5 Gyr.

Table 2.

Sources of the observational constraints on X-ray evolution used in this paper.

3.4. Results: X-ray evolution

Combining our rotation models with the relation for X-ray emission in the previous section allows us to understand the evolution of stellar X-ray emission. In Sect. 2.3, we used our rotational evolution model to evolve the observed 150 Myr rotation distribution from 1 Myr to 5 Gyr, which allows us to derive a model distribution at the ages of the clusters with observed X-ray luminosity distributions. For each star in the model rotation distribution, we predict an average log R_X using Eq. (17) and add to that a random value selected using a normal distribution that represents the spread around the best fit relation. In Fig. 8, blue circles show our model distributions for each observed cluster, and red circles and green triangles show the measured distributions for the young clusters described above. We find very good agreement, providing important validation of our approach. The comparison is difficult however due to observational limitations, especially due to detection thresholds for many clusters being at high X-ray luminosities. This is especially a problem for low mass stars since they are typically less X-ray luminous.

While measurements exist for a large number of older field stars, using these to test the evolution of X-ray emission is difficult given the lack of known ages for these stars. As a demonstration of the later X-ray evolution, we show our model L_X distribution at 1 and 5 Gyr in the lower right panel of Fig. 8 and compare these to ROSAT measurements of nearby field stars. Most of the field star measurements are contained well within the 1 and 5 Gyr distributions and the trend of higher L_X for higher mass stars is clearly visible, though as expected, some outliers are present given the large number of stars, the fact that not all stars in the sample are indeed within this age range, and the likely presence of tight binary systems.

In Fig. 9, we summarise the evolution of rotation (left panels) and X-ray emission (middle panels) for stellar masses between 0.1 and 1.2 M_⊙. In each panel, the red, green, and blue lines show the slow, medium, and fast rotator tracks, defined as the 5th, 50th, and 95th percentiles of the 150 Myr rotation distribution, and the dashed line shows the saturation threshold. At the start of the evolutionary sequence, there is a large spread in rotation rates at all masses, but this does not translate into a large spread in emission since the saturation threshold is at such low rotation rates that there are no unsaturated stars. There is a strong positive dependence of L_X on stellar mass due to the fact that higher mass stars have higher bolometric luminosities. At 10 Myr, the situation is similar, but the L_X values of all stars have decreased due to the decrease in bolometric luminosity on the pre-main-sequence. Some of the most slowly rotating stars with masses similar to that of the Sun have fallen out of saturation and have L_X values slightly below the saturation threshold. This is due to the increase in the saturation threshold rotation rate and not to rotational evolution. By 100 Myr, the L_X of the entire distribution has decreased and most solar mass stars have become unsaturated, leading to an increased spread in L_X values for these masses. These changes in the L_X distribution happen despite the fact that stars spin up on the pre-main-sequence and are due to the evolution of the decrease in the bolometric luminosities of stars and an increase in the rotation rate of the saturation threshold (driven by decreasing convective turnover times). After 100 Myr, the stellar bolometric luminosities and saturation threshold rotation rates remain approximately constant and most changes are a result of rotational evolution.

Fig. 9.

Evolution of the distributions of rotation (left-column), X-ray luminosity (middle-column), and X-ray flux in the habitable zone (right-column) against stellar mass, with each row showing the ages labelled in the left column. In each panel, the black circles show the model distribution, the red, green, and blue lines show the values for the slow, medium, and fast rotator cases, and the dashed black line shows the saturation threshold values. In the middle column, the insets show R_X against Ro for the model distribution, with the red line showing the empirical relation given by Eq. (17). The X-ray flux in the habitable zone is defined as the flux at an orbital distance half-way between the moist and maximum greenhouse limits and these limits are calculated at all ages assuming the stellar properties at 5 Gyr.

As stars spin down, a larger fraction of them enter the unsaturated regime, with higher mass stars tending to become unsaturated earlier than lower mass stars due to a combination of two effects. Firstly, due primarily to the lower rotation rates for the saturation threshold, rapidly rotating lower mass stars feel a smaller spin-down torque than higher mass stars with similar rotation rates, and therefore remain rapidly rotating longer. Secondly, the lower X-ray saturation threshold rotation rate for lower mass stars means that they must spin down more in order to become unsaturated. By 1 Gyr, most stars with masses > 0.4 M_⊙ and almost all stars with masses > 0.6 M_⊙ are unsaturated. By 5 Gyr, almost all stars are unsaturated, except those with masses of ∼0.1 M_⊙, and the convergence of the rotation distribution means that the initial rotation rate is no longer a factor influencing the emission. The dependence of the time that stars remain in the saturated regime on stellar mass and initial rotation rate is shown in Fig. 10.

Fig. 10.

Age at which a star falls below the saturation threshold as a function of stellar mass for slow (red), medium (green), and fast (blue) rotators. The background circles show when each of the stars in our model distribution become unsaturated. The y-axis is logarithmic for ages up to 1000 Myr and linear at later ages.

In Fig. 11, we show slow, medium, and fast rotator evolutionary tracks for L_X for several stellar masses. The shaded areas around each line represent the spread around the best fit relation between Rossby number and X-ray emission and extends one standard deviation of this spread above and below the average lines. Assuming that this spread is a result of activity variability only, stars spend approximately 90% of the time within these shaded regions. In the solar mass case, the X-ray emission of the slow and fast rotator tracks start to diverge in the first 10 Myr and are different by approximately a factor of ∼50 for most of the first ∼500 Myr. For the lower mass cases, this difference is smaller due to the saturation threshold being at a much lower rotation rate. For the 0.5 M_⊙ and 0.25 M_⊙ cases, the slow and fast rotator tracks are almost identical until 200 Myr and 2 Gyr respectively.

Fig. 11.

Evolutionary tracks for stellar X-ray luminosity for slow, medium, and fast rotators with several masses. The shaded areas around each line represents the spread around the best fit relation, showing one standard deviation of the spread above and below the average.

In Fig. 12, we show as functions of mass the ages at which the average L_X decays below values of 10²⁸, 10²⁹, and 10³⁰ erg s⁻¹. These can be seen as measures of the activity lifetimes for stars with different masses and initial rotation rates. Two separate regimes can be seen in the figure: at higher stellar masses, slow rotators cross the thresholds earlier than rapid rotators, and at lower stellar masses, all three rotation tracks cross the thresholds at approximately the same age. The reason for this difference is that higher mass stars decay across these thresholds due to rotational spin-down whereas the less luminous lower mass stars cross the thresholds due to the decay in their bolometric luminosities while still in the saturated regime. Likely the most unexpected feature of these results is the fact that when using L_X to measure activity lifetimes, higher mass stars remain active longer than lower mass stars.

Fig. 12.

Ages of stars when they fall below threshold values for their X-ray luminosity (left panel) and X-ray flux in the habitable zone (right panel) as functions of stellar mass. In both panels, each colour represents a specific choice for the threshold value and the lines from bottom to top for each colour shows the values for slow, medium, and fast rotator cases. In both panels, the y-axes are logarithmic for ages up to 1000 Myr and linear at later ages.

4. Stellar EUV and Ly-α emission

A primary motivation for this study is the need to understand the influences of stellar activity on the formation and evolution of planetary atmospheres. In previous sections, we concentrate only on X-ray emission since it has clear and abundant observational constraints, but in most cases the ultraviolet part of the spectrum contributes more to the heating and expansion of planetary upper atmospheres. A complication is that the range of wavelengths of the stellar X-ray and ultraviolet spectrum relevant for atmospheric escape depends on the chemical composition of the atmosphere (described in Sect. 3 of Johnstone et al. 2019a) since different chemical species have different wavelength-dependent absorption cross-sections. For example, much of the heating in the Earth’s thermosphere and upper mesosphere is due to absorption by O₂ and O₃, which are effective at absorbing radiation at wavelengths up to ∼200 nm for O₂ and beyond for O₃. This is in stark contrast to primordial atmospheres composed primarily of hydrogen since H₂ does not absorb radiation longwards of 112 nm. For highly irradiated atmospheres undergoing strong hydrodynamic escape in the form of a transonic Parker wind, such as short-period planets with hydrogen dominated atmospheres, the escape can be driven mostly by X-rays if EUV does not penetrate below the sonic point (Owen & Jackson 2012). In this section, we discuss the evolution of EUV and Ly-α emission. Discussions of emission at far-ultraviolet and longer wavelengths and its relation to X-ray and EUV emission are studied by France et al. (2013, 2018), Linsky et al. (2014), Shkolnik & Barman (2014), and Peacock et al. (2020).

4.1. EUV emission

We define the term ‘XUV’ here to refer to the wavelength range 0.1 and 92 nm We define EUV as the wavelength range 10–92 nm and to be consistent with the literature, we define X-rays to be 2.4–0.1 keV (0.517–12.4 nm). Although there is slight overlap, the additional 2.4 nm contributes only a few percent of the emission at EUV wavelengths. The 2.4 keV boundary is largely arbitrary and as shown in Table 1 of Telleschi et al. (2005), increasing it to 10 keV (0.124 nm) leads to negligible increases in L_X. Despite interstellar absorption making EUV observable for only very nearby stars, some useful observational constraints on EUV emission are available (e.g. Bowyer et al. 2000; Sanz-Forcada et al. 2011; Drake et al. 2020). The EUV emission of stars follow similar trends with rotation and spectral type as X-rays (Mathioudakis et al. 1995), but since more active stars have hotter coronae (Schmitt 1997; Telleschi et al. 2005), a larger fraction of the emitted energy is at shorter wavelengths, meaning that the X-ray to EUV luminosity ratios are higher and EUV emission decays slower with rotational spin-down (Ribas et al. 2005).

An important question is which quantities should be related: possible choices include the luminosities (L_X and L_EUV), the surface fluxes (F_X and F_EUV), and the luminosities normalised to the bolometric luminosities (R_X and R_EUV). We want to know which of these best correlates with the physical nature of the emitting plasma since this determines the spectral energy distribution. This question was addressed by Johnstone & Güdel (2015) who used measurements of the temperatures of coronal plasma for stars with different spectral types and found a single tight mass-independent dependence on X-ray surface flux, given by

$\begin{matrix} {\bar{T}}_{cor} = 0.11 F_{X}^{0.26}, \end{matrix}$ $\begin{aligned} \bar{T}_\mathrm{cor} = 0.11 F_\mathrm{X} ^{0.26}, \end{aligned}$ (18)

where ${\bar{T}}_{cor}$ $\bar{T}_{\mathrm{cor}}$ is the emission measure weighted average temperature of the X-ray emitting plasma. A similar relation was found by Wood et al. (2018) and no such mass-independent relation exists for L_X and R_X. It is therefore reasonable to assume a single relation between X-ray and EUV surface fluxes for all spectral types. Although we do not use Eq. (18) to derive the X-ray–EUV relation, it provides useful physical understanding of the evolutionary trends. Stars with higher X-ray surface fluxes have hotter coronae, and since coronae dominate emission at wavelengths below ∼40 nm (as shown in Fig. 10 of Fontenla et al. 2009), a larger fraction of their XUV emission is at shorter wavelengths. The evolution of coronal temperature is shown in Fig. 13.

Fig. 13.

Evolution of average coronal temperature as a function of stellar mass. Grey circles show our model distribution and blue, green, and red lines show our slow, medium, and fast rotator tracks.

To constrain the X-ray–EUV relation, we compare measurements of X-ray and EUV emission for a sample of nearby stars observed by the Extreme Ultraviolet Explorer (EUVE) spacecraft¹. We include most F, G, K, and M dwarfs in the sample of Craig et al. (1997), including σ² CrB, which is an RS CVn binary composed of main-sequence F and G dwarfs. We include additionally EK Dra, π¹ UMa, and BL/UV Cet. For each star, we collect masses, radii, and interstellar absorption from the literature, as summarised in Table 3. For binaries, we have two options: firstly, when both stars are expected to contribute to the observed emission, we calculate their surface fluxes by summing the surface areas of the two stars and list both components in Table 3, and secondly when one component is expected to dominate, we consider only that component. Due to ISM extinction, the sample cannot be used to reliably estimate the X-ray–EUV relation at wavelengths beyond 36 nm, so we break the EUV into segment 1 with the range 10–36 nm and segment 2 with the range 36–92 nm.

Table 3.

Sample of stars with EUV constraints.

We first correct the EUVE spectra for ISM absorption, which for many of our stars is significant at longer wavelengths within the EUV wavelength range considered. We compute the interstellar absorption (the factor exp(−τ) by which the unabsorbed spectrum is reduced, where τ is the optical depth) in XSPEC (Arnaud 1996). We use the tbabs absorption model for a standard composition of the interstellar gas with average admixtures of dust, based on Wilms et al. (2000). We then integrate the unabsorbed spectrum between 10 and 36 nm to get L_EUV, 1. To derive L_X, we use the results of the ROSAT All-Sky Survey (RASS) since that provides a consistent and complete set of X-ray measurements for all stars and all RASS measurements are taken from the NEXXUS database (Schmitt & Liefke 2004)², except EK Dra which is not included and instead we use L_X determined by Güdel et al. (1997), correcting for the slightly different distance estimate that they used. The RASS L_X values listed in the NEXXUS database are determined assuming a single count rate to energy flux conversion factor, but since this factor should be dependent on the spectral shape, we recalculate L_X using the count rates and hardness-ratios (listed in NEXXUS as ‘HR1’), the distances listed in Table 3, and the hardness-ratio dependent conversion factor given by Fleming et al. (1995) and Schmitt et al. (1995).

The correlation between X-ray and EUV surface fluxes is shown in Fig. 14 and confirms that a single mass-independent scaling law between F_X and F_EUV, 1 is appropriate. Since AU Mic and AT Mic are likely on the pre-main-sequence, Fig. 14 suggests that this is also valid for pre-main-sequence stars. We include also daily average solar spectra over several solar cycles calculated from the Flare Irradiance Spectrum Model (Chamberlin et al. 2007). At medium and high solar activity, the slope of this distribution is consistent with the stellar sample, but there is a turnoff at low activity, which could indicate this relation changes for very low activity stars, though King et al. (2018) did not find such a dramatic change in the F_X–F_EUV, 1 relation based on solar spectra derived from the TIMED/SEE mission. To include the Sun, we split the daily averages into low, medium, and high activity states based on F_X assuming three equal-sized bins spaced between the minimum and maximum values and then calculate the average X-ray and EUV fluxes for each bin. Using these three points in our fit ensures the Sun is weighted higher than other stars but still does not dominate the fit. Using the OLS Bisector method, our fit to all stars in the sample gives

$\begin{matrix} log F_{EUV, 1} = 2.04 + 0.681 log F_{X}, \end{matrix}$ $\begin{aligned} \log F_\mathrm{EUV,1} = 2.04 + 0.681 \log F_\mathrm{X} , \end{aligned}$ (19)

Fig. 14.

Relation between X-ray (0.517–12.4 nm) and EUV (10–36 nm) surface fluxes, F_X and F_EUV, for the sample of stars given in Table 3. The colours and sizes of the circles represent spectral type and mass and the yellow region and the inset shows the Sun at different times during the solar cycle. The open red circle shows the M dwarf GJ 832, with fluxes derived from the stellar parameters and semi-empirical model spectrum of Fontenla et al. (2016).

which gives

$\begin{matrix} \frac{F_{EUV, 1}}{F_{X}} = 110 F_{X}^{- 0.319}, \end{matrix}$ $\begin{aligned} \frac{F_\mathrm{EUV,1} }{F_\mathrm{X} } = 110 F_\mathrm{X} ^{-0.319} , \end{aligned}$ (20)

where both fluxes are given in erg s⁻¹ cm⁻². We find a steeper dependence of EUV emission on X-ray emission that those derived by Chadney et al. (2015) and King et al. (2018) largely because we do not fit our relation to the Sun only. Our dependence is shallower than the relation Sanz-Forcada et al. (2011) because we relate surface fluxes instead of luminosities.

There are several sources of uncertainty in empirical X-ray–EUV scaling relations. Stellar radii are rarely accurately determined causing uncertainties in the surface fluxes. The lack of simultaneous X-ray and EUV measurements, which combined with activity variability adds significant scatter to Fig. 14, is another source of uncertainty and could explain why there is more scatter at the high flux part of Fig. 14 since M dwarfs are typically more variable on short timescales (Sect. 3.2). For example, EQ Peg is the main outlier and it has a large number of L_X determinations listed on the NEXXUS database, most of which are ∼10^28.9 erg s⁻¹ which would put it very close to our best fit line, whereas the value from RASS that we use is significantly lower. Another issue is that M dwarfs with very low surface fluxes are not included in the EUVE sample, likely because they necessarily have low EUV luminosities and were not detectable due to ISM absorption. We instead show in Fig. 14 (though do not include in our fit) the M dwarf GJ 832, which has an XUV spectrum derived from the semi-empirical models of Fontenla et al. (2016) as part of the MUSCLES Treasury Survey. This star sits close to out best fit relation, which is reassuring, but it is more luminous in X-rays than expected given its EUV flux, possibly suggesting a small mass dependence in the slope of the F_X–F_EUV, 1 relation. This is however difficult to test since the only stars with surface fluxes similar to that of the Sun are α Cen and Procyon and the M dwarf with the smallest surface fluxes is Proxima Centauri which lies two orders of magnitude above the Sun in Fig. 14. Our sample does not contain many low activity stars because ISM absorption makes such stars undetectable at EUV wavelengths for all but the nearest stars. As can be seen in Fig. 7 of Schmitt (1997), there are many nearby stars with X-ray surface fluxes similar to that of the Sun (10⁴–10⁵ erg s⁻¹ cm⁻²), including also K and M dwarfs, but these stars were not detected by EUVE and are therefore not included in our sample (except α Cen and Procyon).

To constrain the relation between F_EUV, 1 (10–36 nm) and F_EUV, 2 (36–92 nm), we use the Sun only. As shown in Fig. 14, the X-ray–EUV relation derived from our sample of stars is consistent with the relation that we get considering only the active Sun. We therefore consider only solar values with L_X above 10²⁷ erg s⁻¹ since with this threshold value, we get a F_X–F_EUV, 1 relation from the Sun only that is consistent with Eq. (19). We find

$\begin{matrix} log F_{EUV, 2} = - 0.341 + 0.920 log F_{EUV, 1}, \end{matrix}$ $\begin{aligned} \log F_\mathrm{EUV,2} = -0.341 + 0.920 \log F_\mathrm{EUV,1} , \end{aligned}$ (21)

which gives

$\begin{matrix} \frac{F_{EUV, 2}}{F_{EUV, 1}} = 0.924 F_{EUV, 1}^{- 0.0798}, \end{matrix}$ $\begin{aligned} \frac{F_\mathrm{EUV,2} }{F_\mathrm{EUV,1} } = 0.924 F_\mathrm{EUV,1} ^{-0.0798}, \end{aligned}$ (22)

where both fluxes are given in erg s⁻¹ cm⁻². This is less reliable than the F_X–F_EUV, 1 relation derived above since it is derived considering only the Sun, which varies over only a small fraction of the parameter space. This part of the EUV spectrum is also less important since for the Sun, it contributes typically between 30 and 45% of the total EUV (10–92 nm) emission, and this contribution is likely much lower for very active stars.

4.2. Ly-α emission

The most important feature in a star’s far ultraviolet spectrum is the Ly-α emission line at 121.6 nm, which is formed in the transition region and upper chromosphere (Avrett & Loeser 2008) and as another manifestation of magnetic activity, Ly-α correlates well with emission at shorter wavelengths (Linsky et al. 2014). The Ly-α line has been used to constrain the properties of stellar winds and planetary atmospheres (Ehrenreich et al. 2008; Wood et al. 2014; Kislyakova et al. 2014) and since it often has a luminosity higher than that of the entire X-ray and EUV, it is important to understand its evolution. Although most of the line is absorbed by the ISM, reconstructions of the intrinsic Ly-α line fluxes for a large number of stars are available in the literature (Wood et al. 2005).

The relation between X-ray and Ly-α was studied by Linsky et al. (2013) who showed that the ratio of the Ly-α to X-ray emission has a powerlaw dependence on the X-ray flux at 1 AU (which is proportional to L_X). They found that for F, G and K dwarfs this dependence is very similar with only a small offset for K dwarfs, but for M dwarfs, the dependence is shifted to lower Ly-α to X-ray ratios at each X-ray flux and has a much larger scatter (shown in their Fig. 7). The reason for these differences is that, similar to coronal temperature and EUV emission, Ly-α emission depends not on L_X, but scales with the X-ray surface flux, F_X, as shown by Wood et al. (2005). In Fig. 15, we show the relation between F_Lyα/F_X (or equivalently L_Lyα/L_X) and F_X using the fluxes given in Table 1 of Linsky et al. (2013). For this, we convert their 1 AU fluxes into F_X using stellar radii determined from the spectral types listed in their Table 1 and the data given by Pecaut & Mamajek (2013). When F_X is used instead of the flux at 1 AU, the relation becomes mass-independent and we no longer see the large scatter in the values for M dwarfs. This scatter is a result of M dwarfs having a large range of radii (the surface area of an M0V star is 30 times larger than that of an M9.5V star). The best fit relation between F_Lyα and F_X is

$\begin{matrix} log F_{Ly α} = 3.97 + 0.375 log F_{X}, \end{matrix}$ $\begin{aligned} \log F_{\mathrm{Ly} \alpha } = 3.97 + 0.375 \log F_\mathrm{X} , \end{aligned}$ (23)

Fig. 15.

Ratio of Ly-α to X-ray emission as a function of X-ray surface flux, F_X, for F, G, K, and M stars derived from the data given in Table 1 of Linsky et al. (2013).

which gives

$\begin{matrix} \frac{F_{Ly α}}{F_{X}} = 1.96 \times 10^{4} F_{X}^{- 0.681}, \end{matrix}$ $\begin{aligned} \frac{F_{\mathrm{Ly} \alpha }}{F_\mathrm{X} } = 1.96 \times 10^{4} F_\mathrm{X} ^{-0.681} , \end{aligned}$ (24)

where both fluxes are given in erg s⁻¹ cm⁻².

4.3. EUV and Ly-α evolution

We show in Fig. 16 the evolution of X-ray surface flux, F_EUV/F_X, and F_Lyα/F_X. In the left column of Fig. 16, circle colours represent which of the three parts of the spectrum considered (X-ray, EUV, and Ly-α) is the most luminous for each star in the model distribution and dark outlines indicate that this part of the spectrum is more luminous than the other two combined. We note that the EUV luminosity is always either the first or second most luminous of these three. It is also important that the spread in activity levels shown in Fig. 16 is likely a result of short-term variability, and at each mass and age, we expect stars to spend some time at each location within this spread.

Fig. 16.

Evolution of X-ray surface flux (left column), the EUV to X-ray ratio (middle column), and the Ly-α to X-ray ratio (right column) as functions of stellar mass. The EUV emission is calculated considering both 10–36 and 36–92 nm. In the left column, circle colours shows which out of X-ray (blue), EUV (green), and Ly-α (red) has the highest luminosity, and circles with dark outlines show that this is more luminous than the other two combined. For example, blue circles with white outlines are used when L_X > L_EUV and L_X > L_Lyα but L_X < L_EUV + L_Lyα, and blue circles with dark outlines are used when L_X > L_EUV + L_Lyα. In the middle and right columns, grey circles show our model distribution and blue, green, and red lines show our slow, medium, and fast rotator tracks. The vertical green line in the lower middle panel shows the range of values for the Sun at activity maximum.

At 2 Myr, most stars are more luminous in X-rays than in EUV and Ly-α, though the distribution does contain some stars that are more EUV luminous, and even some low activity outliers that are more Ly-α luminous, especially at lower masses. At later ages, the decline in activity leads to cooler coronal temperatures and a decay in emission that is more rapid in X-ray than in EUV and Ly-α, which leads to an increase in F_EUV/F_X and an even more rapid increase in F_Lyα/F_X. By 5 Gyr, the EUV to X-ray ratios are mostly distributed between 1.5 and 4. This is lower than the ratio for the Sun, which at activity maximum is typically between 4 and 7, possibly because the Sun appears to be less active than other stars with similar masses and rotation rates (Reinhold et al. 2020). At this age, some low activity stars are more than an order of magnitude more luminous in Ly-α than in X-rays.

5. Fluxes in the habitable zone

It is interesting to consider also the mass dependence and evolution of X-ray and EUV fluxes in the habitable zone. We define the habitable zone orbital distance as being half way between the moist and maximum greenhouse limits calculated using the relations derived by Kopparapu et al. (2013). Although stellar properties evolve significantly on the pre-main-sequence, we are interested in planets that spend billions of years in the habitable zone and therefore calculate time-independent habitable zone boundaries assuming stellar properties at 5 Gyr, shown in Fig. 17. While the bolometric luminosity is the main factor determining these boundaries, the effective temperature, T_eff also has an effect. Cooler stars have photospheric spectra that are shifted to longer wavelengths relative to the Sun’s making them more effective at heating the surfaces of planets, meaning that habitable zone planets orbiting stars with lower T_eff values generally receive lower bolometric fluxes from their host stars. This could cause habitable zone planets orbiting lower mass stars to receive a smaller X-ray flux than those orbiting higher mass stars, so long as both stars are in the saturated regime and on the main-sequence, but we find that this effect is not significant.

Fig. 17.

Evolution of stellar bolometric emission as a fraction of the value at 5 Gyr for different stellar masses from the stellar evolution models of Spada et al. (2013). The values on the y-axis can be interpreted both as the normalised bolometric luminosities and as the normalised bolometric fluxes in the habitable zone. The habitable zone orbital distances that we use are based on the 5 Gyr stellar properties and are shown in the inset as the dashed black line.

In Figs. 9 and 11, we show the evolution of the X-ray flux in the habitable zone, F_X, HZ, for stars of different masses and we see that L_X and F_X, HZ have very different mass dependences at all ages. At 2 Myr, habitable zone planets orbiting lower mass stars receive a much higher X-ray flux than those orbiting higher mass stars despite higher mass stars being more luminous. This might be initially surprising since all stars at this age are saturated and the saturation L_X is proportional to the bolometric luminosity, meaning we would expect similar X-ray fluxes in the habitable zone. The reason can be understood if we consider the ratio of the bolometric luminosity at 2 Myr to the value on the main-sequence, remembering that we define the HZ by the main-sequence stellar properties. This is shown in Fig. 17. Due to their slower pre-main-sequence evolution, this ratio is larger for lower mass stars, meaning that HZ planets receive a larger bolometric flux and therefore a larger X-ray flux.

During the first 100 Myr, F_X, HZ decreases by more than an order of magnitude for very low mass stars due to the decreases in their bolometric luminosities. For solar mass stars, a similarly large decrease in F_X, HZ can be seen for slow and medium rotators, but the effect is caused by the decrease in convective turnover times causing these stars to enter the unsaturated regime. Planets orbiting rapidly rotating solar mass stars receive approximately constant X-ray fluxes. The subsequent spin-down at later ages leads to rapid declines in F_X, HZ for almost all stellar masses and this is especially the case for the higher mass stars, leading to a strong relation between stellar mass and F_X, HZ at later ages. At 5 Gyr, HZ planets orbiting low mass M dwarfs receive X-ray fluxes that are two orders of magnitude higher those received by HZ planets orbiting G dwarfs.

In Fig. 12, we show as functions of mass the ages at which the X-ray fluxes in the HZ decay below values of 10, and 100 erg s⁻¹ cm⁻². For comparison, the X-ray flux received by the modern Earth typically varies between 0.15 and 1.15 erg s⁻¹ cm⁻². This is also a measure of the activity lifetimes of stars, though it is very different to the activity lifetimes for L_X discussed in the previous section and these two different ways to measure activity lifetimes have very different mass dependences. For solar mass stars, the HZ X-ray flux decays below 100 erg s⁻¹ cm⁻² at ages of 10 Myr for initilly slow rotators and 700 Myr for initially fast rotators. For lower mass stars, the activity HZ X-ray fluxes remain above this threshold for much longer and 0.1 M_⊙ stars can take up to 5 Gyr to cross the threshold. These timescales are much longer for the 10 erg s⁻¹ cm⁻² threshold and stars with masses below 0.4 M_⊙ will likely never cross the threshold.

We consider also the total XUV energy emitted by the star and the total energy absorbed by HZ planets integrated over evolutionary timescales. These quantities influence the amount of atmospheric gas planets can lose over their lifetimes, though other factors such as the masses of the planets, the compositions of the atmospheres, which escape mechanisms dominate, and the amount of gas available in the atmospheres at different times also influence total losses. In the upper panels of Fig. 18, we show the stellar XUV luminosity integrated between 1 and 1000 Myr and between 1000 and 5000 Myr as functions of stellar mass for the slow, medium, and fast rotator cases. The integrated luminosity is useful for understanding how stars with different masses and initial rotation rates influence planets with the same orbital distances. For stars with masses above ∼0.5 M_⊙, there is a significant difference in the total energy emitted by initially slow and initially fast rotating stars. For solar mass stars, this difference is approximately an order of magnitude in the first 1000 Myr and a factor of two between 1000 and 5000 Myr. In the first 1000 Myr, there is a very strong dependence on stellar mass, with more massive stars emitting significantly more XUV radiation than less massive stars, and the magnitude of this difference depends on the initial rotation rates of the stars. A 0.2 M_⊙ star emits two orders of magnitude less energy than a rapidly rotating solar mass star and one order of magnitude less than a slowly rotating solar mass star. Between 1000 and 5000 Myr, the dependence on initial rotation is much smaller, but the mass dependence is approximately the same.

Fig. 18.

Total stellar XUV (< 100 nm) emission integrated between 1 and 1000 Myr (upper-left) and between 1000 and 5000 Myr (upper-right) and XUV fluence (time integrated XUV flux) in the habitable zone between 1 and 1000 Myr (lower-left) and between 1000 and 5000 Myr (lower-right) as functions of stellar mass for the slow, medium, and fast rotator cases. The background circles show the values for each star in our model distribution.

In the lower panels of Fig. 18, we show the stellar XUV fluence, defined as the time integrated XUV flux, in the habitable zone in these two time periods. The fluence is useful for understanding how stars with different masses and initial rotation rates influence the atmospheres of habitable zone planets. The dependence on stellar mass for the HZ fluence is very different to the dependence for the integrated luminosity. A planet in the habitable zone of a lower mass star receives significantly more XUV energy over both time intervals. In the first 1000 Myr, a HZ planet orbiting a 0.2 M_⊙ star receives two orders of magnitude more energy than a HZ planet orbiting an initially slowly rotating solar mass star but only a factor of a few more than a HZ planet orbiting an initially rapidly rotating solar mass star. Between 1000 and 5000 Myr, this mass dependence is even larger.

6. The contribution of flares

The influence of flares is implicitly considered throughout this study since we do not distinguish between quiescent and flare states in any of the observational constraints on XUV emission. Our description of variability describes how often stars spend at each activity level for each evolutionary stage taking into account both quiescent and flaring emission. In fact, the distinction between quiescent and flaring states is largely arbitrary since much of the quiescent XUV emission is from flares that are not energetic enough to be distinguished from the background emission and it is even possible that all XUV emission from active stars is from flares happening at such high rates that the superposition of each flare lightcurve creates the appearance of steady non-flaring emission (Kashyap et al. 2002; Güdel et al. 2003; Telleschi et al. 2005), which could explain the correlation between stellar activity and coronal temperature (Güdel et al. 2004). This could even be true for inactive stars like the modern Sun since it is possible that the background heating that maintains the large temperatures in the Sun’s outer atmosphere is provided by very small-scale nanoflares (Parker 1988; Jess et al. 2019). The distinction between truly quiescent emission and flares is in this picture not possible although part of the apparently steady radiation could be from non-flaring active regions. XUV flares with a given radiated energy produce higher contrast in stars with smaller surface areas, such as M dwarfs, and in lower activity stars for a given spectral type because less unrelated active region area contributes to the flaring and non-flaring background emission (shown for example by Reale et al. 2004 for an M dwarf and by Telleschi et al. 2005 for solar analogues at different activity levels) biasing the detectability of individual flares in lightcurves.

Although flares do not cause additional XUV emission not considered in preceding sections, it is interesting to consider flare activity more explicitly to understand the evolution and spectral type dependence of flare rates. In recent years, large samples of optical flares have been observed by the Kepler spacecraft (Shibayama et al. 2013; Davenport 2016). These measurements show that rapid rotators flare more often than slow rotators, leading to a decrease in flare rates on the main-sequence over evolutionary timescales as stars spin down (Davenport et al. 2019). They also find that higher mass stars have much more energetic flares: for example, in the catalogue of Kepler flares compiled by Yang & Liu (2019), the upper bound (99th percentile) of flare energies for 0.1–0.2 M_⊙ stars is ∼10³³ erg, whereas this upper bound for solar mass stars is ∼10³⁶ erg, where these values refer to the total energies emitted by the flares in the Kepler bandpass.

The fraction of stars determined to be flaring and the measured optical flare rates were higher for lower mass stars, which has been interpreted to mean M and K dwarfs flare more often than G and F dwarfs. This interpretation is likely incorrect and the higher flaring rates seen on low mass stars is likely caused by detection biases related to the contrast between photospheric and flare emission (Balona 2015). Although flares are ubiquitous among main-sequence F, G, K, and M dwarfs, flares are typically only seen on a few percent of Kepler targets. It is more difficult to detect flares on higher mass stars since the stars have much higher bolometric luminosities and photospheric temperatures closer to the temperatures of flaring plasma in the photospheric and chromospheric footpoints. The lower limit for the energies of flares that can be detected is therefore much higher for higher mass stars, biasing our flare rate statistics. For example, in the flare catalogue of Yang & Liu (2019), the lower bound (1st percentile) for flares on solar mass stars is ∼10³³ erg, which is similar to the upper bound (99th percentile) for flares detected in that catalogue on late M dwarfs and is above the white light energies of even the most energetic solar flares (Woods et al. 2006; Schrijver et al. 2012).

An unbiased comparison of flare rates between different types of stars requires a common lower energy threshold for all stars in a sample above which all flares are counted. Audard et al. (2000) studied the EUV lightcurves and flare statistics of F, G, K, and M dwarfs and found that the frequencies of flares is almost linearly proportional to X-ray luminosity with stars of all spectral types following the same relation. They found the frequencies of flares with total emitted X-ray and EUV energies above 10³² erg to be given by

$\begin{matrix} N (> 10^{32} erg) = 1.9 \times 10^{- 27} L_{X}^{0.95}, \end{matrix}$ $\begin{aligned} N({>}10^{32}~\mathrm{erg} ) = 1.9 \times 10^{-27} L_\mathrm{X} ^{0.95}, \end{aligned}$ (25)

where N is in day⁻¹ and L_X is in erg s⁻¹. Audard et al. (2000) also showed that among the stars considered, flares with energies above 10³² erg typically contribute ∼10% of the energy emitted at X-ray wavelengths. In Fig. 19, we combine Eq. (25) with our estimates for the evolution of L_X described in previous sections to show how the rates of flares depend on stellar mass and age. As expected, long-term activity decay leads to decreases in flare rates at all masses. Since flare rates scale with X-ray luminosity, low mass stars flare less often at energies above 10³² erg than high mass stars and this trend is visible at all ages.

Fig. 19.

Evolution of the rate of flares with total emitted X-ray and EUV energies above 10³² erg for different stellar masses. Red, green and blue lines refer to fast, medium, and slow rotators and grey circles show our model distribution.

The rate of flares with energy E is given by

$\begin{matrix} \frac{d N}{d E} \propto E^{- α}, \end{matrix}$ $\begin{aligned} \frac{\mathrm{d}N}{\mathrm{d}E} \propto E^{-\alpha } , \end{aligned}$ (26)

where estimates of α for solar flares range from 1.35 to 2.90 (Jess et al. 2019). For stellar flares at XUV wavelengths, Audard et al. (2000) found values typically between 1.6 and 2.4 with no clear dependence on spectral type or activity level. The rates of flares above energy thresholds of E₁ and E₂ are related by

$\begin{matrix} \frac{N (> E_{1})}{N (> E_{2})} = {(\frac{E_{1}}{E_{2}})}^{1 - α} . \end{matrix}$ $\begin{aligned} \frac{N\left({>}E_1\right)}{N\left({>}E_2\right)} = \left( \frac{E_1}{E_2} \right)^{1-\alpha } . \end{aligned}$ (27)

We should expect therefore that the rates of flares above any threshold energy follow similar dependences on mass and age as those shown in Fig. 19, though our poor constraints on α make conclusions of this sort uncertain.

To understand the likely effects of flares on habitable zone planets, it is also interesting to consider flares with XUV fluences in the habitable zone above a given threshold. This fluence is the XUV flux from the flare in the habitable zone integrated over the flare duration, given by $E / 4 π a_{HZ}^{2}$ $E/4 \pi a_\mathrm{HZ}^2$ , and we assume here our definition of the habitable zone orbital distance, a_HZ, discussed in Sect. 5. We assume a threshold fluence of 1.8 × 10⁴ erg cm⁻², equivalent to the fluence received by a habitable zone planet orbiting a G dwarf from a flare with an XUV energy of 10³² erg. For each stellar mass, we calculate the corresponding threshold energy needed to give a fluence in the habitable zone equal to this value, which for 0.75, 0.5, and 0.25 M_⊙ stars is 10^31.4, 10^30.7, and 10^30.1 erg respectively. Using Eq. (27), we then calculate the rates of flares above these threshold energies for each stellar mass and age and the results are shown in Fig. 20. Since Eq. (27) requires an assumption for α, we calculate separately cases for values of 1.6 and 2.4 to demonstrate how this influences the results. From X-ray emission of pre-main-sequence stars in the Taurus molecular cloud, Stelzer et al. (2007) found α = 2.4 ± 0.5 and for solar analogues at different activity levels, Telleschi et al. (2005) estimated α values between 2.2 and 2.8, so it is reasonable to expect our α = 2.4 case to be closer to reality for XUV flares, but it is also uncertain how far below 10³² erg we can extrapolate the flare rate distribution without a significant decrease in α.

Fig. 20.

Evolution of the rate of flares with total emitted X-ray and EUV fluences in the habitable zone above 1.8 × 10⁴ erg cm⁻² for different stellar masses. The left and right columns show results for α (Eq. (26)) of 1.6 and 2.4, respectively, showing that our uncertainties in the flare energy distribution can significantly influence our understanding of the effects of flares on habitable zone planets. Red, green and blue lines refer to fast, medium, and slow rotators and grey circles show our model distribution.

As we go to lower mass stars, the threshold energy goes down so we count a larger fraction of flares in our flare rate statistic. This is true for both α values, and in both cases this effect compensates for the fact that lower mass stars flare less often overall. For α = 1.6, these two effects approximately cancel out and we get a mass-independent distribution, meaning that habitable zone planets receive the same amount of XUV energy from flares above our fluence threshold regardless of their host star mass. For α = 2.4, the first effect is stronger and so we get higher flare rates for lower mass stars, meaning that habitable zone planets orbiting lower mass stars are exposed to significantly more flares above our fluence threshold. If the latter case is closer to reality, we can expect that the atmospheres of habitable zone planets are more strongly influenced by the XUV emission of large flares in M dwarf systems than in systems with higher mass stars despite M dwarfs flaring less often.

7. Discussion and conclusions

In this paper, we develop a comprehensive and empirical description of the rotation and XUV evolution of F, G, K, and M dwarfs. Our model is constrained and validated by an extensive catalogue of stellar rotation and XUV emission measurements from the literature, especially in young stellar clusters. The evolution of stellar rotation and XUV emission can be summarised as follows:-

– At ages of ∼1 Myr, the rotation rates of stars are distributed between approximately 1 and 50 Ω_⊙ and this distribution is mass-independent, at least for masses above 0.4 M_⊙.

– This wide distribution gets wider during the pre-main-sequence spin-up phase and then converges on the main-sequence due to wind driven spin-down. Lower mass stars remain rapidly rotating longer and at later ages spin down to slower rotation rates.

– There is a mass-independent relation between stellar Rossby number and R_X described by separate power-laws in the saturated (low Ro) and unsaturated (high Ro) regimes. This relation has a large scatter due likely to stellar variability. At each evolutionary stage, stars likely spend almost 20% of their time with X-ray luminosties that are a factor of three above or below their long-term averages.

– In the first few million years, the high convective turnover times mean that the saturation threshold rotation rates, Ω_sat, are very low and stellar XUV emission depends not on rotation, but on mass only, with higher mass stars being more XUV luminous. During the pre-main-sequence, Ω_sat increases which causes slowly rotating stars to fall out of saturation despite spinning up.

– Lower mass stars have lower Ω_sat and therefore must spin down more before entering the unsaturated regime. Combined with the longer early period of rapid rotation, this means that lower mass stars remain saturated at their peak activity levels for longer. Late M dwarfs remain saturated for billions of years.

– At young ages, the early spread in rotation causes an additional wide spread in XUV emission, which is especially large for higher mass stars such as G and F dwarfs due to their lower Ω_sat. At all masses, initially rapid rotators remain active longer and emit more XUV energy over their lifetimes than initially slow rotators. A star’s mass and initial rotation rate are the two main parameters for determining its XUV evolution.

– Stellar XUV emission decays over evolutionary timescales due to the decreasing L_bol and increasing Ω_sat on the pre-main-sequence and rotational spin-down on the main-sequence. At all ages, higher mass stars tend to be more XUV luminous than lower mass stars.

– At all ages, the XUV fluxes in the habitable zone (assuming HZ boundaries at 5 Gyr) are higher for lower mass stars due to their closer habitable zones and longer evolutionary timescales. Important are both the longer phases of decreasing L_bol on the pre-main-sequence and the longer saturation times on the main-sequence.

– As activity decays, XUV spectra become more shifted to longer wavelengths causing emission at shorter wavelengths to decay more rapidly. For most inactive stars, the Ly-α emission line is more luminous than both X-ray and EUV.

– At all ages, higher mass stars flare more often at all energies than lower mass stars, but habitable zone planets likely receive more XUV energy from flares when orbiting lower mass stars. Flare rates at all energies decrease as activity decays.

As we demonstrate in this paper, a realistic description of the evolution of stellar activity, and especially of X-ray and EUV emission, must be based on a description of rotational evolution and an understanding of how the rotation distribution as a function of stellar mass evolves with time. Single unique decay laws for stellar activity are unable to describe the range of possible evolutionary tracks that stars with different initial rotation rates can follow. Also important is the fact that stars are variable on short timescales and likely spend much of their lives significantly more and less active than the long-term average. We show that it is important to be clear about which measure of activity is being used when describing the mass dependence of activity: for example, XUV luminosity and XUV flux in the habitable zone lead to very different descriptions.

It is commonly believed, especially in the exoplanet community, that M dwarfs are more X-ray and EUV active than G dwarfs, implying that planets orbiting lower mass stars will be subject to higher rates of atmospheric escape. As we show, M dwarfs are in fact less XUV luminous than higher mass stars at all ages and emit much less XUV radiation over their lifetimes. If we consider planets with similar orbital distances, a planet orbiting an M dwarf receives less radiation over its lifetime and therefore likely experiences less overall atmospheric erosion than a similar planet orbiting a G dwarf. However, if we instead consider planets with similar effective temperatures, such as those in the habitable zone, the idea that M dwarfs are more active is justified due to their much longer evolutionary timescales. These conclusions extend to stellar flares: while the rates of flares with a given total energy are higher for G dwarfs at all ages, the rates of flares with a given habitable zone fluence are likely higher for M dwarfs, at least for the most energetic flares.

In this paper, we do not concentrate on evolution beyond the age of the Sun, largely because of this phase is less interesting for planetary atmosphere formation. Until recently, we have lacked observational constraints on the later evolution of rotation and activity and it has usually been assumed that stars continue to follow the expected Skumanich spin-down (P_rot ∝ t^0.5) until the end of the main-sequence. This assumption has been challenged by recent astroseimological age determinations for older field stars with van Saders et al. (2016) finding evidence for older field stars that rotate faster than expected, suggesting that at approximately the age (or Rossby number) of the modern Sun, changes in magnetic activity cause a reduction in the spin-down rate. If correct, this could indicate that stars in the middle of their main-sequence evolution transition into lower activity states (Metcalfe et al. 2016; Metcalfe & Egeland 2019). Evidence for the Sun being less active than expected given its mass and rotation (Reinhold et al. 2020) could suggest that it has already undergone such a transition, and could explain our need for an additional factor in Eq. (4) to reproduce Skumanich spin-down using solar properties.

Fully understanding stellar activity and its evolution is an important step to interpreting the results of recent and upcoming exoplanetary missions. We make available with this study a large number of evolutionary tracks for rotation and XUV emission for the range of stellar masses that we consider which are intended to be used as essential input in studies into the long-term evolution of atmospheres and their surface conditions. In future studies, we will combine our results with state-of-the-art planetary atmosphere evolution models to explore the different ways that stars with different masses and rotation rates influence the atmospheres of terrestrial planets.

¹

EUVE spectra for can be obtained from the MAST archive at https://archive.stsci.edu/euve/ and particularly for the Craig et al. (1997) sample from https://archive.stsci.edu/prepds/atlaseuve/datalist.html

²

RASS X-ray measurements can be obtained from https://hsweb.hs.uni-hamburg.de/projects/nexxus/index.html

Acknowledgments

This study was carried out with the support by the Austrian Science Fund (FWF) project S11601-N16 “Pathways to Habitability: From Disk to Active Stars, Planets and Life” and the related subproject S11604-N16. This work was also scientifically supported by the ExoplANETS-A project (Exoplanet Atmosphere New Emission Transmission Spectra Analysis; http://exoplanet-atmosphere.eu) funded from the EU’s Horizon-2020 programme (grant agreement no. 776403).

References

Affer, L., Micela, G., Favata, F., Flaccomio, E., & Bouvier, J. 2013, MNRAS, 430, 1433 [NASA ADS] [CrossRef] [Google Scholar]
Agüeros, M. A., Bowsher, E. C., Bochanski, J. J., et al. 2018, ApJ, 862, 33 [NASA ADS] [CrossRef] [Google Scholar]
Allain, S. 1998, A&A, 333, 629 [NASA ADS] [Google Scholar]
Amard, L., & Matt, S. P. 2020, ApJ, 889, 108 [Google Scholar]
Argiroffi, C., Caramazza, M., Micela, G., et al. 2016, A&A, 589, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Arnaud, K. 1996, ASP Conf., 17 [Google Scholar]
Audard, M., Güdel, M., Drake, J. J., & Kashyap, V. L. 2000, ApJ, 541, 396 [NASA ADS] [CrossRef] [Google Scholar]
Audard, M., Güdel, M., & Skinner, S. L. 2003, ApJ, 589, 983 [NASA ADS] [CrossRef] [Google Scholar]
Aufdenberg, J. P., Ludwig, H. G., & Kervella, P. 2005, ApJ, 633, 424 [NASA ADS] [CrossRef] [Google Scholar]
Avrett, E. H., & Loeser, R. 2008, ApJS, 175, 229 [Google Scholar]
Balona, L. A. 2015, MNRAS, 447, 2714 [NASA ADS] [CrossRef] [Google Scholar]
Barnes, S. A., Weingrill, J., Fritzewski, D., Strassmeier, K. G., & Platais, I. 2016, ApJ, 823, 16 [NASA ADS] [CrossRef] [Google Scholar]
Belcher, J. W., & MacGregor, K. B. 1976, ApJ, 210, 498 [NASA ADS] [CrossRef] [Google Scholar]
Bessell, M. S. 1991, AJ, 101, 662 [NASA ADS] [CrossRef] [Google Scholar]
Boro Saikia, S., Marvin, C. J., Jeffers, S. V., et al. 2018, A&A, 616, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bouvier, J., Matt, S. P., Mohanty, S., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, T. Henning, et al., 433 [Google Scholar]
Bowyer, S., Drake, J. J., & Vennes, S. 2000, ARA&A, 38, 231 [NASA ADS] [CrossRef] [Google Scholar]
Boyajian, T. S., McAlister, H. A., van Belle, G., et al. 2012, ApJ, 746, 101 [NASA ADS] [CrossRef] [Google Scholar]
Briggs, K. R., Güdel, M., Telleschi, A., et al. 2007, A&A, 468, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Caballero, J. A. 2009, A&A, 507, 251 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Chadney, J. M., Galand, M., Unruh, Y. C., Koskinen, T. T., & Sanz-Forcada, J. 2015, Icarus, 250, 357 [NASA ADS] [CrossRef] [Google Scholar]
Chamberlin, P. C., Woods, T. N., & Eparvier, F. G. 2007, Space Weather, 5, S07005 [Google Scholar]
Claire, M. W., Sheets, J., Cohen, M., et al. 2012, ApJ, 757, 95 [NASA ADS] [CrossRef] [Google Scholar]
Cohen, O. 2017, ApJ, 835, 220 [NASA ADS] [CrossRef] [Google Scholar]
Craig, N., Abbott, M., Finley, D., et al. 1997, ApJS, 113, 131 [NASA ADS] [CrossRef] [Google Scholar]
Cranmer, S. R. 2009, ApJ, 706, 824 [Google Scholar]
Curtis, J. L., Agüeros, M. A., Douglas, S. T., & Meibom, S. 2019, ApJ, 879, 49 [NASA ADS] [CrossRef] [Google Scholar]
Davenport, J. R. A. 2016, ApJ, 829, 23 [Google Scholar]
Davenport, J. R. A., Covey, K. R., Clarke, R. W., et al. 2019, ApJ, 871, 241 [Google Scholar]
Davison, C. L., White, R. J., Henry, T. J., et al. 2015, AJ, 149, 106 [NASA ADS] [CrossRef] [Google Scholar]
Delfosse, X., Forveille, T., Ségransan, D., et al. 2000, A&A, 364, 217 [NASA ADS] [Google Scholar]
Douglas, S. T., Agüeros, M. A., Covey, K. R., & Kraus, A. 2017, ApJ, 842, 83 [NASA ADS] [CrossRef] [Google Scholar]
Douglas, S. T., Curtis, J. L., Agüeros, M. A., et al. 2019, ApJ, 879, 100 [NASA ADS] [CrossRef] [Google Scholar]
Drake, J. J., Kashyap, V. L., Wargelin, B. J., & Wolk, S. J. 2020, ApJ, 893, 137 [Google Scholar]
Eggenberger, P., Miglio, A., Carrier, F., Fernandes, J., & Santos, N. C. 2008, A&A, 482, 631 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ehrenreich, D., Lecavelier Des Etangs, A., Hébrard, G., et al. 2008, A&A, 483, 933 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ekström, S., Meynet, G., Maeder, A., & Barblan, F. 2008, A&A, 478, 467 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Esselstein, R., Aigrain, S., Vanderburg, A., et al. 2018, ApJ, 859, 167 [NASA ADS] [CrossRef] [Google Scholar]
Favata, F., Micela, G., & Reale, F. 2000, A&A, 354, 1021 [NASA ADS] [Google Scholar]
Fernandes, J., Lebreton, Y., Baglin, A., & Morel, P. 1998, A&A, 338, 455 [NASA ADS] [Google Scholar]
Finley, A. J., & Matt, S. P. 2018, ApJ, 854, 78 [NASA ADS] [CrossRef] [Google Scholar]
Fleming, T. A., Molendi, S., Maccacaro, T., & Wolter, A. 1995, ApJS, 99, 701 [NASA ADS] [CrossRef] [Google Scholar]
Fontenla, J. M., Curdt, W., Haberreiter, M., Harder, J., & Tian, H. 2009, ApJ, 707, 482 [Google Scholar]
Fontenla, J. M., Linsky, J. L., Witbrod, J., et al. 2016, ApJ, 830, 154 [NASA ADS] [CrossRef] [Google Scholar]
France, K., Froning, C. S., Linsky, J. L., et al. 2013, ApJ, 763, 149 [NASA ADS] [CrossRef] [Google Scholar]
France, K., Arulanantham, N., Fossati, L., et al. 2018, ApJS, 239, 16 [CrossRef] [Google Scholar]
Gallet, F., & Bouvier, J. 2013, A&A, 556, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gallet, F., & Bouvier, J. 2015, A&A, 577, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Garraffo, C., Drake, J. J., & Cohen, O. 2016, A&A, 595, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gatewood, G., & Han, I. 2006, AJ, 131, 1015 [NASA ADS] [CrossRef] [Google Scholar]
Gonzalez, G. 2016a, MNRAS, 463, 3513 [Google Scholar]
Gonzalez, G. 2016b, MNRAS, 459, 1060 [NASA ADS] [CrossRef] [Google Scholar]
Gray, R. O., Napier, M. G., & Winkler, L. I. 2001, AJ, 121, 2148 [NASA ADS] [CrossRef] [Google Scholar]
Gregory, S. G., Donati, J. F., Morin, J., et al. 2012, ApJ, 755, 97 [NASA ADS] [CrossRef] [Google Scholar]
Gregory, S. G., Adams, F. C., & Davies, C. L. 2016, MNRAS, 457, 3836 [NASA ADS] [CrossRef] [Google Scholar]
Güdel, M., Guinan, E. F., & Skinner, S. L. 1997, ApJ, 483, 947 [NASA ADS] [CrossRef] [Google Scholar]
Güdel, M., Audard, M., Kashyap, V. L., Drake, J. J., & Guinan, E. F. 2003, ApJ, 582, 423 [NASA ADS] [CrossRef] [Google Scholar]
Güdel, M., Audard, M., Reale, F., Skinner, S. L., & Linsky, J. L. 2004, A&A, 416, 713 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Güdel, M., Briggs, K. R., Arzner, K., et al. 2007, A&A, 468, 353 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Guirado, J. C., Martí-Vidal, I., Marcaide, J. M., et al. 2006, A&A, 446, 733 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Guirado, J. C., Marcaide, J. M., Martí-Vidal, I., et al. 2011, A&A, 533, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hartman, J. D., Gaudi, B. S., Pinsonneault, M. H., et al. 2009, ApJ, 691, 342 [NASA ADS] [CrossRef] [Google Scholar]
Hartman, J. D., Bakos, G. Á., Kovács, G., & Noyes, R. W. 2010, MNRAS, 408, 475 [NASA ADS] [CrossRef] [Google Scholar]
Henderson, C. B., & Stassun, K. G. 2012, ApJ, 747, 51 [NASA ADS] [CrossRef] [Google Scholar]
Herbst, W., Bailer-Jones, C. A. L., Mundt, R., Meisenheimer, K., & Wackermann, R. 2002, A&A, 396, 513 [EDP Sciences] [Google Scholar]
Holzwarth, V., & Jardine, M. 2005, A&A, 444, 661 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Houk, N., & Smith-Moore, M. 1988, Michigan Catalogue of Two-dimensional Spectral Types for the HD Stars. Volume 4, Declinations -26.0 to -12.0, 4 [Google Scholar]
Irwin, J., Hodgkin, S., Aigrain, S., et al. 2007, MNRAS, 377, 741 [NASA ADS] [CrossRef] [Google Scholar]
Irwin, J., Hodgkin, S., Aigrain, S., et al. 2008, MNRAS, 383, 1588 [NASA ADS] [CrossRef] [Google Scholar]
Irwin, J., Aigrain, S., Bouvier, J., et al. 2009, MNRAS, 392, 1456 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Irwin, J., Berta, Z. K., Burke, C. J., et al. 2011, ApJ, 727, 56 [Google Scholar]
Isobe, T., Feigelson, E. D., Akritas, M. G., & Babu, G. J. 1990, ApJ, 364, 104 [Google Scholar]
James, D. J., & Jeffries, R. D. 1997, MNRAS, 292, 252 [NASA ADS] [CrossRef] [Google Scholar]
Jardine, M. 2004, A&A, 414, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Jeffries, R. D., Evans, P. A., Pye, J. P., & Briggs, K. R. 2006, MNRAS, 367, 781 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Jeffries, R. D., Jackson, R. J., Briggs, K. R., Evans, P. A., & Pye, J. P. 2011, MNRAS, 411, 2099 [NASA ADS] [CrossRef] [Google Scholar]
Jess, D. B., Dillon, C. J., Kirk, M. S., et al. 2019, ApJ, 871, 133 [Google Scholar]
Johnstone, C. P. 2017, A&A, 598, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Johnstone, C. P. 2020, ApJ, 890, 79 [Google Scholar]
Johnstone, C. P., & Güdel, M. 2015, A&A, 578, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Johnstone, C. P., Jardine, M., Gregory, S. G., Donati, J. F., & Hussain, G. 2014, MNRAS, 437, 3202 [NASA ADS] [CrossRef] [Google Scholar]
Johnstone, C. P., Güdel, M., Stökl, A., et al. 2015a, ApJ, 815, L12 [Google Scholar]
Johnstone, C. P., Güdel, M., Brott, I., & Lüftinger, T. 2015b, A&A, 577, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Johnstone, C. P., Pilat-Lohinger, E., Lüftinger, T., Güdel, M., & Stökl, A. 2019a, A&A, 626, A22 [EDP Sciences] [Google Scholar]
Johnstone, C. P., Khodachenko, M. L., Lüftinger, T., et al. 2019b, A&A, 624, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Joy, A. H., & Abt, H. A. 1974, ApJS, 28, 1 [NASA ADS] [CrossRef] [Google Scholar]
Kashyap, V. L., Drake, J. J., Güdel, M., & Audard, M. 2002, ApJ, 580, 1118 [NASA ADS] [CrossRef] [Google Scholar]
Keenan, P. C., & McNeil, R. C. 1989, ApJS, 71, 245 [NASA ADS] [CrossRef] [Google Scholar]
Kervella, P., Thévenin, F., Morel, P., et al. 2004, A&A, 413, 251 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kervella, P., Bigot, L., Gallenne, A., & Thévenin, F. 2017, A&A, 597, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
King, G. W., Wheatley, P. J., Salz, M., et al. 2018, MNRAS, 478, 1193 [NASA ADS] [Google Scholar]
Kirkpatrick, J. D., Henry, T. J., & McCarthy, D. W. 1991, ApJS, 77, 417 [NASA ADS] [CrossRef] [Google Scholar]
Kislyakova, K. G., Holmström, M., Lammer, H., Odert, P., & Khodachenko, M. L. 2014, Science, 346, 981 [Google Scholar]
Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [Google Scholar]
Kubyshkina, D., Fossati, L., Mustill, A. J., et al. 2019, A&A, 632, A65 [CrossRef] [EDP Sciences] [Google Scholar]
Lépine, S., Hilton, E. J., Mann, A. W., et al. 2013, AJ, 145, 102 [NASA ADS] [CrossRef] [Google Scholar]
Levato, H., & Abt, H. A. 1978, PASP, 90, 429 [NASA ADS] [CrossRef] [Google Scholar]
Linsky, J. L., France, K., & Ayres, T. 2013, ApJ, 766, 69 [NASA ADS] [CrossRef] [Google Scholar]
Linsky, J. L., Fontenla, J., & France, K. 2014, ApJ, 780, 61 [NASA ADS] [CrossRef] [Google Scholar]
Luger, R., Barnes, R., Lopez, E., et al. 2015, Astrobiology, 15, 57 [NASA ADS] [CrossRef] [Google Scholar]
MacGregor, K. B., & Brenner, M. 1991, ApJ, 376, 204 [NASA ADS] [CrossRef] [Google Scholar]
Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars (Springer Berlin Heidelberg) [Google Scholar]
Magaudda, E., Stelzer, B., Covey, K. R., et al. 2020, A&A, 638, A20 [CrossRef] [EDP Sciences] [Google Scholar]
Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [NASA ADS] [CrossRef] [Google Scholar]
Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64 [CrossRef] [Google Scholar]
Marino, A., Micela, G., & Peres, G. 2000, A&A, 353, 177 [NASA ADS] [Google Scholar]
Marino, A., Micela, G., Peres, G., & Sciortino, S. 2002, A&A, 383, 210 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Marino, A., Micela, G., Peres, G., & Sciortino, S. 2003, A&A, 406, 629 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass, G. G., & Worley, C. E. 2001, AJ, 122, 3466 [Google Scholar]
Mathioudakis, M., Fruscione, A., Drake, J. J., et al. 1995, A&A, 300, 775 [NASA ADS] [Google Scholar]
Matt, S., & Pudritz, R. E. 2008, ApJ, 681, 391 [NASA ADS] [CrossRef] [Google Scholar]
Matt, S. P., MacGregor, K. B., Pinsonneault, M. H., & Greene, T. P. 2012, ApJ, 754, L26 [NASA ADS] [CrossRef] [Google Scholar]
Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJ, 799, L23 [NASA ADS] [CrossRef] [Google Scholar]
McQuillan, A., Mazeh, T., & Aigrain, S. 2014, ApJS, 211, 24 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Meibom, S., Mathieu, R. D., & Stassun, K. G. 2009, ApJ, 695, 679 [NASA ADS] [CrossRef] [Google Scholar]
Meibom, S., Barnes, S. A., Platais, I., et al. 2015, Nature, 517, 589 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Messina, S., Millward, M., Buccino, A., et al. 2017, A&A, 600, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Metcalfe, T. S., & Egeland, R. 2019, ApJ, 871, 39 [NASA ADS] [CrossRef] [Google Scholar]
Metcalfe, T. S., Egeland, R., & van Saders, J. 2016, ApJ, 826, L2 [Google Scholar]
Micela, G., Sciortino, S., Harnden, F. R., et al. 1999, A&A, 341, 751 [Google Scholar]
Monsignori Fossi, B. C., Landini, M., Fruscione, A., & Dupuis, J. 1995, ApJ, 449, 376 [NASA ADS] [CrossRef] [Google Scholar]
Montes, D., López-Santiago, J., Gálvez, M. C., et al. 2001, MNRAS, 328, 45 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Moraux, E., Artemenko, S., Bouvier, J., et al. 2013, A&A, 560, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Morin, J., Donati, J. F., Petit, P., et al. 2008, MNRAS, 390, 567 [Google Scholar]
Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23 [Google Scholar]
Newton, E. R., Irwin, J., Charbonneau, D., et al. 2016, ApJ, 821, 93 [NASA ADS] [CrossRef] [Google Scholar]
Newton, E. R., Irwin, J., Charbonneau, D., et al. 2017, ApJ, 834, 85 [NASA ADS] [CrossRef] [Google Scholar]
Nielsen, M. B., Gizon, L., Schunker, H., & Karoff, C. 2013, A&A, 557, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Núñez, A., & Agüeros, M. A. 2016, ApJ, 830, 44 [Google Scholar]
Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931 [Google Scholar]
Pantolmos, G., & Matt, S. P. 2017, ApJ, 849, 83 [Google Scholar]
Parker, E. N. 1988, ApJ, 330, 474 [Google Scholar]
Peacock, S., Barman, T., Shkolnik, E. L., et al. 2020, ApJ, 895, 5 [Google Scholar]
Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9 [Google Scholar]
Petit, P., Donati, J. F., Aurière, M., et al. 2005, MNRAS, 361, 837 [NASA ADS] [CrossRef] [Google Scholar]
Pillitteri, I., Micela, G., Sciortino, S., & Favata, F. 2003, A&A, 399, 919 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pillitteri, I., Micela, G., Damiani, F., & Sciortino, S. 2006, A&A, 450, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Plavchan, P., Werner, M. W., Chen, C. H., et al. 2009, ApJ, 698, 1068 [NASA ADS] [CrossRef] [Google Scholar]
Prosser, C. F., Stauffer, J. R., Caillault, J. P., et al. 1995, AJ, 110, 1229 [Google Scholar]
Prosser, C. F., Randich, S., Stauffer, J. R., Schmitt, J. H. M. M., & Simon, T. 1996, AJ, 112, 1570 [NASA ADS] [CrossRef] [Google Scholar]
Prosser, C. P., Randich, S., & Simon, T. 1998, Astron. Nachr., 319, 215 [NASA ADS] [CrossRef] [Google Scholar]
Raghavan, D., McAlister, H. A., Torres, G., et al. 2009, ApJ, 690, 394 [NASA ADS] [CrossRef] [Google Scholar]
Reale, F., Güdel, M., Peres, G., & Audard, M. 2004, A&A, 416, 733 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rebull, L. M., Wolff, S. C., & Strom, S. E. 2004, AJ, 127, 1029 [NASA ADS] [CrossRef] [Google Scholar]
Rebull, L. M., Stauffer, J. R., Bouvier, J., et al. 2016, AJ, 152, 113 [NASA ADS] [CrossRef] [Google Scholar]
Redfield, S., & Linsky, J. L. 2008, ApJ, 673, 283 [NASA ADS] [CrossRef] [Google Scholar]
Reiners, A., & Basri, G. 2008, ApJ, 684, 1390 [NASA ADS] [CrossRef] [Google Scholar]
Reiners, A., Schüssler, M., & Passegger, V. M. 2014, ApJ, 794, 144 [NASA ADS] [CrossRef] [Google Scholar]
Reinhold, T., Shapiro, A. I., Solanki, S. K., et al. 2020, Science, 368, 518 [Google Scholar]
Réville, V., Brun, A. S., Matt, S. P., Strugarek, A., & Pinto, R. F. 2015, ApJ, 798, 116 [NASA ADS] [CrossRef] [Google Scholar]
Ribas, I., Guinan, E. F., Güdel, M., & Audard, M. 2005, ApJ, 622, 680 [NASA ADS] [CrossRef] [Google Scholar]
Roble, R. G., Ridley, E. C., Richmond, A. D., & Dickinson, R. E. 1988, Geophys. Res. Lett., 15, 1325 [Google Scholar]
Rodríguez-Ledesma, M. V., Mundt, R., & Eislöffel, J. 2009, A&A, 502, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rucinski, S. M., Mewe, R., Kaastra, J. S., Vilhu, O., & White, S. M. 1995, ApJ, 449, 900 [NASA ADS] [CrossRef] [Google Scholar]
Santos, A. R. G., García, R. A., Mathur, S., et al. 2019, ApJS, 244, 21 [Google Scholar]
Sanz-Forcada, J., Brickhouse, N. S., & Dupree, A. K. 2003, ApJS, 145, 147 [NASA ADS] [CrossRef] [Google Scholar]
Sanz-Forcada, J., Micela, G., Ribas, I., et al. 2011, A&A, 532, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schmitt, J. H. M. M. 1997, A&A, 318, 215 [NASA ADS] [Google Scholar]
Schmitt, J. H. M. M., & Liefke, C. 2004, A&A, 417, 651 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schmitt, J. H. M. M., Fleming, T. A., & Giampapa, M. S. 1995, ApJ, 450, 392 [NASA ADS] [CrossRef] [Google Scholar]
Schrijver, C. J., Beer, J., Baltensperger, U., et al. 2012, J. Geophys. Res. (Space Phys.), 117, A08103 [Google Scholar]
See, V., Matt, S. P., Finley, A. J., et al. 2019, ApJ, 886, 120 [Google Scholar]
Shibayama, T., Maehara, H., Notsu, S., et al. 2013, ApJS, 209, 5 [Google Scholar]
Shkolnik, E. L., & Barman, T. S. 2014, AJ, 148, 64 [NASA ADS] [CrossRef] [Google Scholar]
Spada, F., Lanzafame, A. C., Lanza, A. F., Messina, S., & Collier Cameron, A. 2011, MNRAS, 416, 447 [NASA ADS] [Google Scholar]
Spada, F., Demarque, P., Kim, Y.-C., & Sills, A. 2013, ApJ, 776, 87 [Google Scholar]
Stauffer, J. R., Caillault, J. P., Gagne, M., Prosser, C. F., & Hartmann, L. W. 1994, ApJS, 91, 625 [NASA ADS] [CrossRef] [Google Scholar]
Stelzer, B., Neuhäuser, R., & Hambaryan, V. 2000, A&A, 356, 949 [NASA ADS] [Google Scholar]
Stelzer, B., Flaccomio, E., Briggs, K., et al. 2007, A&A, 468, 463 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Stern, R. A., Schmitt, J. H. M. M., & Kahabka, P. T. 1995, ApJ, 448, 683 [Google Scholar]
Strassmeier, K. G. 1994, A&AS, 103, 413 [NASA ADS] [Google Scholar]
Telleschi, A., Güdel, M., Briggs, K., et al. 2005, ApJ, 622, 653 [NASA ADS] [CrossRef] [Google Scholar]
Thévenin, F., Provost, J., Morel, P., et al. 2002, A&A, 392, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tian, F., Toon, O. B., Pavlov, A. A., & De Sterck, H. 2005, ApJ, 621, 1049 [Google Scholar]
Torres, G., & Ribas, I. 2002, ApJ, 567, 1140 [NASA ADS] [CrossRef] [Google Scholar]
Torres, C. A. O., Quast, G. R., da Silva, L., et al. 2006, A&A, 460, 695 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, A&A, 577, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
van Saders, J. L., Ceillier, T., Metcalfe, T. S., et al. 2016, Nature, 529, 181 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014, MNRAS, 441, 2361 [NASA ADS] [CrossRef] [Google Scholar]
Vigan, A., Bonavita, M., Biller, B., et al. 2017, A&A, 603, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Waite, I. A., Marsden, S. C., Carter, B. D., et al. 2017, MNRAS, 465, 2076 [NASA ADS] [CrossRef] [Google Scholar]
Weber, E. J., & Davis, L. 1967, ApJ, 148, 217 [NASA ADS] [CrossRef] [Google Scholar]
West, A. A., Hawley, S. L., Bochanski, J. J., et al. 2008, AJ, 135, 785 [NASA ADS] [CrossRef] [Google Scholar]
Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914 [Google Scholar]
Wood, B. E., Redfield, S., Linsky, J. L., Müller, H.-R., & Zank, G. P. 2005, ApJS, 159, 118 [NASA ADS] [CrossRef] [Google Scholar]
Wood, B. E., Müller, H.-R., Redfield, S., & Edelman, E. 2014, ApJ, 781, L33 [NASA ADS] [CrossRef] [Google Scholar]
Wood, B. E., Laming, J. M., Warren, H. P., & Poppenhaeger, K. 2018, ApJ, 862, 66 [NASA ADS] [CrossRef] [Google Scholar]
Woods, T. N., Kopp, G., & Chamberlin, P. C. 2006, J. Geophys. Res. (Space Phys.), 111, A10S14 [Google Scholar]
Wright, N. J., & Drake, J. J. 2016, Nature, 535, 526 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48 [Google Scholar]
Wright, N. J., Newton, E. R., Williams, P. K. G., Drake, J. J., & Yadav, R. K. 2018, MNRAS, 479, 2351 [NASA ADS] [CrossRef] [Google Scholar]
Yang, H., & Johns-Krull, C. M. 2011, ApJ, 729, 83 [NASA ADS] [CrossRef] [Google Scholar]
Yang, H., & Liu, J. 2019, ApJS, 241, 29 [CrossRef] [Google Scholar]

Appendix A: Grid of rotation and XUV evolution tracks

With this paper, we provide a comprehensive set of evolutionary tracks for stellar rotation and XUV emission. For each case, we provide two files: files labelled ‘basic’ include only the evolution of the surface rotation rate and the X-ray and EUV luminosities, and files labelled ‘extended’ include more detailed information about rotation and XUV emission including envelope and core rotation rates, Rossby number, the various torques calculated in the model, wind mass loss rates, dipole field strengths, and the break-up rotation rates. The XUV emission quantities included in the extended data files are X-ray and EUV luminosities, X-ray surface fluxes, R_X, emission measure weighted average coronal temperature, and the X-ray and EUV fluxes in the habitable zones. We also include X-ray and EUV luminosities at one and two standard deviations of the spread in the Ro–R_X relation above and below the averages, representing the likely variability of stars on non-evolutionary timescales. For each case, we also include the star’s percentile in the rotation distribution at 150 Myr for its mass and the habitable zone boundaries for that mass calculated assuming the 5 Gyr stellar properties.

We provide a basic grid of models for stars with masses between 0.1 and 1.2 M_⊙ and with initial (1 Myr) rotation rates between 0.1 Ω_⊙ and the break-up rotation rate at 1 Myr. Although we there is no observational support for stars at 1 Myr rotating slower than ∼1 Ω_⊙, we include slower rotating tracks for completeness. The tracks start at an age of 1 Myr and end at 12 Gyr for stars with masses below 0.95 M_⊙ and near the end of the main-sequence for higher mass stars. In addition, we include also a second grid of models where the initial rotation rates have been binned instead by the percentile of our model rotation distribution based on the rotation distribution at 150 Myr. We include all integer percentiles between the 2nd and the 98th. As described above, this distribution is based on the rotation distribution of the 150 Myr combined cluster with stars that we find to be rotating unrealistically slow or fast removed. Finally, we also include also tracks for each of the stars in our model distribution and the distributions at each age from 1 Myr to 10 Gyr. We do not include in these distributions stars that have evolved off the main-sequence.

Appendix B: Fitting procedure for rotation model

The rotational evolution model described in Sect. 2.1 contains five unconstrained parameters. These are a_w and b_w from Eq. (7), and a_ce, b_ce, and c_ce from Eq. (12). We constrain these parameters empirically using the observational constraints described in Sect. 2.2 and an MCMC method described here. While our results could reveal useful information about the physical properties of stellar winds and the angular momentum transport in stellar interiors, this is not our aim in this paper.

The goodness-of-fit parameter that we use in our fitting procedure is the likelihood, given by

$\begin{matrix} log L = \frac{1}{N_{obs}} \sum_{i = 1}^{N_{mass}} γ_{i} \sum_{j = 1}^{N_{age, i}} {log}_{10} L_{ij}, \end{matrix}$ $\begin{aligned} \log L = \frac{1}{N_\mathrm{obs} } \sum _{i=1}^{N_\mathrm{mass} } \gamma _i \sum _{j=1}^{N_\mathrm{age,i} } \log _{10} L_{ij} , \end{aligned}$ (B.1)

where N_obs is the total number of observed stars used for all mass bins and ages, N_mass is the number of mass bins, γ_i is a parameter used to adjust the importance of individual mass bins in the fitting procedure, N_age, i is the number of age bins considered in the ith mass bin, and log L_ij is the log (base 10) of the likelihood for the ith mass and jth age bin. This is given by

$\begin{matrix} log L_{ij} = \sum_{k = 1}^{N_{⋆ i j}} log f_{ij} (Ω_{⋆ i j k}), \end{matrix}$ $\begin{aligned} \log L_{ij} = \sum _{k=1}^{N_{\star ij}} \log f_{ij} \left( \Omega _{\star ijk} \right) , \end{aligned}$ (B.2)

where N_⋆ij is the number of stars observed at this mass and age bin, f_ij is the probability density function (PDF) describing how likely it is for a given star at this mass and age to have a rotation rate given by Ω_⋆ijk. The first and second sums in Eq. (B.1) are over all mass and age bins and the sum in Eq. (B.2) is over all observed stars at this mass and age. The values of Ω_⋆ijk are taken from the observational constraints described in Sect. 2.2.

Calculating f_ij for each mass and age requires the initial PDF describing the distribution of rotation rates at 1 Myr and the specification of the fit parameters in our rotational evolution model. Our main assumption is that at each mass and age, the underlying PDF for the rotation rates has a double normal distribution dependence on log₁₀Ω_⋆, which is reasonable given the forms of the observational constraints. For a given mass and age, the PDF is given by

$\begin{matrix} f (Ω_{⋆}) = \sum_{1}^{2} A_{i} \frac{1}{σ_{i} \sqrt{2 π}} exp [\frac{{({log}_{10} Ω_{⋆} - μ_{i})}^{2}}{2 σ_{i}^{2}}] \end{matrix}$ $\begin{aligned} f (\Omega _\star ) = \sum _1^2 A_{i} \frac{1}{\sigma _i \sqrt{2 \pi }} \exp \left[ \frac{\left( \log _{10} \Omega _\star - \upmu _i \right)^2}{2 \sigma _i^2} \right] \end{aligned}$ (B.3)

where μ_i and σ_i are the means and standard deviations of the two normal distributions and A₁ and A₂ controls how important the two normal distributions are relative to each other and are related by A₁ + A₂ = 1. For the starting distribution, we use A₁ = 0.618, μ₁ = 0.714, σ₁ = 0.255, A₂ = 0.382, μ₂ = 1.342, and σ₂ = 0.180. These values are derived from the rotation distributions measured in the Orion Nebula Cluster (Rodríguez-Ledesma et al. 2009) and NGC 6540 (Henderson & Stassun 2012) assuming that Ω_⋆ is given in units of Ω_⊙ = 2.67 × 10⁻⁶ rad s⁻¹. To evolve this function forwards in time for a given stellar mass, we start with a distribution of 500 stars at 1 Myr with rotation rates chosen randomly based on this starting PDF. We then evolve each star forwards in time to the observed ages using our physical model and at each of the observed ages we fit a double normal distribution (the parameters in Eq. (B.3)) using an MCMC method. Doing this separately for each mass bin gives us f_ij(Ω_⋆ijk) from Eq. (B.2) needed to calculate the likelihood.

The method described above allows us to calculate our goodness-of-fit parameter for a given set of our unconstrained model parameters, which we describe collectively as X. Our aim is to find the value of X that corresponds to the maximum value of our goodness-of-fit parameter, which we do using an MCMC method based on the Metropolis-Hastings algorithm. Starting from an initial X chosen randomly between reasonable limits, we iteratively evolve X by making small changes to each of the parameters. The sizes and directions of the changes to the parameters in each step are determined randomly, making our exploration of the parameter space into a random walk. To best explore the parameter space, we perform 1600 separate fit attempts and our final set of parameters are those that give the largest likelihood from all 1600.

Appendix C: Fitting X-ray relation

To use the relation between Ro and R_X given in Eq. (17), we need to empirically determine the parameters β₁, β₂, Ro_sat, and R_X, sat. To fit this relation, we use the OLS(Y|X) method, where Y is log R_X and X is log Ro. This method involves finding the set of parameters that minimised the squares of the vertical distances in Fig. 4 between the observed data points and the power law relation. This is the method recommended by Isobe et al. (1990) when Y has a causal dependence on X or when we want to use the relation to predict Y from a measured X, which most likely best represents our situation. In this case, the cause is the star’s magnetic dynamo, the strength of which is characterised by Ro, and the effect is the resulting magnetic surface field and X-ray emission, measured by R_X. Similarly, our aim in this paper is to use our empirical relation to predict evolutionary tracks for X-ray emission based on observationally constrained tracks for rotation. It therefore seems that OLS(Y|X) is the most appropriate for our purposes, but we cannot rule out that other methods such as the OLS Bisector method could be more appropriate. Given the amount of spread in the relation between Ro and R_X, the OLS Bisector method would lead to a steeper power-law fit in the unsaturated regime (Reiners et al. 2014).

Our fitting procedure involves minimising S, given by

$\begin{matrix} S = \sum_{i} {[log R_{X, i} - log R_{X} (R o_{i})]}^{2}, \end{matrix}$ $\begin{aligned} S = \sum _i \left[ \log R_{\mathrm{X} ,i} - \log R_\mathrm{X} (Ro_i) \right]^2 , \end{aligned}$ (C.1)

where the sum is over all stars in our observed sample, Ro_i and R_X, i are the measured values for the ith star in the sample, and R_X(Ro_i) is the corresponding R_X calculated from Eq. (17) with a Rossby number Ro_i. For a given set of Ro_sat, and R_X, sat, we can easily minimise S separately for all stars with Ro_i < Ro_sat to get β₁ and Ro_i > Ro_sat to get β₂. The only constraint on these fits is the necessity that the lines pass through the saturation point, meaning that the fit must be done numerically, which is trivial given that dS/dβ is a linear function of β and the fit can therefore be performed by finding where dS/dβ = 0. We find the values of Ro_sat and R_X, sat that minimise S using a simple gradient descent method and uncertainties on these parameters are estimated using the bootstrapping method.

All Tables

Table 1.

Sources of the observational constraints on rotational evolution used in this paper.

	Fig. 1. Upper panel: factor f in Eq. (7), as a function of rotation with the blue circles showing the results of MHD simulations for rapidly rotating stars by Johnstone (2017) and the dashed black line showing our fit formula given by Eq. (8). Lower panel: evolution of the break-up rotation rate for different stellar masses.
In the text

	Fig. 10. Age at which a star falls below the saturation threshold as a function of stellar mass for slow (red), medium (green), and fast (blue) rotators. The background circles show when each of the stars in our model distribution become unsaturated. The y-axis is logarithmic for ages up to 1000 Myr and linear at later ages.
In the text

	Fig. 11. Evolutionary tracks for stellar X-ray luminosity for slow, medium, and fast rotators with several masses. The shaded areas around each line represents the spread around the best fit relation, showing one standard deviation of the spread above and below the average.
In the text

	Fig. 13. Evolution of average coronal temperature as a function of stellar mass. Grey circles show our model distribution and blue, green, and red lines show our slow, medium, and fast rotator tracks.
In the text

	Fig. 15. Ratio of Ly-α to X-ray emission as a function of X-ray surface flux, F_X, for F, G, K, and M stars derived from the data given in Table 1 of Linsky et al. (2013).
In the text

	Fig. 19. Evolution of the rate of flares with total emitted X-ray and EUV energies above 10³² erg for different stellar masses. Red, green and blue lines refer to fast, medium, and slow rotators and grey circles show our model distribution.
In the text

The active lives of stars: A complete description of the rotation and XUV evolution of F, G, K, and M dwarfs⋆

1. Introduction

2. Rotational evolution

2.1. Rotational evolution model

2.2. Observational constraints

2.3. Results: The evolution of rotation

3. X-ray evolution

3.1. X-ray relations

3.2. Variability and the nature of the spread

3.3. Observations of X-ray evolution in young clusters

3.4. Results: X-ray evolution

4. Stellar EUV and Ly-α emission

4.1. EUV emission

4.2. Ly-α emission

4.3. EUV and Ly-α evolution

5. Fluxes in the habitable zone

6. The contribution of flares

7. Discussion and conclusions

Acknowledgments

References

Appendix A: Grid of rotation and XUV evolution tracks

Appendix B: Fitting procedure for rotation model

Appendix C: Fitting X-ray relation

All Tables

All Figures

The active lives of stars: A complete description of the rotation and XUV evolution of F, G, K, and M dwarfs^⋆