HAZMAT. VIII. A Spectroscopic Analysis of the Ultraviolet Evolution of K Stars: Additional Evidence for K Dwarf Rotational Stalling in the First Gigayear

The near-UV (NUV, 1700-3200 Å), far-UV (FUV, 912-1700 Å), and extreme-UV (EUV, 100-912 Å) each affect planetary atmospheres differently. The NUV and FUV radiation will photo-dissociate the molecules in the atmosphere, including those most important for the development and detection of life (e.g. H₂O, O₂, OCS, and NH₃) (e.g. Lichtenegger et al., 2010; Segura et al., 2010; Hu et al., 2012). The FUV to NUV flux ratio determines the lifetimes of biogenic species such as CH₄, N₂O, and CH₃Cl, allowing for a greater probability of detection with a greater FUV to NUV ratio (Segura et al., 2005), but will also increase abiotic oxygen and ozone, thus producing higher chances of false-positive biosignatures (e.g. Domagal-Goldman et al., 2014; Tian et al., 2014; Harman et al., 2015). Richey-Yowell et al. (2019) measured this ratio to be up to 10 times smaller for K stars than for M stars, which have been measured photometrically to be between $\approx 0.1-0.3$ (Schneider & Shkolnik, 2018) and spectroscopically to be between $\approx 0.5-3$ (France et al., 2013), giving K stars an added advantage for the detection of real biosignatures. The EUV values are crucial since the EUV will heat and ionize the atmosphere, powering thermal atmospheric escape (Lammer et al., 2003, 2007; Koskinen et al., 2010). However, the EUV flux from stars other than the Sun is predominantly absorbed by the interstellar medium. Limited data exist in the EUV and must therefore be estimated from empirical relationships, differential emission measure models, or stellar atmospheric models constrained by data at FUV, NUV, or X-ray wavelengths (e.g., Pagano, 2009; Linsky & Güdel, 2015; Fontenla et al., 2016; Peacock et al., 2019, 2020; Tilipman et al., 2021).

The photo-dissociation of the molecules in planetary atmospheres can also lead to increased production of hazes in reducing atmospheres of sub-Neptunes and Archean Earth (Earth 4.0 – 2.5 Myr ago) analogs, thus drastically altering the spectrum and the detectability of biosignatures (Zerkle et al., 2012; Arney et al., 2017). For future observations such as with the James Webb Space Telescope (JWST) that will focus on planetary spectra, hazes could limit the detectability of spectral features in the atmosphere (e.g., Fauchez et al., 2019). While no HZ planets around K dwarfs are expected to be seen by TESS due to orbital period constraints (Barclay et al., 2018), $\sim$30 HZ K dwarf terrestrial planets are expected to be discovered in the next decade with ESA’s PLATO mission (PLATO Definition Study Report, 2017).

Several studies have recently focused on K stars as optimal hosts for habitable planets. Both Cuntz & Guinan (2016) and Heller & Armstrong (2014) found that early-type K stars may offer the most suitable conditions for life and may even be “super-habitable,” i.e., even more habitable than an Earth-sized planet around a G2V star. From a photochemistry perspective, studies have demonstrated that K6V stars ($\approx$0.7 M_⊙) may be the most suitable laboratories for producing detectable CH₄ and O₂ (Arney et al., 2017; Arney, 2019), while the increased O₂ of planets around these stars would permit complex life to develop much faster (Lingam & Loeb, 2017).

Using the Galaxy Evolution Explorer (GALEX) broadband UV photometry of K stars ranging in age from 10 Myr to $\approx$5 Gyr, Richey-Yowell et al. (2019) compared the photometric UV evolution of K stars to that of early- and late-type M stars (Shkolnik & Barman, 2014; Schneider & Shkolnik, 2018). They reported that the intrinsic NUV radiation of K stars only decreased slightly (half an order of magnitude) between 10 Myr – $\approx$5 Gyr (the average age for the field stars), and the FUV radiation decreased similarly from 10 – 650 Myr before falling off more steeply (an order of magnitude) by 5 Gyr. This trend is different for late M stars, which are active longer because of their slower spin-down rate, and early M stars, which remain constant through 150 Myr before declining sharply (1.5 orders of magnitude) in both the NUV and FUV. However, GALEX observations only give us a narrow photometric window into the overall evolution of K star UV flux, as some of the strongest and thus most important emission features of these stars are outside of the GALEX bandpasses (e.g. C III, Si III, N V, C II, and Mg II).

In this paper, we present an examination of the evolution of upper atmosphere UV emission lines of 39 K stars, building upon the photometric K star evolution work of Richey-Yowell et al. (2019), as well as expanding the M star UV spectral analysis of Loyd et al. (2021) into the K star regime. Using the Cosmic Origins Spectrograph (COS) on the Hubble Space Telescope (HST), we obtained spectra of K stars at ages 40 Myr (Kraus et al., 2014), 650 Myr (Martín et al., 2018), and $\approx$5 Gyr (the average age of the field stars) in both the NUV (2000-3100 Å) and FUV (1100-1800 Å) in order to measure the UV temporal and rotation evolution of K stars. In §2, we discuss the target selection strategy followed by a description of the methods for analysis in §3. Finally, we present the results in §4 and discussion the implications for these conclusions in §5.

To explore the UV evolution of K stars, we selected stars at three distinct ages: young, intermediate, and old. The young star sample of 10 K stars was taken from the Tucana-Horologium young moving group (40 $\pm$ 5 Myr, Kraus et al. 2014) and confirmed to be bona-fide members through their 3D space velocity motions and lithium depletion by Kraus et al. (2014). Similarly, our intermediate stellar sample of 15 K stars was selected from members of the Hyades moving group (650 $\pm$ 70 Myr, Martín et al. 2018), confirmed through 3D space velocity kinematics by Goldman et al. (2013). Finally, we identified 14 field stars within 30 pc using the SIMBAD database and assumed an age of 5 Gyr for our field age sample, as this is the average age of the stellar field.

K stars with the same spectral type will have different masses at different ages. Therefore, we must consider the evolution of stars of similar mass in time. We estimated the masses of the stars using Table 1 of Richey-Yowell et al. (2019), which were calculated from the spectral type to effective temperature transformations of Pecaut & Mamajek (2013) and the model isochrones of Baraffe et al. (2015). We select stars only from 0.6 – 0.8 M_⊙ to be included in our samples. A summary of the targets is presented in Table 1.3.

The 39 targets were observed in the NUV and/or FUV with the Hubble Space Telescope (HST) Cosmic Origins Spectrograph (COS) instrument. Of these, 16 are new observations from GO-15955 (PI: Richey-Yowell) and 14 are new observations from SNAP-15091 (PI: Agüeros); the remaining are publicly archived observations from GO-14767 (PI: Sing), GO-14784 (PI: Shkolnik), GO-13650 (PI: France), SNAP-14633 (PI: France), and SNAP-11687 (PI: Ayres). The HST program number is identified with each target in Table 2.

For the new targets from GO-15955 presented in this paper, each was visited once by HST and was exposed in the COS G230L (NUV), G160M (FUV), and G130M (FUV) gratings, with central wavelengths 2950 Å, 1577 Å, and 1291 Å  respectively. Total exposure times ranged from $\approx$250 s to $\approx$12,000 s depending on the brightness of the object. Time between exposures were typically $\approx$100 s or $\approx$45 minutes if the Earth was occulting the target during HST’s orbit. The individual exposures were dithered on the detector to reduce noise that may be introduced into the spectra from observations at a single location. For the new targets from SNAP-15091, each was visited once by HST and exposed in the COS G230L grating with central wavelength 2950 Å. Exposure times ranged from $\approx$90 s to $\approx$1,200 s. The additional archived targets were observed in at least one of the COS G230L, G160M, or G130M gratings, although the central wavelengths varied from program to program. A summary of the observations can be found in Table 2.

Initial data reduction was carried out through the CALCOS pipeline version 3.3.10. Since we are interested in the evolution of the quiescent UV flux of these stars, we removed any flares in the data. If a flare was found, we removed the flare by discarding the time region where the flare occurred using the python program costools. The remaining data were then re-run through the CALCOS pipeline to produce flux-calibrated 1D spectra. While these spectra ideally represents a quiescent stellar state, unresolved flares may still be present.

To identify flares, we created light curves by extracting the spectrum (in terms of [counts s^-1 Å^-1]) from the time-tagged (corrtag) files and subtracting the background region (see Figure 3 of Loyd & France 2014). After removing regions contaminated by Earth’s geocoronal emission (such as Ly$\alpha$ and O I), we then applied the flux calibration from the original x1d spectrum to the time-tagged spectrum, to create a conversion from [counts s^-1 Å^-1] to [erg s^-1 cm^-2 Å^-1]. After this, we integrated over the bandpass to create a flux-calibrated light curve for each target binned in 5 second intervals. A dozen flares above 3$\sigma$ were identified and will be analyzed separately in a future paper. No additional flare candidates were identified by eye.

For spectra that were observed during multiple visits, the data were coadded. We used the coaddx1d¹¹1https://github.com/cosmonaut/python-coaddx1d (Danforth et al., 2010, 2016) code originally in IDL and translated into python by Nicholas Nell to coadd the data, weighting by error.

The final spectra were shifted to the rest frame of the star and reviewed individually. Examples representative of the spectra at each age are shown in Figure 1.

Two of the observed targets from GO-15955, STKM 1-495 and UCAC4 518-008871 showed clear binarity in our spectroscopic observations. Additionally, archival Keck NIRC2 images (Proposal I.D. #H246N2, PI: Gaidos) confirm that STKM 1-495 is in fact a binary system. These two targets were thus not considered in our study due to both blending and potential increased activity due to tidal and/or magnetic interactions between their companions.

We analyze the strongest chromospheric and transition region emission features in the COS bandpasses: C III, Si III, N V, C II, Si IV, C IV, He II, and Mg II. These lines were also analyzed in Loyd et al. (2021) and thus we can make direct comparisons with the UV evolution of early M stars. We list the rest wavelength and formation temperature in Table 3, as well as the wavelength range over which we integrate to calculate the surface fluxes.

For more accurate line fluxes, we used the astropy module specutils to fit the continuum immediately outside the spectral regions identified in Table 3 with a cubic Chebyshev polynomial and subtract the fit. The accuracy of these fits were confirmed by eye for each line of each target.

We additionally analyze the NUV continuum as well as the FUV continuum of each spectrum. The continua in these regions are dominated by multiple weak absorption and emission lines and therefore are more accurately considered to be "psuedo-continua". To calculate the pseudo-continua, we consider the same regions as Loyd et al. (2021), which cumulatively cover 251 Å. The wavelengths of the regions over which we integrated are additionally found in Table 3. These regions exclude areas of strong emission lines not considered in this study and non-stellar background contamination such as geo-coronal emission features. For the archived field-age targets which only had data in one of the two FUV gratings, we take the ratio of the target pseudo-continuum available to the same region in HD 266611 (the target with the FUV pseudo-continuum flux closest to the average of the FUV pseudo-continuum flux for the entire sample), and multiply the FUV pseudo-continuum flux of HD 266611 by that ratio. This assumes that the shape of the continuum does not vary significantly between targets, as is the case in Figure 2.

In order to calculate the surface flux from the observed flux, we estimated the radius of the star first using the spectral type to effective temperature relationship in Pecaut & Mamajek (2013) and then using the effective temperature to mass and radius relationships in Baraffe et al. (2015).

While the ages of the moving group members are well constrained, the ages of field stars are typically estimated from isochrones and gyrochrones. We utilized the stardate package by Angus et al. (2019) to determine these ages, providing the $G,BP,$ and $RP$ magnitudes and parallaxes from $Gaia$ EDR3 (Brown et al., 2021), the $J,H,$ and $K$ magnitudes from 2MASS (Skrutskie et al., 2006), and the rotation period if available. A summary of these data are seen in Figure 3.

The rotation periods were gathered from the literature where available. For stars with no reported error in their rotation period, we assume 10%. For those without measured rotation periods, we carried out a Lomb-Scargle (Lomb, 1976; Scargle, 1982) search of $TESS$ lightcurves using the program lightkurve (Lightkurve Collaboration et al., 2018). The lightkurve periodograms test frequencies up to the Nyquist frequency, using grid steps of $1/\textrm{observation time}$. UCAC3 13-31217, HD 266611, HD 85512, BD+57 1274, BD+49 2126, and BD+05 2767 had TESS observations; however, no clear variability was identified, which is not uncommon for older stars.

The Rossby number is defined as the ratio of a star’s rotation period to its convective turnover time (i.e. the time it takes for circulation within a convective cell): $R_{o}=P_{rot}/\tau$. This number was first shown by Noyes et al. (1984) to be correlated to the magnetic activity of stars, and therefore should relate directly to the chromospheric, transition region, and coronal surface fluxes of stars which are driven by this magnetic activity. To calculate the convective turnover time, we use the empirical relationship between $V-K_{s}$ color and the turnover timescale of convection by Wright et al. (2018). For two stars in our data without $V-K_{s}$ colors, we instead use the mass – convective turnover time relationship, again established by Wright et al. (2018). These were found to produce consistent Rossby numbers with the rest of the data set.

The complete list of estimated stellar ages, measured rotation periods, and derived Rossby numbers can be found in Table 1.3.

To fit the surface flux, $F$, as a function of age, rotation period, or Rossby number, we adopt a piecewise model that includes a period of saturated (constant) activity, followed by a power law decline, as utilized in other chromospheric studies (e.g. Wright et al., 2011; Jackson et al., 2012; Loyd et al., 2021; Pineda et al., 2021):

Using GALEX broadband UV photometry of 400 K stars ranging in age from 10 Myr to $\approx$5 Gyr, Richey-Yowell et al. (2019) saw a decrease (half an order of magnitude) in median flux by 650 Myr; however, the interquartiles (the middle 50% of the population) continued to overlap with the younger stars. Figure 7 shows the relative GALEX UV surface fluxes of the stars in this work compared to their complete sample. The stars in this work are typical compared to other K stars and fall within the interquartiles of the GALEX sample.

It should also be noted that the X-ray evolution of K stars (e.g. Richey-Yowell et al., 2019; Jackson et al., 2012) starts to decline after 10 Myr, much earlier than seen in the UV, unlike for M stars where both the UV and X-ray decline at similar ages (Shkolnik & Barman, 2014). In comparing the fractional X-ray (coronal) and Ly$\alpha$ (chromospheric) flux to bolometric luminosity, Linsky et al. (2020) observed that the relative emission from the chromosphere and corona show different distributions between M and K stars, potentially representing different heating mechanisms between the chromospheres and coronae of these stars. The corona can also heat the transition region through thermal conduction in a process known as back-heating (e.g., Kopp & Kuperus, 1968; Athay, 1990). Therefore, the FUV lines in this study that map the transition region may also experience a different evolution than the X-rays or other UV chromospheric lines.

Several studies (e.g. Tu et al., 2015; Gondoin, 2018; Magaudda et al., 2020) have shown that the initial rotation period of low-mass stars affects the evolution of the X-ray flux. Stars with faster initial rotation periods are saturated for longer and decline at a faster rate; conversely, stars with slower initial rotation periods decline in flux at younger ages, but decline at a lower rate. This trend most likely continues into the UV and accounts partially for the large scatter seen in Richey-Yowell et al. (2019) among K dwarfs. However, the rotation rates of our sample span a range of periods, yet we see constant flux for those at shorter or longer periods within the same age group. Therefore, it is unlikely that the prolonged period of saturated activity observed is due to a selection of only stars with faster original rotation periods.

Using K2 data from open clusters, Curtis et al. (2019, 2020) show that current gyrochronological models fail to predict the rotation period of K stars beyond 650 Myr. Instead, the K stars appear to stall their spin-down, potentially maintaining their rotation period for over 1.3 Gyr before resuming to spin down. This effect is not seen in the same rotation data from Curtis et al. (2019) for G stars and it remains unknown if this occurs in M stars.

Spada & Lanzafame (2020) were able to closely recreate these observations using a two-zone model for the star, where differential rotation between the convective envelope and the radiative core develops as magnetic wind braking reaches a maximum. The stalling may occur due to a balance between the spin-down torque from the magnetic wind braking and a spin-up torque from the transfer of angular momentum from the more quickly rotating radiative core. The redistribution of angular momentum is highly mass dependent, therefore the lower mass stars (i.e., K stars) are stalled for longer periods. Curtis et al. (2020), however, do not rule out other theories such as reduced magnetic braking. If the core-envelope recoupling is the underlying method for the stalling, then this would imply that M dwarfs would not experience this effect, since stars become more fully convective with lower mass, and that there is perhaps a dominant relative size of the convection zone at which stalling is more pronounced.

This phenomenon could explain the prolonged saturation period seen in the results presented above. If the K stars remain at a faster rotation period through a couple of Gyrs, than we would expect the evolution with age to follow a similar trend to that of the evolution of rotation period with age, maintaining a faster rotation period and thus a higher flux for longer periods of time. Similarly, flux as a function of rotation period would not be dissimilar from a case in which the spin-down of the star did not stall, which indeed is what we see when comparing the K stars with the M stars.