Accurate Modeling of Lyα Profiles and Their Impact on Photolysis of Terrestrial Planet Atmospheres

Sarah Peacock; Travis S. Barman; Adam C. Schneider; Michaela Leung; Edward W. Schwieterman; Evgenya L. Shkolnik; R. O. Parke Loyd

doi:10.3847/1538-4357/ac77f2

1. Introduction

The H I Lyα (Lyα; λ1215.67 Å) line comprises ∼37%–75% of the total UV flux (1150–3100 Å) from most late-type stars (France et al. 2013), but is difficult to analyze because in >99% of cases it is severely affected by absorption from neutral hydrogen in the interstellar medium (ISM) and if observed from low-earth orbit, contamination from geocoronal airglow. Due to this, model-dependent reconstructions are needed to approximate intrinsic Lyα fluxes. While significant effort goes into these reconstructions, they depend on assumptions about the shape of the line core and the structure of the ISM. Youngblood et al. (2016) estimate that both could independently yield ∼30% inaccuracies, with one of the main sources of uncertainty for Lyα reconstructions being whether or not self-reversal is assumed in the line profiles. In the case of low resolution observations, Youngblood et al. (2022) found that neglecting self-reversal from Lyα reconstructions of G and K dwarfs can result in overestimations in flux of up to 100%, and worse for M dwarfs, up to 180%.

One of the few stars with noncontaminated Lyα observations is the Sun. Early modeling efforts recognized that departures from local-thermodynamic equilibrium (LTE) are relevant because the line forms across a region where collisions transition from being important to less important, thereby allowing the line source function to deviate from the Planck function. This non-LTE effect is seen most dramatically in the form of a self-reversal in the line core. Recent modeling of late-type stars with the PHOENIX (Hauschildt 1993; Hauschildt & Baron 2006; Baron & Hauschildt 2007) atmosphere code has predicted deep self-reversals in the Lyα line cores with up to 2× lower flux than reconstructed profiles of the same stars (Peacock et al. 2019a, 2019b).

Though it cannot be directly observed for nearly all late-type stars, there does exist one complete Lyα spectrum, that of Kapteyn's star, an M1 subdwarf with RV = 254 km s⁻¹ (Guinan et al. 2016; Youngblood et al. 2022). There are two additional stars, Ross 1044 (M0) and Ross 825 (K3) that have the majority of their Lyα flux shifted out of the geocoronal line core and contamination from the ISM because of their exceptionally large RVs (>150 km s⁻¹; Schneider et al. 2019). These observations reveal either no self-reversals or very slight ones in the line cores, indicating that previous PHOENIX models underpredict the Lyα line core in these stars by a factor of ∼2 (Figure 1).

**Figure 1.** Hubble Space Telescope (HST) Space Telescope Imaging Spectrograph (STIS) spectra of Ross 825 (K3; left), Ross 1044 (M0; middle), and Kapteyn's Star (M1; right) plotted as thick black lines reproduced from Schneider et al. (2019) and Youngblood et al. (2022). An altered version of the HST Cosmic Origins Spectrograph (COS) spectrum of Kapteyn's Star from Guinan et al. (2016) is plotted as the gray curve in the third panel. This line profile was produced by mirroring the right-hand side of the stellar component of the Lyα observation. The large RVs of Ross 1044 (−169.6 km s⁻¹), Ross 825 (−340.2 km s⁻¹), and Kapteyn's Star (+245.2 km s⁻¹) allow their stellar Lyα emission lines to be well separated from the contaminating geocoronal airglow emission, with little Lyα flux scattered by the ISM. Model ISM transmittance curves are plotted as black dotted lines (Schneider et al. 2019; Youngblood et al. 2022). The blue profiles are the intrinsic PHOENIX model profiles, which when multiplied by the ISM transmittance curves yield the red profiles. These red profiles should reproduce the solid black curves, but severely underestimate the core flux.
Download figure:
Standard image High-resolution image

Accurately measuring and modeling the Lyα emission lines from low-mass stars is vital for our ability to build predictive stellar spectra. The hydrogen spectrum, and in particular Lyα, plays an important role in the modeling of low-mass stellar atmospheres for two main reasons: (1) the ionization balance of H i/ H ii in the outer layers partly determines the atmospheric structure, and (2) the core of Lyα is very sensitive to these layers (chromosphere and transition region (TR)) where poorly understood nonradiative heating processes affect the atmospheric structure (Short & Doyle 1997). The bulk of the extreme ultraviolet (EUV) flux spanning the Lyman continuum (≤912 Å) is formed in these same outer layers, and therefore, improving the modeling of the Lyα line directly translates to improvements in modeling the EUV spectrum. Further, accurate stellar Lyα flux estimates are required by photochemical models assessing exoplanet atmospheric abundances as Lyα controls the photodissociation of important molecules, including H₂O and CH₄ (e.g., Rugheimer et al. 2015).

For these reasons, Lyα has been extensively studied over the past several decades. Much of the previous observational work has focused on Lyα line wings, which are more generally accessible than the complete intrinsic profile since optically thick hydrogen absorption located in the intervening ISM attenuates photons in the Lyα core. Analyses of observed Lyα wings have informed us that the line width is directly connected to chromospheric heating, effective temperature, surface gravity, and elemental abundance (Ayres 1979; Linsky 1980). Ayres (1979) found that the widths of chromospheric lines are controlled by the stellar temperature distribution more than chromosphere dynamics or magnetic heating. Specific to Lyα, Gayley (1994) determined that the width of this line is largely controlled by the electron density in the chromosphere, where the broad lines of Lyα form. Recent work by Youngblood et al. (2022) analyzed Lyα observations of high RV stars (RV ∼ ±84–245 km s⁻¹) and confirmed that this particular chromospheric line width is correlated with surface gravity and effective temperature in addition to the depth of the self-reversal, which decreases with increasing surface gravity.

Early modeling efforts investigated the formation of Lyα as a function of atmospheric conditions and chromospheric structure. However, the approximations that were typically assumed and the inability to compare to intrinsic Lyα line profiles left an incomplete understanding of the line formation and sets of constraints on the chromospheric structure.

In a series of papers, the Armagh group (Houdebine & Doyle 1994a, 1994b; Houdebine et al. 1995; Houdebine & Stempels 1997) explored the detailed line formation physics of the hydrogen spectrum in an extensive grid of chromospheric models of early M dwarfs. They found that the temperature break between the chromosphere and transition region (T_b) is the main parameter that influences the line width and wings of Lyα, but that this temperature break zone is only important for high pressure atmospheres. The densest models produce broad profiles with no self-reversal and strong wings, but as the transition region pressure decreases, these lower pressure models present with extended wings, and a self-reversal appears and strengthens. They also found that Lyα line profiles and fluxes are not sensitive to changes in the temperature minimum, turbulent velocity, nor rotational broadening.

These initial Armagh group studies used the non-LTE radiative transfer code of Carlsson (1986) to reproduce hydrogen emission profiles of M dwarf stars, including high-resolution Hα and Hβ profiles and the ratio of Lyα to Hα surface fluxes. Their computations used a 16 level hydrogen atom and assumed complete frequency redistribution (CRD). Short & Doyle (1997) refined these studies using the PHOENIX model atmosphere code and included additional physics, specifically adding line blanketing (via the contribution of numerous chromospheric atomic and molecular lines) in the background radiation field. They also treated additional species in non-LTE. In examining the effects of line blanketing on the computation of the hydrogen spectrum, they found that the predicted equivalent width of Lyα can be reduced by ∼0.3 dex and the corresponding line flux increased by a factor of five, concluding that careful treatment of background opacity is important when modeling the hydrogen spectrum.

Falchi & Mauas (1998) used the Pandora program (Avrett & Loeser 1984) to further explore the impact of common approximations assumed in chromospheric models: the omission of line blanketing to the opacity, assuming CRD in the Lyα line, and the assumption that minority species are in LTE. They found that neglecting line blanketing in the background opacities strongly effects the Lyα line center intensity and also that partial redistribution (PRD) effects are very important for computing Lyα in cool stars. When assuming minority species in LTE versus non-LTE, there was noticeable change in continuum emission shortward of 2600 Å, but they did not find any effect on the Lyα line profile. The increased non-LTE species set used for this study included eight levels of H and 91 levels of 11 ionization stages of He, C, Fe, Si, Ca, Na, Al, and Mg. This set was found to be incomplete by Fuhrmeister et al. (2005), who concluded that the non-LTE treatment of C, N, and O has a significant influence on the amplitude of hydrogen emission lines in M dwarfs. Recently, Peacock et al. (2019a) further increased the non-LTE species set to a total of 62 ionization stages for most light elements up to Ni and added PRD to PHOENIX, confirming that PRD effects are important for reproducing the wings of Lyα line profiles of old late-type M stars.

These previous studies have made great strides in understanding the formation of the Lyα line in low-mass stars, but none have been able to draw concrete conclusions about the intrinsic profiles of these stars due to the lack of observational/empirical constraints. Now equipped with new observations of high radial velocity stars, over 70% of the intrinsic Lyα flux is visible. We use these data to further improve our understanding of the microphysics in the upper atmospheres of low-mass stars and, in the process, improve the predictive capabilities of stellar atmosphere models computing EUV spectra. In Section 2 we explain the specifications in our model setup and the stellar observations used as empirical guidance. In Section 3 we describe a variety of analyses conducted to test the sensitivity of the Lyα line core to different components of the model. We present our results in Section 4. In Section 5 we discuss the effects of age and activity on the computation of the Lyα line and how correctly estimating Lyα flux in synthetic stellar spectra impacts the modeling of observable chemistry in terrestrial planet atmospheres. We summarize our findings and describe future work in Section 6.

2. Model Setup

We computed models with the PHOENIX atmosphere code to reproduce the intrinsic Lyα profiles of Kapteyn's Star, Ross 1044, and Ross 825. For each modeled star, we built photospheric structures defined by literature values of effective temperature (T_eff), mass (M_⋆), surface gravity (log(g)), and metallicity ([Fe/H]) (Table 1). To the photosphere, we added ad hoc chromospheric and transition region structures that have linear temperature rises with log(column mass) up to 2 × 10⁵ K. This maximum temperature is above the range at which Lyα forms (2 × 10³–8 × 10⁴ K) and is consistent with previous similar modeling efforts (Peacock et al. 2019b, 2020). The temperature at the top of each chromosphere ranges from 7000 to 8000 K, determined by the point at which hydrogen becomes fully ionized and the atmosphere becomes thermally unstable. The upper atmospheric structures (Figure 2) were adjusted to simultaneously reproduce high-resolution Hubble Space Telescope (HST) Lyα observations and Galaxy Evolution Explorer (GALEX) near-ultraviolet (NUV) photometry of each star (Table 2, Section 2.1).

**Figure 2.** Temperature-column mass structures for models of Ross 825 (red), Ross 1044 (green), and Kapteyn's Star (blue).
Download figure:
Standard image High-resolution image

Table 1. Stellar Properties

Property	Ross 825	References	Ross 1044	References	Kapteyn's Star	References
Spectral Type	K3	1	M0	2	sdM1.0	3
Age (Gyr)	>10	2	>10	2	${11.5}_{-1.5}^{+0.5}$	4
RV (km s⁻¹)	−340.17 ± 0.67	5	−169.55 ± 1.80	5	+244.99 ± 0.18	5
Distance (pc)	98.26 ± 0.18	5	37.44 ± 0.03	5	3.93 ± 0.01	5
T_eff (K)	4680 ± 177	6	3754 ± 95	7	3570 ± 80 (3650)	8
log(g) (cm s⁻²)	4.75	9	4.99	9	4.96 ± 0.13	10
R_⋆ (R_☉)	0.83	11	0.53	11	0.29 ± 0.025	10
M_⋆ (M_☉)	0.55	9	0.3	9	0.28 ± 0.01	10
[Fe/H]	−1.28	12	−1.01 ± 0.21	13	−0.86 ± 0.05	14

Note. The BHAC97 models (Baraffe et al. 1997) were used to obtain the surface gravity and mass for Ross 825 and Ross 1044. We took a rounded average between two rows (corresponding to T_eff = 4377 and 5068 K) to determine the values for Ross 825. We scaled the models by ${R}_{\star }^{2}/{d}^{2}$ using literature values for radius for Ross 825 (0.57 ${}_{0.07}^{+0.15}$ R_☉, Gaia Collaboration et al. 2018) and Ross 1044 (0.38 ± 0.03 R_☉, Newton et al. 2015), but this resulted in uniform offsets between measured and synthetic optical and IR photometry for the stars. In order to align the models with the observations, the radii had to be increased to the values listed in this table.

References. (1) Bidelman (1985); (2) Schneider et al. (2019); (3) Hawley et al. (1996); (4) Kotoneva et al. (2005), Wylie-de Boer et al. 2010; (5) Gaia Collaboration et al. (2021); (6) Stassun et al. (2018); (7) Newton et al. (2015); (8) Anglada-Escude et al. (2014); (9) Baraffe et al. (1997); (10) Ségransan et al. (2003); (11) this work; (12) Schuster et al. (2006); (13) Newton et al. (2014); (14) Woolf & Wallerstein (2005).

Download table as: ASCII Typeset image

Table 2. Observed and Modeled Lyα and GALEX NUV Fluxes

Star	Observed		Modeled
	Lyα	NUV	Lyα	NUV
Ross 825	${2.65}_{-0.76}^{+0.38}$ × 10⁻¹⁵ ^a	(1.37 ± 0.12) × 10⁻¹²	2.79 × 10⁻¹⁵	1.45 × 10⁻¹²
Ross 1044	${2.54}_{-0.11}^{+0.14}$ × 10⁻¹⁴ ^a	(1.47 ± 0.46) × 10⁻¹³	1.59 × 10⁻¹⁴	1.53 × 10⁻¹³
Kapteyn's Star	${2.88}_{-0.08}^{+0.16}$ × 10⁻¹³ ^b	(1.09 ± 0.01) × 10⁻¹²	3.25 × 10⁻¹³	1.3 × 10⁻¹²
	5.27 ± 0.26 × 10⁻¹³ ^c

Notes. Fluxes are in units of (erg s⁻¹ cm⁻²); F_Lyα was calculated from 1214.7 to 1216.7 Å and F_NUV was calculated over the GALEX NUV filter bandpass (1687–3010 Å). All values of F_Lyα are the intrinsic line fluxes without ISM absorption.

^aSchneider et al. (2019)^bYoungblood et al. (2022)^cGuinan et al. (2016)

Download table as: ASCII Typeset image

The Lyα profiles were computed assuming Voigt functions. We also adopted the findings from the aforementioned studies by including line blanketing in the background opacities, computing Lyα with PRD, and including a robust non-LTE species set. For the line blanketing, we included a total of 15,332 bound–free transitions and 233,871 bound–bound transitions. For the non-LTE calculations, we took into account 15,355 levels for 73 specifically considered ionization stages of the most abundant elements in the Sun,⁹ including a 30 level hydrogen atom.

We tested the effects of varying the microturbulent velocity distribution throughout the atmosphere and confirmed the findings of Falchi & Mauas (1998) that changes in this parameter do not affect the Lyα line profile. This is because Lyα is a broad line and other broadening mechanisms are much more important than the Doppler effect. We ultimately used a microturbulent velocity distribution of 2 km s⁻¹ in the photosphere and a slope in the chromosphere and transition region that is a fraction of the sound speed (0.35 × v_sound) in each layer, with a maximum velocity capped at 10 km s⁻¹, as originally done in Fuhrmeister et al. (2005) and continued in our previous modeling of low-mass stars (Peacock et al. 2019a, 2019b, 2020).

With this initial model setup, we were able to reproduce the observed Lyα line widths but not the line core (Figure 1). In Section 3 we detail various analyses performed to quantify their effects on the intensity of the Lyα core flux and ultimately match the complete observed line profile.

2.1. Empirical Guidance

To constrain the temperature structure in the chromosphere and transition region, we guided the models to simultaneously match HST Lyα observations and GALEX NUV photometry.

The HST observations for Ross 825 and Ross 1044 were taken with the G140M grating on the Space Telescope Imaging Spectrograph (STIS; Schneider et al. 2019). Kapteyn's Star has been observed twice, both with the STIS/E140M grating (Youngblood et al. 2022) and with the lower resolution G130M grating on the Cosmic Origins Spectrograph (COS; Guinan et al. 2016). There is a large difference in the computed Lyα fluxes from two these measurements, the lower resolution COS observation yielding a line flux that is 1.85× that of the higher resolution STIS observation. Youngblood et al. (2022) explain that this difference could be astrophysical in nature or a result of a STIS flux calibration issue. Since the source of the flux inconsistency between measurements is not definitive and the main purpose of this work is to accurately model the line shape of Lyα, we choose to match our model to the higher resolution STIS observation. The high RVs of the stars correspond to a Lyα shift of 1.37 Å, 0.68 Å, and 0.99 Å for Ross 825, Ross 1044, and Kapteyn's Star, respectively. These wavelength shifts expose between 63% and 95% of the intrinsic profiles. To compare to these measurements, we convolved our models to the resolution of the observations, accounted for the radial velocity shifts, and multiplied by the ISM transmittance curves shown in Figure 1 (dotted lines).

The GALEX NUV photometry for these stars were not taken contemporaneously with the HST observations; however, the sample stars are old (>10 Gyr), metal deficient ([Fe/H] < −0.8), and optically inactive (Hα equivalent widths <1.0 Å; Reid et al. 1995; Melbourne et al. 2020), suggesting that the stars were likely emitting similar levels of UV flux during both observations (Houdebine & Stempels 1997). Recent work by See et al. (2021) found that photometric variability amplitude and metallicity are positively correlated, meaning metal-poor stars, such as those in our sample, are less magnetically active. We therefore are not concerned about the noncontemporaneity of the UV observations used as empirical guidance for our models.

3. Analysis

As described in the previous section, our initial models produced Lyα profiles that match the wings of the observations well, but underpredicted the core flux by almost a factor of 2. Here we describe a series of analyses conducted to quantify the effect of various model components on the computation of the line profile.

3.1. Effects of Temperature Structure

First, we explored the impact of the simplifying assumption made in our ad hoc temperature structure, where the chromospheric temperature rise increases linearly with log(column mass). Houdebine & Doyle (1994a) found that model atmospheres of low-mass stars with constant temperature gradients in the chromosphere as a function of log(column mass) best reproduce the observed profiles of hydrogen lines. This finding has been confirmed in many studies including those by Andretta et al. (1997), Short & Doyle (1998), Fuhrmeister et al. (2005), and Peacock et al. (2019b) where synthetic spectra constructed with linear chromospheric structures simultaneously reproduce the observed UV continuum and many line profiles of cool stars. Other chromosphere models by, e.g., Fontenla et al. (2016) and Tilipman et al. (2021) have a steep temperature rise in the lower chromosphere followed by a near-constant temperature plateau in the upper chromosphere similar to the solar structure. The temperature plateau results from the balance of radiative losses with nonradiative heating, and is where singly ionized metals are the dominant stages of ionization (Linsky 2017).

We identified the linearity of the chromospheric temperature structure as a potential source of the discrepancy because the core of Lyα forms at temperatures near the upper chromosphere and lower transition region (Sim & Jordan 2005) and so, by increasing the temperature in the upper chromosphere, there could be a change in the computed line profile. For this test, we computed models with identical prescriptions in the photosphere and transition region, maintaining the same temperature minimum and temperature at the base of the transition region, but varying profiles in the chromosphere (top panel, Figure 3). We used the solar chromospheric structure as guidance, initiating a steep temperature gradient in the lower chromosphere that transitions to a shallower gradient (green) or near-constant temperature plateau (blue). In the second panel of Figure 3, we plot radiative quantities for the center wavelength of Lyα. In each case, the source function is increasing with temperature in most of the chromospheric layers and the atmosphere is relatively close to thermal equilibrium. These are the layers over which the wings of Lyα are forming and so the line profiles present with wings in emission (bottom panel). The differences between the models in the lower-to-mid-chromosphere results in increasingly inflated wings as the initial chromospheric temperature rise extends to higher temperatures.

**Figure 3.** Top: varying temperature structures in the chromospheric layers: Our initial setup with a linear rise (red) is modified by adding a steep rise in the lower chromosphere (green, blue). Top middle: the source (S_λ, solid) and Planck (B_λ, dashed) functions for each model at the center wavelength of Lyα are plotted in corresponding colors. Bottom middle: the ratios of non-LTE to LTE number density (departure coefficients) for hydrogen for the n = 1 (asterisk) and n = 2 (solid) states. Adding nonlinearity to the chromospheric temperature rise results in the source function being more closely tied to the Planck function in the midchromosphere and larger deviations toward the upper chromosphere. In all three cases, there is a downturn in both the source function and departure coefficients near the chromosphere–TR boundary. Bottom: normalized Lyα line profiles. Adding a temperature plateau to the chromosphere increases the Lyα line width, but does not impact the line core flux.
Download figure:
Standard image High-resolution image

The Lyα core forms near the chromosphere–TR boundary, where departures from LTE are large (third panel) and the emerging photons are no longer coupled to the local temperature. In all three models, the source function in these layers turns over, decreasing with the rising temperature and resulting in deeply self-reversed line cores. With these results, we conclude that the simplifying assumption of a linear rise with log(column mass) in the chromosphere does not impact the core of Lyα and we maintain that this general temperature prescription produces spectra with the best overall match to both spectroscopic and photometric UV observations.

Next, we examined the temperature structure over the narrow range where the chromosphere and transition region are joined. We investigated this narrow range because it is the specific region where the core of Lyα forms and where the previous models have a sharp downwards turn in the source function. We considered that there may be a numerical issue resulting from the dramatic change in temperature gradient between the chromosphere (∇T_ch ≃ 10⁶ K g⁻¹ cm²) and the transition region (∇T_TR = 10⁸ K g⁻¹ cm²).

For this test, we took a model with our initial setup (linear temperature rises in the chromosphere and TR; top panel, Figure 4, red) and smoothed the layers around the chromosphere–TR boundary, adding curvature over increasingly broader ranges of temperature and column mass (Figure 4, green and blue). The alterations to the temperature structure over these boundary layers do not impact the source function (top-middle panel, Figure 4, solid curves). In each model, the source function traces the Planck function through the chromosphere, but deviates in the upper chromosphere and TR, decreasing with increasing temperature and resulting in the self-reversed core. The increased curvature results in larger departures from non-LTE in the ground state of hydrogen (bottom-middle panel, Figure 4), but minimal change to the n = 2 state. Ultimately, adding this curvature to the chromosphere–TR boundary does not affect the computed Lyα line profile.

In order to see any change in Lyα, we had to greatly extend the curvature, smoothing the layers between 6000 and 20,000 K. In doing so, the entire UV spectrum shortwards of 2500 Å increased in flux by ∼2 orders of magnitude, but the Lyα core flux only showed modest changes in the depth of the central reversal, increasing the line flux by 10%.

3.2. Departure Coefficients

The departure coefficients, or the ratios of non-LTE to LTE number density, are a proxy for the population of each level and are required for the calculation of the emissivity and absorption coefficients. The Lyα line is formed by electron transitions in the hydrogen atom from the level 2 state to ground (n = 1). The Lyman continuum, which largely shapes the EUV spectrum shortwards of 912 Å, results from bound–free transitions from the ground state, while the Balmer series are transitions to the n = 2 state, including the Hα line (n = 3 → 2, 6562.5 Å) and the Balmer continuum (n = 2 → ∞ , 3647–6563 Å).

A commonality observed in all of the models in the previously described temperature structure analyses was that the level 2 state of hydrogen is underpopulated at the base of the transition region. Falchi & Mauas (1998) also found that this second level of hydrogen was underpopulated for M dwarf chromosphere models computed with the Pandora program (Avrett & Loeser 1984), but for CRD models of inactive stars, only. Our stars are similarly inactive; however the underpopulation is still apparent when including PRD for the computation of Lyα.

In our original model setup, the departure coefficients for the n = 1 state of hydrogen slightly decrease at the chromopshere–TR boundary, while the n = 2 state drops by several orders of magnitude (Figure 5, black). We have found that this underpopulation is a direct result of an upward jump in the radiative rates at the onset of the transition region. When manually increasing the ratio for n = 1 over all layers in the upper chromosphere where previous downturns occurred (Figure 5, red), the flux in the Lyman continuum decreases and the self-reversal in the Lyα core deepens. When the same adjustment is made to the n = 2 state (Figure 5, blue), the Lyman continuum is unchanged and the flux in the core of Lyα increases such that the line is nearly in full emission, displaying only a slight self-reversal. Simultaneously increasing the ratios for both n = 1 and 2 yields additive results from adjusting n = 1 or n = 2 alone (Figure 5, purple): the same decreased flux in the Lyman continuum from adjusting just the n = 1 state, and a slightly weaker Lyα profile than that resulting from adjusting just n = 2.

**Figure 5.** Adjusting the ratios of the non-LTE to LTE number density (n_non−LTE/n_LTE) for hydrogen in the n = 1 and n = 2 states in the layers surrounding the chromosphere–TR boundary affect the flux in the Lyman continuum (n = 1 → ∞ ), Lyα (n = 2 → 1), and Hα (n = 3 → 2). Top row: n_non−LTE/n_LTE for H i(n = 1) (*asterisks*) and H i(n = 2) (solid) in the transition region (log(cm) < −5.5) and upper chromosphere (log(cm) > −5.5). These departure coefficients for the initial `PHOENIX` model setup are shown in the top left panel (black). Subsequent panels in the top row have a minimum of n_non−LTE/n_LTE = 3 set for n = 1 (red), n = 2 (blue), and both n = 1 and 2 (purple). Bottom row: output `PHOENIX` spectra corresponding with each departure coefficient adjustment are plotted in matching colors. Bottom left: the Lyman continuum decreases as H i(n = 1) is increased at the chromosphere–TR boundary (red, purple). Bottom middle: the core of Lyα changes in all three scenarios, with the reversal deepening as H i(n = 1) is increased and filling in as H i(n = 2) is increased. Bottom right: Hα is sensitive to changes in both n = 1 and 2, increasing flux in the line core as the values are increased at the boundary.
Download figure:
Standard image High-resolution image

These findings follow as the bound–bound source function can be expressed with departure coefficients as:

$\begin{eqnarray}&&{S}_{\nu }^{l}\approx \frac{{b}_{u}}{{b}_{l}}{B}_{\nu }=\frac{2h{\nu }^{3}}{{c}^{2}}\frac{\psi /\phi }{\frac{{b}_{l}}{{b}_{u}}{e}^{h\nu /{kT}}-\frac{\chi }{\psi }},\end{eqnarray} \tag{ 1 }$

where b_u and b_l are the departure coefficients, n_non−LTE/n_LTE, for the upper and lower level, respectively, as defined in (Rutten 2003). B_ν is the Boltzmann distribution, given in full form on the right-hand side of the equation, where ψ, ϕ, and χ are the profile functions for extinction, emission, and induced emission. Increasing the departure coefficients in the upper level, n = 2, increases the source function and results in emission. Increasing the departure coefficients in the lower level, n = 1, has the inverse effect.

Hα is sensitive to changes to both the n = 1 and 2 departure coefficients, increasing flux in the line core as the ratios are increased in the upper chromosphere. This emission line does not directly correspond to electron transitions with the ground state of hydrogen, but adjusting the population in the n = 1 state affects the availability of electrons for other transitions, which is why there is still a slight change in flux. In these stars, the Balmer continuum is very weak and these adjustments to the departure coefficients change the total flux in those wavelengths by <1%.

In Figure 6, we analyze the sensitivity of the Lyα line to the ratio of n_non−LTE/n_LTE for hydrogen in the n = 2 state. Our original model is plotted in black with asterisks denoting each model layer. At the chromosphere–TR boundary, there is a sharp decrease of three orders of magnitude in n_non−LTE/n_LTE. We have found that the strength of the self-reversal in the core of Lyα is highly sensitive to the minimum set in the layers surrounding the chromosphere–TR boundary. By setting the minimum to 1.0, the sharp decrease between Model Layer 42 in the chromosphere and Model Layer 41 in the transition region reduces to approximately one order of magnitude and the central reversal in the line profile becomes less severe than in the original model (red and black curves). The difference in Lyα line flux between these two models is ∼10%. Increasing the minimum to 5.0 (purple curve), the difference in n_non−LTE/n_LTE in these layers decreases to a factor of 3 and the Lyα line profile is nearly in full emission, with 1.5× the line flux in the original model.

4. Results

As the computed Lyα profiles are highly sensitive to the changes in the departure coefficients of H i, we were able to reproduce the observations by specifying a tailored minimum for the n = 2 state in the layers surrounding the chromosphere–TR boundary for each star. The final computed Lyα and NUV fluxes for each star are listed in Table 2, reproducing the observed quantities. To match these observations, the minimum is set to three for Ross 825 and Kapteyn's Star, and set to 10 for Ross 1044. We compare the models before and after these minima are set in Figure 7. The difference in Lyα flux for each star ranges from 35% to 150%. The difference in Hα flux ranges from 5% to 20%. There is no discernible change to any other part of the spectrum besides these two emission lines. This result contradicts a finding in Houdebine & Doyle (1994a) that Lyman and Balmer series lines can be modeled almost separately; however, it should be noted that by modifying the departure coefficients manually, we are manipulating the occupation numbers and no longer conserving the total number density of hydrogen. This means that the model is no longer physically consistent; however, 99.99% of the total number density of hydrogen is in the n = 1 state (throughout the entire atmosphere and specifically across the chromosphere–TR boundary), so increasing the occupation numbers in the n = 2 state does not significantly affect the total number density of hydrogen.

**Figure 7.** The corrected `PHOENIX` spectra for Ross 825 (top), Ross 1044 (middle), and Kapteyn's Star (bottom) are plotted in red and compared to both the initial `PHOENIX` models (blue) and the Hubble Space Telescope (HST), Galaxy Evolution Explorer (GALEX), and Hα measurements of each star (black). The Hα observation for Ross 1044 comes from the Palomar/Michigan State University (PMSU) survey (Reid et al. 1995) and yields a flat profile. The Kapteyn's Star Hα observation comes from a HARPS observation taken within 6 months of the Lyα profile (Anglada-Escudé et al. 2016). There are no available Hα observations for Ross 825. The final model spectra have good agreement with the measured Lyα profiles, including core flux, and the NUV photometry. When multiplying the final model spectra by the ISM transmittance curves from Figure 1, the Lyα profiles overlay the HST observations. The comparison with initial `PHOENIX` models show that the adjustments made to the H i(n = 2) departure coefficients result in negligible changes to EUV or NUV wavelengths, only increasing flux in the cores of Lyα and Hα.
Download figure:
Standard image High-resolution image

We have found that the cause for the underpopulation of hydrogen in the n = 2 state is related to an upward jump in the radiative rates, indicating an issue with the mean intensity caused by missing or incorrect opacities in the code. While we are using a large non-LTE species set, there are still some residual transitions computed in LTE. Removing these lines from the model calculation does not affect the radiative rates, therefore, this is not the cause and we suspect there may be some other missing UV opacity source.

Earlier modeling efforts utilized a Lyα to Hα flux ratio in low-mass stars as a diagnostic for determining the thickness of the transition region. Doyle et al. (1990) found that the ratio of observed Lyα to Hα fluxes for M dwarfs is close to one, using low resolution International Ultraviolet Explorer (IUE) satellite data to extract the Lyα profiles. While they separated the geocoronal Lyα from stellar contributions, they did not correct for interstellar absorption. As a result, these early modeling efforts (e.g., Houdebine & Doyle 1994a; Short & Doyle 1997) struggled to reproduce the flux ratio, consistently producing models with Lyα to Hα flux ratios between two and five. Now equipped with Lyα observations that reveal the intrinsic line flux to guide our models, we find that the continuum normalized Lyα to Hα flux ratio is between 5 and 20. This discrepancy is likely due to a combination of the initial modeling efforts not accounting for the ISM absorption in addition to our sample stars being very old and optically inactive, resulting in flat Hα profiles.

5. Discussion

5.1. Effects of Stellar Inactivity

All three of our sample stars are old and inactive, with subsolar metallicity. The low metallicity of these stars results in low electron density at chromospheric temperatures resulting in weak Lyα emission compared to more solar-like stars; however, the effects of metallicity have not been extensively studied in chromospheric modeling previously. Because these stars are not representative of the general stellar population, we consider the effects of both activity and [Fe/H] on the Lyα line and the population of hydrogen in the n = 2 state.

As FGKM stars age, they increase in effective temperature, but decrease in radius and gravity. They spin more slowly, have reduced dynamo production of magnetic fields, and become less UV active as their chromospheres shift outward to lower pressures (Peacock et al. 2020). Older stars also have lower metallicity than their younger counterparts because they were born in an environment where less metals were available. In Figure 8, we compare our initial >10 Gyr model from Section 3 (red) to a star with the same photosphere, but more active upper atmosphere (green) and to a young star, defined by lower T_eff, higher gravity, and solar metallicity, plus a more active upper atmosphere (blue). This young star is representative of the initial model star at 10 Myr, obtaining the photospheric parameters from BHAC15 models (T_eff = 4142 K, log(g) = 4.184, [Fe/H] = 0.0; Baraffe et al. 2015) and the upper atmospheric structure from a 10 Myr early M star (Peacock et al. 2020).

**Figure 8.** Parameterizations for stellar inactivity (including metallicity and the onset temperature for the chromospheric temperature rise) do not affect the severity of the decrease in n_non−LTE/n_LTE at the chromosphere–TR boundary for hydrogen in the n = 2 state (middle panel). We compare our initial model (>10 Gyr, [Fe/H] = −1.28, *red*) to a young, high activity model (10 Myr, [Fe/H] = 0.0, blue) and old, high activity model (>10 Gyr, [Fe/H] = −1.28, green). The old, high activity model has the same photospheric parameters as the initial model, but with it is chromosphere initiated at a higher column mass, similar to that of the young model. Right: normalized Lyα line profiles. The corresponding Lyα profiles present with similarly deep central reversals, but increase in line flux with increasing chromospheric activity and metallicity.
Download figure:
Standard image High-resolution image

In both the >10 Gyr model with increased chromospheric activity and the 10 Myr model with solar metallicity, the n = 2 state of hydrogen is similarly underpopulated at the the boundary between the chromosphere and transition region. The sharp decrease between the ratios of non-LTE to LTE number density at these layers is comparable to that of the initial model regardless of the choice of chromospheric parameterization or [Fe/H]. As a result, the Lyα profiles for all three cases present with deep self-reversals in the line core. The total line flux does increase with activity, increasing by 25% between the initial model and the old, high activity model, and by a factor of ∼10 between the initial and 10 Myr models.

While the underpopulation of hydrogen in the n = 2 state is not significantly affected by including the effects of metallicity in the modeling, it is important to include [Fe/H] in these calculations. Keeping all parameters equal except for [Fe/H], there are noticeable changes throughout the computed spectrum. Considering just Lyα, an increase in [Fe/H] by 0.3 increases the Lyα flux by ∼25%, similar to the changes seen by shifting the initial chromospheric temperature rise inwards toward a higher column mass.

5.2. Photochemical Impact on Terrestrial Planetary Atmospheres

Recent work by Teal et al. (2022) quantified the effects of using different quality stellar UV inputs when modeling the photochemical response of exoplanet atmospheres. They found that comparable results are obtained when using either scaling relations or reconstructed proxy spectra as UV input for modeling modern Earth-like planets, but that observable differences occur in the modeled transmission spectra of hazy planets depending on the choice of UV input. They concluded that the biggest effect comes from the UV continuum. Here we take a closer look specifically at the impact of Lyα on the photochemical response of terrestrial planet atmospheres, and quantify the relative importance of modeling this line correctly.

For this analysis, we used the atmos photochemical model, previously described in Arney et al. (2016), Lincowski et al. (2018), and Teal et al. (2022), to examine the impact of different fills in the core of Lyα (similar to those shown in Figure 6). Atmos is a 1D photochemical model that has been used for a variety of terrestrial exoplanet applications, including to simulate planetary atmospheres similar to those of the modern and early Earth (Arney et al. 2016; Teal et al. 2022). We used the most recent publicly available version of atmos,¹⁰ including all of the modifications described in Teal et al. (2022). These modifications include a high-resolution internal grid that prevents Lyα from being spread over too wide a wavelength range. This is particularly important for our analysis in which we are fine tuning the Lyα profile.

We explored two photochemical reaction templates: one that replicates the modern O₂-rich Earth around the Sun (21% O₂), and another for an anoxic low-methane atmosphere. The "anoxic" atmosphere is not entirely without oxygen; however, its O₂ content is determined self-consistently by photochemical production and destruction, and loss to deposition at the surface, rather than a large fixed mixing ratio. The oxic atmosphere contains 239 reactions and 50 gaseous species while the anoxic case has 392 reactions and 74 gaseous species. The oxic atmosphere case assumes a 288 K surface with the empirical temperature-pressure profile of the modern Earth and water vapor consistent with a surface relative humidity of 80%. Fixed surface flux boundary conditions are assumed for CH₄ and N₂O consistent with the biological productivity of the modern Earth. The low-methane, anoxic atmosphere assumes a 275 K surface with a 180 K stratosphere and the same surface relative humidity of 80%. The fixed surface flux of CH₄ is based on a reasonable maximum abiotic flux from volcanism and water–rock reactions (Guzmán-Marmolejo et al. 2013). The oxic and anoxic cases also vary in their assumed CO₂ concentrations, 400 ppm and 1000 ppm, respectively. In both cases the surface pressure is 1 bar and N₂ is used as a filler gas.

Figure 9 shows the differences in an oxic atmosphere exposed to a stellar spectrum with a deeply self-reversed Lyα profile (Figure 6, black) and the same spectrum with Lyα in full emission (Figure 6, purple). There are minor changes in the overall mixing ratio profiles, specifically located at the top of the atmosphere. In general, the full emission scenario shows slightly lower concentrations for most gases in the mesosphere (>50 km) due to more robust photolysis (with the exception of photolysis products such as CO and O, which are enhanced for the same reason). Gas concentrations in the stratosphere (14–50 km) and troposphere (surface to 14 km) are essentially unchanged. For all atmospheric constituents, the change is less than a factor of two, an amount that does not significantly impact the detectability of spectral features. This is because the locations in the atmosphere with the greatest divergence are of the lowest densities and would contribute minimal opacity in transit transmission observations.

The anoxic case, shown in Figure 10, displays similar trends when comparing the deep self-reversal and full emission Lyα line scenarios. The larger CO₂ concentration for this case results in more far-ultraviolet (and specifically Lyα) shielding at the top of the atmosphere, muting the photolysis impact for most altitudes. The most notable difference between the deep self-reversal and full emission scenarios is in the mesospheric CH₄ concentrations. This is due to more robust photolysis at the top of the atmosphere for the full emission scenario, where there is a higher total amount of stellar UV flux. Similar to the oxic case, these changes are less than a factor of two and would not impact the detectability of spectral features in the planetary spectrum due to the low densities at the most impacted altitudes.

When stellar Lyα is in full emission, the atmospheric profiles for both the anoxic and oxic planets show very slight amounts of increased photolysis at altitudes near the top of the atmosphere. Overall the atmospheric mixing ratios and potentially detectable spectral features are not significantly impacted by the differing Lyα profiles explored here. These findings are consistent with those of Teal et al. (2022) who find that their results are most sensitive to the UV continuum, which is fixed between the scenarios presented here.

6. Summary/Conclusions

The Lyα observations of Ross 825, Ross 1044, and Kapteyn's Star have allowed for mostly unobstructed access to the line cores of Lyα, revealing the intrinsic profiles of old, inactive, low-mass stars. Our modeled profiles show more severe self-reversals with earlier spectral types, consistent with previous findings that the Lyα core reversal depth correlates with surface gravity (e.g., Youngblood et al. 2022).

Modeling the Lyα lines of these stars has demonstrated the sensitivity of the core flux to the departures from LTE in the n = 2 state of H i at the boundary between the chromosphere and transition region. These departure coefficients have revealed that there are missing or incorrect opacities in the code and require different minimums to be set in these layers in order to match each Lyα profile. When the models are adjusted to match the observations, there is an increase in Lyα flux of 35%–150% and minimal effect on the Balmer series or EUV spectrum. The Hα flux decreases by ≤20%, while the Balmer continuum and EUV spectrum both decrease by <1%.

We analyzed the impact of these stellar spectral differences when modeling the photochemical response of high molecular weight terrestrial exoplanet atmospheres. We found that corresponding changes to spectrally active gases are very small compared to other uncertainties and will not impact the modeling of planetary spectra.

We also confirmed that using a simplifying assumption of a linear temperature rise with log(column mass) when modeling the stellar chromosphere produces spectra with the best overall match to UV observations. We found that smoothing the temperature structure at the boundary of the chromosphere and transition region yielded negligible changes to the computed UV spectrum.

6.1. Future Work

The results from this work reemphasize that the depth of the self-reversal in the core of Lyα increases with spectral type; however, it is unclear how exactly the line profile changes with uniform steps in T_eff of less than 1000 K. This is important to know for improving the accuracy of modeling stellar atmospheres for a broad range of stars. In order to match the observed profiles of the three target stars in this work, tailored adjustments had to be made. This limited sample size did not reveal a clear trend for what the minimum value for n_non−LTE/n_LTE needs to be based on the T_eff or gravity of the star, although there may be a correlation with the location of the TR. Increasing this sample size would potentially reveal a correlation between the n = 2 departure coefficients and T_eff or other stellar parameters that can be used to apply an accurate correction to all models.

With the upcoming HST program 16646, we are tripling the observational sample of low-mass stars for which Lyα can be measured directly (RV > 100 km s⁻¹), uniformly sampling stars with T_eff from 3400 to 5500 K in steps of ∼500 K. This will enable us to better understand how the Lyα profile changes with spectral type. Ayres (1979), Linsky (1980), and Youngblood et al. (2022) have each connected chromospheric emission line widths as a function of chromospheric heating, T_eff, surface gravity, and elemental abundance. With these new observations we will determine if the total Lyα line strength is connected as well.

Reconstructed Lyα fluxes have been used extensively to produce correlations with other spectral emission lines (Linsky et al. 2013; Youngblood et al. 2017; Melbourne et al. 2020), broadband UV photometry (Shkolnik & Barman 2014), and X-ray fluxes (Linsky et al. 2020). Such correlations have been used to evaluate the life-supporting capabilities of different types of stars (Cuntz & Guinan 2016; Teal et al. 2022), but may have systematic offsets due to the assumptions made in the Lyα reconstructions. With a better understanding of what the intrinsic Lyα fluxes are for main-sequence stars, we can assess the accuracy of these correlations and subsequent analyses conducted with these inputs.

This research was supported by an appointment to the NASA Postdoctoral Program at the NASA Goddard Space Flight Center, administered by Universities Space Research Association through a contract with NASA. The material includes work performed as part of the CHAMPs (Consortium on Habitability and Atmospheres of M dwarf Planets) team, supported by the National Aeronautics and Space Administration (NASA) under grant No. 80NSSC21K0905 issued through the Interdisciplinary Consortia for Astrobiology Research (ICAR) program. M.L. and E.W.S. acknowledge additional support from the Alternative Earths Team, supported by NASA under grant No. 80NSSC21K0594 issued through the ICAR program. This work made use of tools developed by the Virtual Planetary Laboratory, which is a member of the NASA Nexus for Exoplanet System Science and funded via NASA Astrobiology Program grant No. 80NSSC18K0829.

Appendix

As a supplement to Figures 9 and 10 illustrating the photochemical impact of differing stellar Lyα fluxes on oxic and anoxic terrestrial exoplanet atmospheres, we provide the corresponding altitude-dependent photolysis rates in Figures 11 and 12. Photolysis rates are maximized at the top of the atmosphere (mesosphere) and are greatest for the full emission case (Figure 6, purple profile), followed by the fill 2 (Figure 6, green profile), fill 1 (Figure 6, red profile), and deep self-reversal (Figure 6, black profile) cases. The photolysis rates are essentially identical in the stratosphere and troposphere for all species in both atmospheric scenarios. The only exception is the O₃ photolysis rate for the deep self-reversal case in the oxic atmosphere, though this has a minimal effect on the O₃ profile as shown in Figure 9.

**Figure 11.** Comparison of photolysis rates in an oxic atmosphere for the Lyα fill cases considered in this work. In general, photolysis rates are maximized at the top of the atmosphere and are greatest for the full emission scenario, followed by the fill 1, fill 2, and deep self-reversal cases. A linear x-axis is chosen to best illustrate the differences between the various Lyα fill cases at the top of the atmosphere. These results correspond to the atmospheric mixing ratios shown in Figure 9.
Download figure:
Standard image High-resolution image

**Figure 12.** Same as Figure 11 but for our anoxic atmosphere scenario. These results correspond to the atmospheric mixing ratios shown in Figure 10.
Download figure:
Standard image High-resolution image

Accurate Modeling of Lyα Profiles and Their Impact on Photolysis of Terrestrial Planet Atmospheres

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction