Accurate Modeling of Lyman-alpha Profiles and their Impact on Photolysis of Terrestrial Planet Atmospheres

Accurately measuring and modeling the Lyman-$\alpha$ (Ly$\alpha$; $\lambda$1215.67 \AA) emission line from low mass stars is vital for our ability to build predictive high energy stellar spectra, yet interstellar medium (ISM) absorption of this line typically prevents model-measurement comparisons. Ly$\alpha$ also controls the photodissociation of important molecules, like water and methane, in exoplanet atmospheres such that any photochemical models assessing potential biosignatures or atmospheric abundances require accurate Ly$\alpha$ host star flux estimates. Recent observations of three early M and K stars (K3, M0, M1) with exceptionally high radial velocities (>100 km s$^{-1}$) reveal the intrinsic profiles of these types of stars as most of their Ly$\alpha$ flux is shifted away from the geocoronal line core and contamination from the ISM. These observations indicate that previous stellar spectra computed with the PHOENIX atmosphere code have underpredicted the core of Ly$\alpha$ in these types of stars. With these observations, we have been able to better understand the microphysics in the upper atmosphere and improve the predictive capabilities of the PHOENIX atmosphere code. Since these wavelengths drive the photolysis of key molecular species, we also present results analyzing the impact of the resulting changes to the synthetic stellar spectra on observable chemistry in terrestrial planet atmospheres.


INTRODUCTION
The H I Lyman-α (Lyα; λ1215.67 Å) line comprises ∼37%-75% of the total UV flux (1150-3100 Å) from most late-type stars ), 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 non-contaminated 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 selfreversal in the line core. Recent modeling of late-type stars with the PHOENIX (Hauschildt 1993;Hauschildt & Baron 2006; Baron & Hauschildt 2007) atmosphere code have 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,b).
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 −1 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 −1 ) (Schneider et al. 2019). These observations reveal either no self-reversals or very slight ones in the line cores, indicating that previous PHOENIX models under predict the Lyα line core in these stars by a factor of ∼2 (Figure 1).
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 mod-eling 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 non-radiative 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 2 O and CH 4 (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 −1 ) 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,b;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  Youngblood et al. (2022). An altered version of the HST COS spectrum of Kapteyn's Star from Guinan et al. (2016) is plotted as the grey 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 −1 ), Ross 825 (-340.2 km s −1 ), and Kapteyn's Star (+245.2 km s −1 ) 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. 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 re-duced 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 8 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 set up 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.

MODEL SET UP
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 5 K. This maximum temperature is above the range at which Lyα forms (2 × 10 3 -8 × 10 4 K) and is consistent with previous similar modeling efforts (Peacock et al. 2019b(Peacock et al. , 2020. The temperature at the top of each chromosphere ranges from 7000 -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) NUV photometry of each star (Table 2, Section 2.1).
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 in- cluding 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 1 , 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 −1 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 −1 , as originally done in Fuhrmeister et al. (2005) and continued in our previous modeling of low mass stars (Peacock et al. 2019a(Peacock et al. ,b, 2020. With this initial model set up, 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.  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.

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  and with the lower resolution G130M grating on the Cosmic Origins Spectrograph (COS) ).
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 -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 non-contemporaneity of the UV observations used as empirical guidance for our models.

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:

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 non-radiative heating, and is where singly  (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. Bottommiddle: 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 mid-chromosphere and larger deviations towards 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.
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.
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 6 K g −1 cm 2 ) and the transition region (∇T TR = 10 8 g −1 cm 2 ).
For this test, we took a model with our initial set-up (linear temperature rises in the chromosphere and TR; Figure 4. Top: A zoomed in view of models with the same general temperature structure with linear temperature rises in the chromosphere and TR, but with increasing degrees of curvature at the chromosphere-TR boundary. Our initial model set-up is shown in red. Top-middle: Corresponding: source (S λ , solid ) and Planck (B λ , dashed ) functions at the central wavelength of Lyα and Bottom-middle: departure coefficients for the n=1 (asterisk ) and n=2 (solid ) states of hydrogen. Adding increasing curvature to the temperature structure at the boundary region reduces the severity of change between the temperature gradients in the chromosphere and transition region. This change results in an increase the n=1 departure coefficients over the boundary layers, but does not significantly impact the source function nor the n=2 departure coefficients. Bottom: Adding curvature to the boundary between the chromosphere and transition region does not affect the computed Lyα profiles. 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 6,000 and 20,000 K. In doing so, the entire UV spectrum short-ward 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%.

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α.

Transition Region
Chromosphere Photosphere In our original model set-up, 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.
These findings follow as the bound-bound source function can be expressed with departure coefficients as: where b u and b l are the departure coefficients, n non−LT E /n LT E , 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 effects 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−LT E /n LT E 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−LT E /n LT E . 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−LT E /n LT E 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.

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 surround- 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 HST, 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α.
ing 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 3 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 -150%. The difference in Hα flux ranges from 5 -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 modelled 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 occupa- tion numbers in the n=2 state does not significantly affect the total number density of hydrogen.
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 modelling 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 1, using low resolution International Ultraviolet Explorer (IUE) satellite data to extract the Lyα profiles. While they separated the geocoronal 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 2 -5. 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 -20. This discrepancy is likely due to a combination of the initial modelling efforts not accounting for the ISM absorption in addition to our sample stars being very old and optically inactive, resulting in flat Hα profiles.

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 modelling 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 same photosphere, but more active upper atmosphere (green) and 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).
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 a 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 modelling, 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 towards higher column mass.

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), andTeal 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 1-D 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 Figure 9. Comparison of gas mixing ratio profiles in an oxic atmosphere for the endmember stellar spectra explored here (full emission and deep self-reversal Lyα fill cases). At the top of the atmosphere, some species including CH4, H2O, and N2O show slightly diminished concentrations in the full emission case as a result of increased photolysis. Photolysis products CO, O, and OH, show slightly increased concentrations for the same reason. Below ∼ 55 km, the profiles are essentially identical. the most recent publicly available version of atmos 2 , 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 2 -rich Earth around the Sun (21% O 2 ), and another for an anoxic low methane atmosphere. The "anoxic" atmosphere is not entirely without oxygen; however, its O 2 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 Figure 10. Anoxic atmosphere comparison of full emission and deep self-reversal Lyα fill cases. Small differences in the CH4 and H2O profiles can be seen in the mesosphere, but the profiles are equal below 50 km. Similar to the oxic scenario shown in Figure 9, the concentrations of CH4 and H2O are slightly lower at the top of the atmosphere in the full emission case resulting from higher photolyzing flux from Lyαfill, while photolysis products (CO, O) are slightly greater in concentration.
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 4 and N 2 O consistent with the biological productivity of the modern Earth. The lowmethane, anoxic atmosphere assumes a 275 K surface with a 180 K stratosphere and the same surface relative humidity of of 80%. The fixed surface flux of CH 4 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 2 concentrations, 400 ppm and 1,000 ppm, respectively. In both cases the surface pressure is 1 bar and N 2 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 (>50km) 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 2 concentration for this case results in more FUV (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 4 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-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 sensitivity to the UV continuum, which is is fixed between scenarios presented here.

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 type, 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 the 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 boundary of the chromosphere and transition region yielded negligible changes to the computed UV spectrum.

Future Work
The results from this work re-emphasize 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−LT E /n LT E 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 −1 ), uniformly sampling stars with T eff from 3400-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),  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 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.

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 3 photolysis rate for the deep self-reversal case in the oxic atmosphere, though this has a minimal effect on the O 3 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.  Figure 11 but for our anoxic atmosphere scenario. These results correspond to the atmospheric mixing ratios shown in Figure 10.