# Reconstructing the Extreme Ultraviolet Emission of Cool Dwarfs Using
Differential Emission Measure Polynomials

Girish M. Duvvuri Department of Astrophysical and Planetary Sciences,
University of Colorado, Boulder, CO 80309, USA  Laboratory for Atmospheric and
Space Physics, University of Colorado, 600 UCB, Boulder, CO 80309  Center for
Astrophysics and Space Astronomy, University of Colorado, 593 UCB, Boulder, CO
80309  J. Sebastian Pineda Department of Astrophysical and Planetary Sciences,
University of Colorado, Boulder, CO 80309, USA  Laboratory for Atmospheric and
Space Physics, University of Colorado, 600 UCB, Boulder, CO 80309  Zachory K.
Berta-Thompson Department of Astrophysical and Planetary Sciences, University
of Colorado, Boulder, CO 80309, USA  Center for Astrophysics and Space
Astronomy, University of Colorado, 593 UCB, Boulder, CO 80309  Alexander Brown
Center for Astrophysics and Space Astronomy, University of Colorado, 593 UCB,
Boulder, CO 80309  Kevin France Department of Astrophysical and Planetary
Sciences, University of Colorado, Boulder, CO 80309, USA  Laboratory for
Atmospheric and Space Physics, University of Colorado, 600 UCB, Boulder, CO
80309  Center for Astrophysics and Space Astronomy, University of Colorado,
593 UCB, Boulder, CO 80309  Adam F. Kowalski Department of Astrophysical and
Planetary Sciences, University of Colorado, Boulder, CO 80309, USA  Laboratory
for Atmospheric and Space Physics, University of Colorado, 600 UCB, Boulder,
CO 80309  National Solar Observatory, University of Colorado Boulder, 3665
Discovery Drive, Boulder, CO 80303, USA.  Seth Redfield Astronomy Department
and Van Vleck Observatory, Wesleyan University, Middletown, CT 06459, USA
Dennis Tilipman Department of Astrophysical and Planetary Sciences, University
of Colorado, Boulder, CO 80309, USA  National Solar Observatory, University of
Colorado Boulder, 3665 Discovery Drive, Boulder, CO 80303, USA.  Mariela C.
Vieytes  Instituto de Astronomía y Física del Espacio, CC 67 Suc. 28, 1428
Buenos Aires, Argentina  David J. Wilson McDonald Observatory, University of
Texas at Austin, Austin, TX 78712  Allison Youngblood Laboratory for
Atmospheric and Space Physics, University of Colorado, 600 UCB, Boulder, CO
80309  Cynthia S. Froning McDonald Observatory, University of Texas at Austin,
Austin, TX 78712, USA  Jeffrey Linsky JILA, 440 UCB, University of Colorado,
Boulder, CO 80309  R. O. Parke Loyd School of Earth and Space Exploration,
Arizona State University, Tempe, AZ 85287, USA  Pablo Mauas  Instituto de
Astronomía y Física del Espacio, CC 67 Suc. 28, 1428 Buenos Aires, Argentina
Dpto. de Fisica, Facultad de Ciencias Exactas y Naturales (FCEN), Universidad
de Buenos Aires (UBA), Buenos Aires, Argentina  Yamila Miguel  Leiden
Observatory, P.O. Box 9500, 2300 RA Leiden, The Netherlands  Elisabeth R.
Newton Department of Physics and Astronomy, Dartmouth College, Hanover, NH
03755, USA  Sarah Rugheimer University Oxford, Atmospheric, Oceanic, and
Planetary Physics Department, Clarendon Laboratory, Sherrington Road, Oxford
OX1 3PU, UK  P. Christian Schneider  Hamburder Sternwarte, Gojenbergsweg 112,
D-21029 Hamburg, Germany

###### Abstract

Characterizing the atmospheres of planets orbiting M dwarfs requires
understanding the spectral energy distributions of M dwarfs over planetary
lifetimes. Surveys like MUSCLES, HAZMAT, and FUMES have collected
multiwavelength spectra across the spectral type’s range of $T_{\textrm{eff}}$
and activity, but the extreme ultraviolet flux (EUV, 100 to 912 Å) of most of
these stars remains unobserved because of obscuration by the interstellar
medium compounded with limited detector sensitivity. While targets with
observable EUV flux exist, there is no currently operational facility
observing between $150$ and $912\,\textrm{\AA}$. Inferring the spectra of
exoplanet hosts in this regime is critical to studying the evolution of
planetary atmospheres because the EUV heats the top of the thermosphere and
drives atmospheric escape. This paper presents our implementation of the
differential emission measure technique to reconstruct the EUV spectra of cool
dwarfs. We characterize our method’s accuracy and precision by applying it to
the Sun and AU Mic. We then apply it to three fainter M dwarfs: GJ 832,
Barnard’s Star, and TRAPPIST-1. We demonstrate that with the strongest far
ultraviolet (FUV, 912 to 1700 Å) emission lines, observed with _Hubble Space
Telescope_ and/or _Far Ultraviolet Spectroscopic Explorer_ , and a coarse
X-ray spectrum from either _Chandra X-ray Observatory_ or _XMM-Newton_ , we
can reconstruct the Sun’s EUV spectrum to within a factor of 1.8, with our
model’s formal uncertainties encompassing the data. We report the integrated
EUV flux of our M dwarf sample with uncertainties between a factor of 2 to 7
depending on available data quality.

††journal: ApJ††facilities: HST(STIS, COS), FUSE, XMM-Newton, Chandra,
EUVE.††software: astropy (Astropy Collaboration et al., 2013), CHIANTI &
ChiantiPy (Dere et al., 1997; Del Zanna et al., 2015), emcee (Foreman-Mackey
et al., 2013), matplotlib (Hunter, 2007), numpy (Oliphant, 2006).

##  1 Introduction

The discovery and characterization of exoplanets has been accompanied by an
increased interest in the properties of M dwarf stars as potential hosts for
habitable planets. M dwarf planetary systems are abundant, not only because M
dwarfs are $\gtrsim 70\%$ of all stars in the Milky Way (Henry et al., 2006;
Winters et al., 2015), but also because M dwarfs also have an intrinsically
high planet occurrence rate compared to their hotter and more massive siblings
(Dressing & Charbonneau, 2015). Moreover, these systems’ physical properties
benefit their detection and characterization: once for their abundance, twice
for the large transit depths of terrestrial planets projected against small
stellar radii, and thrice for the short orbital periods of planets with Earth-
comparable instellation. The _Transiting Exoplanet Survey Satellite_ (_TESS_)
is predicted to find 1300 planets orbiting M dwarfs (Ballard, 2019), roughly
10 of which will be terrestrial worlds suitable for atmospheric
characterization with the _James Webb Space Telescope_ (Barclay et al., 2018).

These M dwarf planetary systems provide a useful sample for study, but some
caution that the “habitability” potential of these systems may be poor or non-
existent (Scalo et al., 2007; Shields et al., 2016). An M dwarf is a
tempestuous host, prone to flaring (Hawley, 1993; Kowalski et al., 2009; Loyd
et al., 2018b, a), particularly when young, where a mid-to-late M dwarf’s
definition of “young” lasts for billions of years (West et al., 2008).
Compared to the Sun, M dwarfs emit a much higher fraction of their bolometric
flux in the ultraviolet regime (West et al., 2004; Jones & West, 2016). The
extreme ultraviolet region (EUV, defined here as 100 - 912 Å) is particularly
responsible for heating and ionizing the upper atmosphere of planets, dumping
energy into the system and potentially driving atmospheric escape (e.g. Sekiya
et al., 1980; Sanz-Forcada et al., 2010; Owen & Jackson, 2012; Tian & Ida,
2015; Zahnle & Catling, 2017). Any attempt to study an exoplanet atmosphere’s
evolution must be informed by the radiation field it is subject to over the
entirety of its lifetime (Penz & Micela, 2008; Claire et al., 2012; Peacock et
al., 2020). But directly measuring the EUV flux is impeded by the same
mechanism that makes it important for planet atmospheres: its interactions
with atomic hydrogen and helium mean that the interstellar medium blocks some
of the flux from this spectral region for most stars (Cox & Reynolds, 1987;
France et al., 2019). This problem is exacerbated for M dwarfs since the
closest M dwarfs with observable EUV flux either have noisy data from the
_Extreme Ultraviolet Explorer_ (_EUVE_) or no data at all (Craig et al., 1997;
Linsky et al., 2014; France et al., 2016), and there is no presently available
dedicated EUV observatory to remedy the situation (France et al., 2019).

In the absence of direct observation, we must turn to theoretical inference.
Peacock et al. (2019a) and Peacock et al. (2019b) use the PHOENIX 1D stellar
atmosphere code (Hauschildt, 1993; Hauschildt & Baron, 2006; Baron &
Hauschildt, 2007) to model the non-LTE radiative transfer through the
chromospheres and transition regions of M dwarfs but do not include a corona.
Fontenla et al. (2016) adjusts the temperature and pressure profiles of a 1D
stellar atmosphere until the model agrees with the available spectral data,
but this takes time to do well and has to be specific to each star. These
semi-empirical methods require quasi-simultaneous observations from optical to
X-ray wavelengths. All known atomic and molecular processes and species have
to be taken into account in each layer, solving the NLTE coupled system to
match observations of many lines and continua across the spectral range. This
requires a reliable atomic database and laborious fine-tuning to be
successful.

Taking a more empirical approach, Linsky et al. (2014) identifies correlations
between Lyman-$\alpha$ and EUV flux, while Youngblood et al. (2017) identifies
correlations between far-ultraviolet (FUV, 912 to 1700 Å) lines and the
Lyman-$\alpha$ flux, chaining these correlations to the Linsky et al. (2014)
relations to predict the EUV flux in turn. A drawback of this method is that
the uncertainty on each correlation introduces scatter into the predicted EUV
flux while the sample is insufficiently large to investigate the effects of
both effective temperature $T_{\textrm{eff}}$ and stellar activity. France et
al. (2018) correlates certain FUV lines with the EUV flux between 90 and 360 Å
directly, leading to much less scatter in the predicted flux and accounting
for both $T_{\textrm{eff}}$ and stellar activity in their sample, but this
still leaves us with $\sim 600$ Å of EUV flux to estimate.

These limitations of existing methods lead us to use the differential emission
measure (DEM), a technique for EUV spectral synthesis adapted from an earlier
technique called the emission measure distribution. Pottasch (1963) defined
the emission measure distribution
$\equiv\frac{n_{\textrm{O}}}{n_{\textrm{H}}}\int n_{e}^{2}\,ds$ as the
integral of the electron number density squared ($n_{e}^{2}$) along the line
of sight $s$ weighted by the relative abundance of oxygen to hydrogen
($\frac{n_{\textrm{O}}}{n_{\textrm{H}}}$), to describe the plasma environment
of the upper layers of the Sun’s atmosphere. This assumed that the Sun’s upper
atmosphere could be approximated as a series of spherical shells of increasing
temperature, and all emission lines were produced by collisional excitation
and spontaneous radiative decay within restricted spatial regions.

As this picture of spherical symmetry broke down, the differential emission
measure was developed to keep the same 1-dimensional simplification to
temperature but account for the spatial ambiguity of a photon’s origin (see
Mariska 1992 for a detailed overview of the method’s history). The
differential emission measure uses a similar integral expression over a
limited temperature range to estimate the density and temperature environment
of ions emitting an observed line, allowing one to then use those
environmental conditions to estimate the flux from emission lines that cannot
be observed but should be emitted by the same parcel of plasma. When UV
detector technology was in its infancy and instruments had poor flux
calibration, differing by factors of $\gtrsim 2$ in different wavelength
regimes, the DEM could be used to estimate the subset of solar emission line
fluxes with poor data from other lines that were thought to have more accurate
and precise data (Warren et al., 1998). While the state of solar EUV data has
improved, the opacity of the interstellar medium and low sensitivities of
previous and current EUV-capable observatories present a similar spectral
synthesis problem for distant stars. Variations of the DEM have been applied
to other stars like AU Mic by Pagano et al. (2000), $\alpha$ Centauri A and B
by Ayres (2014), and HD 209458 by Louden et al. (2017) to infer the EUV flux
from these stars. Sanz-Forcada et al. (2011) developed scaling relations
between X-ray and EUV fluxes by applying the DEM method to a large sample of
stars, but the paper’s sample had few M dwarfs and lacked enough UV data to
constrain the lower temperature end of the DEM for most of their stars.

In this paper we characterize our uncertainties in fitting the DEM and
propagate them to our predictions of the EUV flux from M dwarfs. Our physical
assumptions and setup are similar to the method described and used by Warren
et al. (1998) to model the EUV irradiance of the Sun, described in Section §2.
The specifics of our implementation are described in Section §3 and we test
our method against data from the Sun in Section §4. In Section §5 we apply our
method to AU Mic, a $\sim 10-20$ Myr old M1 star at a distance of 9.979 pc
(MacGregor et al., 2013; Plavchan et al., 2020). We compare our DEMs of the
Sun and AU Mic to previous DEMs published in the literature and available in
the CHIANTI atomic database (Dere et al., 1997; Del Zanna et al., 2015) in
Section §6.

We compare our predicted spectra for the Sun and AU Mic to data in detail in
Section §7 and in Section §8 we apply our method to different case studies: GJ
832, a planet-hosting M2 V that has predicted EUV fluxes from Linsky et al.
(2014) and semi-empirical models from both Fontenla et al. (2016) and Peacock
et al. (2019b); Barnard’s Star, a $\sim 10$ Gyr old M4 with a candidate planet
(Ribas et al., 2018), with contemporaneous X-ray and FUV data during
quiescence and a flare (France et al., 2020); and TRAPPIST-1, an ultracool
dwarf which hosts at least seven planets (Gillon et al., 2017) and tests our
ability to fit the DEM in an extremely low S/N regime (Wilson et al. submitted
). Our work shows that with _Hubble Space Telescope_ (_Hubble_ or _HST_)
measurements of a few FUV emission line fluxes and a coarse X-ray spectrum
from _Chandra_ or _XMM-Newton_ , we can estimate the EUV spectrum with
meaningful uncertainties for any star whose EUV flux is dominated by emission
lines from the optically thin regions of the star’s upper atmosphere.

##  2 Differential Emission Measure

The following description of the DEM is adapted from Warren et al. (1998).
Many other formulations of the DEM and similar techniques exist, and Mariska
(1992) explains them in more detail. Given an optically thin plasma in a
collisionally dominated time-independent equilibrium with negligible
collisional de-excitation, the radiance of a wavelength transition is given by

$\displaystyle I_{ul}$ $\displaystyle=\frac{1}{4\pi}\int_{\mathrm{line-of-
sight}}n_{u}A_{ul}\frac{hc}{\lambda_{ul}}\,ds\;[\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{sr}^{-1}],$
(1)

where $ul$ signifies a transition from an upper state $u$ to a lower state
$l$, $A_{ul}$ is the Einstein rate coefficient of the transition,
$\lambda_{ul}$ is the wavelength of the transition, $h$ is Planck’s constant,
and $c$ is the speed of light in a vacuum. This quantity is not a spectral
density because it captures all of the emission from the spontaneous radiative
decay without describing a line profile. We can rewrite this integral as

$\displaystyle I_{ul}$
$\displaystyle=\int_{T}G_{ul}(T)\cdot\Psi(T)\,dT\;[\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{sr}^{-1}],$
(2)

where

$\displaystyle G_{ul}(T)$
$\displaystyle=\frac{n_{u}}{n_{\mathrm{ion}}}\frac{n_{\mathrm{ion}}}{n_{\mathrm{element}}}\frac{n_{\mathrm{element}}}{n_{\mathrm{H}}}\frac{1}{n_{e}}\frac{A_{ul}hc}{4\pi\lambda_{ul}}\;[\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{3}\,\mathrm{sr}^{-1}]$
(3)

is the transition’s emissivity contribution function and the differential
emission measure is

$\Psi(T)=n_{e}n_{\mathrm{H}}\frac{ds}{dT}\;[\mathrm{cm}^{-5}\,\mathrm{K}^{-1}].$
(4)

$G_{ul}(T)$, the emissivity contribution function, describes the volume
integrated power of a parcel of gas as a function of temperature. The function
can be computed with a few ingredients: a stellar abundance to give us the
ratio of the number density of any particular element’s atoms to the number
density of hydrogen atoms $\frac{n_{\mathrm{element}}}{n_{\mathrm{H}}}$, the
assumption of collisionally dominated equilibrium (i.e. coronal equilibrium)
and CHIANTI to give us the population fraction of any particular upper state
of an ion $\frac{n_{u}}{n_{\mathrm{ion}}}$ and the population fraction of each
ion per element $\frac{n_{\mathrm{ion}}}{n_{\mathrm{element}}}$, an assumed
local density $n_{e}$, and laboratory measurements or theoretical calculations
of the atomic data $A_{ul}$ and $\lambda_{ul}$. We follow Del Zanna et al.
(2002) in using a constant electron pressure $P_{e}$ to define
$n_{e}(T)=\frac{P_{e}}{k_{B}T}$, where $k_{B}$ is the Boltzmann constant. This
single pressure will not be applicable to the entire temperature domain, but
errors in the $G(T)$ function can be partially compensated for by the
$\Psi(T)$ function as long as the errors are largely a function of temperature
and do not vary significantly across lines formed at the same temperature.

The differential emission measure, $\Psi(T)$, describes the density and
temperature structure along the line of sight, common to all transitions we
observe from the chromosphere, transition region, and corona. Under our
assumptions that the ions are predominantly populated by collisions and
depopulated by spontaneous emission, the flux observed is proportional to the
collision rate. The differential emission measure resembles a reaction rate,
$n_{e}\cdot n_{\textrm{H}}$, weighted by $\frac{ds}{dT}$ which measures how
much of the path length $s$ is at a temperature $T$. In emission measure
studies of other stars, a volume emission measure is commonly employed that
predicts a flux and includes factors of the stellar radius and solid angle
filling factor of the emitting plasma. We adopt the line-of-sight approach to
be able to compare the DEMs of very different stars to each other and to solar
surface features.

For each emission line there is a formation temperature $T_{f}$ that maximizes
the product $G_{{ul}}(T)\cdot\Psi(T)$, and since the emissivity function
$G_{ul}(T)$ tends to be very narrowly peaked, the bulk of of the observed line
flux is emitted by plasma at $\approx T_{f}$. By measuring the observed line
intensities of transitions with a known $G_{ul}(T)$, we can constrain the
value of $\Psi(T)$ within the vicinity of the lines’ formation temperatures
$T_{f}$. Amassing a list of observable transitions over a sufficiently wide
range of $T_{f}$ allows us to fit for the parameters of an assumed functional
form describing $\Psi(T)$ across the temperature domain of the upper stellar
atmosphere. With $\Psi(T)$ in hand and atomic data to construct $G_{ul}(T)$
for the transitions we have not observed but seek to estimate, we can
reconstruct the optically thin emission of the chromosphere, transition
region, and corona. With the exception of the recombination continua addressed
in Section §4.3, optically thin emission lines contribute the majority of the
EUV flux from an M dwarf.

##  3 Implementation

Figure 1: The left panel is a heatmap where the shading represents the value
of the emissivity contribution function $G_{\lambda}(T)$ as a function of
temperature along the horizontal axis and as a function of wavelength along
the vertical axis. This matrix is generated from the CHIANTI atomic database
by specifying the atomic abundances and the electron density as a function of
temperature $n_{e}(T)$. The right panel multiplies the $G_{\lambda}(T)$ matrix
by the DEM, $\Psi(T)$, of the Sun shown as a blue line in Figure 6. The
colorbars of both panels are cut off at the 60${}^{\textrm{th}}$ and
95${}^{\textrm{th}}$ percentiles to highlight the strongest lines. While the
$G_{\lambda}(T)$ matrix shows the atomic data, the right panel represents the
temperature integrand for which lines are actually emitted by stars. EUV lines
are largely formed at temperatures between $10^{5}$ and $10^{6.5}$ K,
requiring FUV measurements to constrain the low temperature end and X-ray
measurements and/or coronal FUV semi-forbidden transitions of highly ionized
iron to constrain the high temperature end.

Figure 2: A sample of the emissivity functions $G_{ul}(T)$ and
$G_{\lambda}(T)$ described in Sections §2 and §3. The solid lines show the
summed emissivity functions for some of the strongest transitions of the
labeled ion observable by _Hubble_ with STIS and COS. The dotted lines show
the summed emissivity functions $G_{\lambda}(T)$ for typical wavelength bins
in a _Chandra_ spectrum of an M dwarf, adding up the $G_{ul}(T)$ functions of
all emission lines that fall within the wavelength bin. The ion emissivities
are constructed by adding together the emissivities of all lines emitted by
that ion, for example both lines of the N V 1239/1243 doublet. The X-ray
wavelength bin emissivities are constructed similarly by adding together the
emissivities of lines from multiple ions that are emitted at similar
wavelengths. This lets us use coarser resolution X-ray spectra for fainter
stars but weakens the temperature specificity of individual flux measurements.
Note the contribution function for a semi-forbidden transition of Fe XXI at
1354 Å with a formation temperature at $\sim 10^{7}$ K. Its emissivity peaks
at a factor of $10^{3}$ times less than other typically observable FUV lines,
making it unobservable for most quiescent M dwarfs. Observing or obtaining an
upper limit for the flux of this line provides an additional constraint on
$\Psi(T)$ at high temperatures.

We use CHIANTI 8.0.7 (Dere et al., 1997; Del Zanna et al., 2015) to calculate
the $G_{ul}(T)$ functions for all the transitions in the database assuming the
elements have a solar coronal abundance (Schmelz et al., 2012). We calculate
these emissivity contribution functions across a temperature range from
$10^{4}$ to $10^{8}$ K for multiple assumed electron pressures ranging from
$P_{e}=10^{12}$ to $10^{25}$ $k_{B}$ cm-3 K. The majority of coronal emission
lines are not strongly sensitive to density, but we test the variation in the
predicted flux as a function of the $P_{e}$ used to calculate $G_{ul}(T)$ in
Section §4.1. We assume that $\Psi(T)$ is well-described by

$\log_{10}\Psi(T)=\sum_{n=0}^{5}c_{n}\mathbf{T_{n}}\left(\frac{\log_{10}T-6}{2}\right)$
(5)

where $\mathbf{T_{n}}(x)$ are the Chebyshev polynomials of the first kind, and
their argument is shifted and scaled to transform the domain to the interval
$[-1,1]$. We use the Chebyshev polynomials following the previous work of
Louden et al. (2017) and because they form an orthonormal basis. Given the
coefficients $c_{n}$, and the list of emissivity contribution functions
$G_{ul}(T)$, we can generate a full spectrum by summing the contribution
functions of all emission lines within a wavelength bin centered on a
wavelength $\lambda$ with a width $\Delta\lambda$ to get a wavelength-specific
contribution function

$G_{\lambda}(T)=\sum_{\mathrm{transitions}}G_{ul}(T)\
\mathbf{if}\left[\lvert\lambda_{ul}-\lambda\rvert\leq\Delta\lambda\right].$
(6)

We then scale the temperature integral of Equation 2 to predict the observed
flux density in each wavelength bin, assuming the bin is wide enough to
contain the entire line profile

$F_{\lambda}=\frac{\pi}{\Delta\lambda}\left(\frac{R^{2}_{\star}}{d^{2}}\right)\int_{T}G_{\lambda}(T)\cdot\Psi(T)\,dT\;[\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\text{\AA}].$
(7)

The scaling factor assumes that the solid angle emitting the flux is
$\approx\pi\left(\frac{R_{\star}}{d}\right)^{2}$ steradians, which is
approximate because the corona extends beyond the stellar radius. We create a
matrix of $G_{\lambda}(T)$ with the wavelength axis at a constant resolving
power $R=\frac{\lambda}{\Delta\lambda}=500$ between 1 and 2000 Å, and the
temperature axis as 2000 logarithmically spaced points between $10^{4}$ and
$10^{8}$ K (see Figure 1). Fitting for the coefficients $c_{n}$ and combining
the $\Psi(T)$ model with this matrix allows us to generate a high-resolution
spectrum, but since the DEM makes no prescription for line shape, the line
profiles are all Dirac-$\delta$ functions, which is why we then divide by the
wavelength bin width $\Delta\lambda$ to get the observed flux density.
Comparing this model to a real stellar spectrum is only reasonable at a low
enough resolution such that the entirety of the line profile is contained
within each resolution element. The $R=\frac{\Delta\lambda}{\lambda}=500$
$G_{\lambda}(T)$ matrix can be downsampled to whatever resolution is required
to contain the line widths of any spectral data used for comparison.

By combining either the $G_{ul}(T)$ functions or the $G_{\lambda}(T)$ matrix
with the polynomial coefficients $c_{n}$, we have a generative model for a
list of integrated line fluxes or a low-resolution spectrum respectively. To
get a usefully constrained model, we need data that covers the full
temperature domain. Figure 2 shows the $G_{\lambda}(T)$ functions for the
wavelength bins of a typical _Chandra_ spectrum, where each bin peaks at a
slightly different temperature but spans $10^{6}$ and $10^{7.5}$ K, and the
$G_{ul}(T)$ functions for the strongest optically thin FUV lines accessible in
a _Space Telescope Imaging Spectrograph_ (STIS) and _Cosmic Origins
Spectrograph_ (COS) spectrum from _Hubble_. There is significant overlap in
the contribution functions of these transitions near $10^{5}$ K, but they
spread out far enough to constrain $\Psi(T)$ between $10^{4}$ and $10^{5.5}$
K. Each line flux measurement can be used to derive an average value of
$\Psi(T)$ near the formation temperature $T_{f}$ of the transition

$\overline{\Psi}(T_{f})=\frac{F_{\textrm{line}}}{\pi\left(\frac{R^{2}_{\star}}{d^{2}}\right)\int_{T}G_{ul}(T)dT}$
(8)

and any individual wavelength bin’s flux density can provide a similar
constraint by substituting $G_{\lambda}(T)$ in for $G_{ul}(T)$ and dividing by
the wavelength bin width $\Delta\lambda$. We do not fit to these averages
because we can directly compare our predicted fluxes to the data, but the
averages are useful for visualizing how an individual flux measurement
constrains the DEM.

Using the affine-invariant Markov Chain Monte Carlo sampler implemented in the
Python package emcee (Foreman-Mackey et al., 2013), we fit the coefficients
$c_{n}$ (see Equation 5) using a combination of the available X-ray data and
integrated FUV line fluxes. Since the uncertainties on the emissivities are
unknown and we have little a priori information on how to characterize the
systematic uncertainties associated with this method, we assume that the
variance is boosted by a scaled multiple of the predicted flux, making our
log-likelihood

$\ln{\mathcal{L}}=\sum_{i}\ln\left(\frac{1}{\sqrt{2\pi\left(\sigma_{y_{i}}^{2}+\left(s\cdot
f(x_{i})\right)^{2}\right)}}\right)-\left(\frac{y_{i}-f(x_{i})}{2\sqrt{\sigma_{y_{i}}^{2}+\left(s\cdot
f(x_{i})\right)^{2}}}\right)^{2}$ (9)

where $f(x_{i})$ is the model prediction, $y_{i}$ is the data,
$\sigma_{y_{i}}$ is the Gaussian uncertainty of the data, and $s$ is the free
parameter that characterizes these unknown systematic uncertainties (which are
assumed to be independent of the data and temperature). Some contributions to
$s$ are likely to be errors in stellar parameters like the stellar abundance,
deviations of level populations from true collisional equilibrium or
variations in the relative abundances along the line of sight, the departure
from being perfectly optically thin $\tau=0$, and the spatial inhomogeneity of
the emitting plasma. This form of the likelihood is independent of the two
types of data described above, allowing us to mix together combinations of
line fluxes and spectra in different wavelength regimes, so long as we ensure
that these do not overlap to count the same data twice.

We incorporate Bayesian priors on individual parameters to modify the
likelihood evaluated in Equation 9. We sample $\log_{10}(s)$ uniformly between
$-2$ and $2$. The mean value of the DEM is set by $c_{0}$, which is sampled
uniformly between 20 to 26 and then exponentially cut off beyond those bounds.
These boundary values were chosen to limit the DEM to physical expectations
for $10^{5}\lesssim n_{e}\lesssim 10^{17}$ cm-3, and path-length
$10^{8}\lesssim ds\lesssim 10^{11}$ cm. The remaining coefficients $c_{n}$ are
sampled uniformly within the bounds $\pm 100$, and then we also require that
the base-10 logarithm of the final polynomial be positive at $\log_{10}T=6$ to
prevent unphysically small DEMs and that the derivative be negative at the
lower bound $T=10^{4}$ K to reflect the higher amount of material in the
photosphere compared to the chromosphere. These priors extend generously
beyond physically realistic DEM shapes, for example they do not require the
DEM to go to 0 at high temperatures, allowing for an infinitely extended
corona. Data constrain the parameter distributions to factors of a few at
most, with the $s-$factor systematic uncertainty typically restricted to the
interval $0.1<s<1$.

##  4 Testing the DEM Method Against the Sun

To test our implementation of the DEM on solar data, we use the Solar
Irradiance Reference Spectra (SIRS) published by Woods et al. (2009). This is
a disk-integrated spectrum of the quiescent Sun assembled from measurements
collected during the 2008 minimum of the solar activity cycle at 1 Å
resolution. Referring to a list of the lines used for the DEM fitting in
Warren et al. (1998), and making a point to select the FUV lines most likely
to be detected in _Hubble_ observations of M dwarfs, listed in Table 1, we
measure their fluxes in this spectrum by subtracting the continuum and
integrating line profiles. Then we selected the X-ray data between $5$ and
$50$ Å, comparable to the regions observed by the _Chandra_ -_X-ray_
_Observatory_ and _XMM_ -_Newton_ , and left the spectrum at its original
resolution of 1 Å wavelength bins, $R=\frac{\lambda}{\Delta\lambda}\leq 50$.
This combination of line fluxes and an X-ray spectrum is the same type of data
we use for M dwarfs discussed later in this work. Table 1 also lists the
integrated fluxes of EUV lines measured from the Woods et al. (2009) spectrum,
used in the test described in Section §4.3. The SIRS did not provide error
bars, but we assigned errors such that we had three versions of the data with
S/N $=$ 1, 10, and 100 to test the sensitivity of the fitting to S/N. The true
errors vary across the observations from different instruments and wavelength
ranges assembled by Woods et al. (2009), but never exceed 10% at instrument
native resolutions which are much finer than the 1 Å bins used here.

Table 1: Integrated fluxes of optically thin lines measured in the Solar
Irradiance Reference Spectrum (Woods et al., 2009) compared to the DEM
predictions.

Ion | Wavelengths | $\log_{10}T_{f}$ | Observed Flux | FUV/X-ray DEM | W/o Anomalous Ions DEM
---|---|---|---|---|---
| [Å] | $\log_{10}(\textrm{[K]})$ | [$10^{-2}$ erg s-1 cm-2] | [$10^{-2}$ erg s-1 cm-2] | [$10^{-2}$ erg s-1 cm-2]
C II | 1335.7 | 4.4 | 14.9 | 11.0 | 11.9
C III | 1175.7 | 4.8 | 4.73 | 7.78 | 3.36
C IV aaThese FUV transitions were not used to fit the “Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 1548.2, 1550.7 | 5.0 | 12.0 | 7.95 | 2.17
N V aaThese FUV transitions were not used to fit the “Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 1238.8, 1242.8 | 5.3 | 1.56 | 1.56 | 0.419
Ne VIIbbThese EUV transitions were used to fit the “‘Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 465.2 | 5.7 | 1.47 | 3.07 | 1.60
Ne VIIIbbThese EUV transitions were used to fit the “‘Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 770.4, 780.3 | 5.8 | 1.80 | 2.83 | 1.92
O IIIbbThese EUV transitions were used to fit the “‘Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 508.2, 525.8, 599.6, 703.9 | 4.9 | 2.71 | 6.63 | 2.09
O IVbbThese EUV transitions were used to fit the “‘Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 554.5, 787.7, 790.2 | 5.2 | 5.85 | 19.9 | 4.95
O VbbThese EUV transitions were used to fit the “‘Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 629.7, 760.4 | 5.4 | 6.47 | 21.4 | 5.77
Si III | 1206.5 | 4.5 | 6.83 | 22.6 | 18.4
Si IVaaThese FUV transitions were not used to fit the “Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 1393.8 | 4.9 | 3.72 | 2.24 | 0.827
Si XIIbbThese EUV transitions were used to fit the “‘Fit with EUV Lines and without Anomalous Ions” model shown in Figure 6 as a red solid line. | 499.4 | 6.3 | 0.699 | 0.921 | 1.03

Note. — In cases where multiple transitions are listed for the same ion, the
reported flux is the summed flux across all listed transitions.

###  4.1 Pressure Sensitivity

Across this broad range of temperatures, no single electron density or
pressure will accurately describe the environmental conditions of the plasma
emitting the observed flux we are using to fit the DEM or the unobserved EUV
flux we are trying to predict. However we must assume some function for the
electron density $n_{e}(T)$ to calculate emissivities if we want to fit the
DEM at all. Updating the emissivity calculation iteratively would be
computationally prohibitive and still fail to accurately describe detailed
non-equilibrium physics. By generating multiple emissivity matrices across a
broad range of electron pressures, $P_{e}=10^{12}$ to $10^{25}$ $k_{B}$ cm-3
K, and fitting a DEM to the solar data with each matrix, we test the
sensitivity of the DEM shape and calculate the variation in the predicted EUV
flux as a function of assumed pressure. Figure 3 shows a representative sample
of these DEMs, which vary only slightly for pressures lower than $10^{20}$
$k_{B}$ cm-3 K and are consistent with each other to within 1$\sigma$
variations of the DEM shape. The horizontal lines are the average
$\overline{\Psi}(T_{f})$ values. To test if any particular model is a
statistically significant improvement over the others, we compare the models’
values of the Bayesian Information Criterion (BIC, Schwarz 1978). Kass &
Raftery (1995) demonstrates that the BIC is related to the natural logarithm
of the Bayes factor, such that a $\Delta BIC=1$ implies the more negative
model is $e$ times more likely than the higher one. The BIC is evaluated with
the equation

$\textrm{BIC}=k\ln(n)-2\ln(\hat{\mathcal{L}})$ (10)

where $k$ is the number of model parameters, $n$ is the number of datapoints,
and $\hat{\mathcal{L}}$ is the maximum-likelihood of the model. This criterion
penalizes a higher number of parameters, and the model significance increases
as the BIC decreases. All models in this comparison have the same number of
parameters, but we also use the BIC later in Section §4.5 to test our method’s
sensitivity to polynomial degree.

Table 2 compares each pressure model’s BIC, estimated systematic uncertainty
characterized by the $s$-factor, and EUV flux integrated from 100 to 912 Å. In
the middle of our pressure range, from $10^{17}$ to $10^{20}$ $k_{B}$ cm-3 K,
the predicted integrated fluxes are consistent with each other to within
1$\sigma$, but the BIC clearly favors the $10^{19}$ model. We adopt the
$P_{e}=10^{19}$ $k_{B}$ cm-3 K emissivity matrix for other tests of the Sun
DEM model moving forward. At pressures higher than $10^{21}\ k_{B}$ cm-3 K,
the DEM shape and predicted fluxes change drastically, likely because the
plasma is optically thick and collisional de-excitation can no longer be
ignored. The base of the solar chromosphere is at a pressure of $\sim 10^{20}$
$k_{B}$ cm-3 K (Mariska, 1992), so a model DEM that assumes the entire upper
atmosphere is at photospheric pressure is bound to be unphysical.

For all other stars, we adopt the same approach of fitting the star’s DEM with
each pressure separately and choosing the model with the best likelihood. We
caution that it is unphysical to interpret these “best” pressures as
representative of a specific region in the stellar atmosphere, and that they
should be seen as the most useful average for implementing the DEM and nothing
more. Future work could involve testing the DEM with temperature-pressure
profiles from stellar atmosphere models to see if this improves the accuracy
and precision of the estimated spectrum.

Figure 3: The DEMs of the Sun fit to the same data, a combination of FUV
emission lines and X-ray spectra, but with emissivity functions calculated at
different electron pressures $P_{e}$. The horizontal bars represent
constraints on the average DEM value $\overline{\Psi}(T_{f})$ in the vicinity
imposed by the emissivity functions $G(T)$ for ion line fluxes or wavelength
bins as described in Section §3, with the horizontal extent of the lines
representing the full width at half maximum for the emissivity function. These
constraints transformed from a flux to an average DEM value are only
approximately illustrative of the fit quality. For a true comparison of the
DEM prediction to the data used to fit the model, see plots of the predicted
line fluxes and X-ray spectra in Figures 7 and 8 respectively. The solid lines
represent the median DEM value as a function of temperature from the posterior
samples for $\Psi(T)$ while the shaded region encloses the
16${}^{\textrm{th}}$ to 84${}^{\textrm{th}}$ percentile values of the DEM. The
red, blue, and brown models correspond to DEM models with emissivities
evaluted at electron pressures $P_{e}=$ $10^{13}$, $10^{19}$, and $10^{25}$
$k_{B}$ cm-3 K respectively. The green dashed bars and black solid bars
represent the flux constraints from the FUV lines and X-ray spectra used to
fit these models, while the dotted purple bars represent the flux constraints
from EUV lines listed in Table 1 but were not used to fit these models. For
visualizing these constraints we use the emissivities calculated by assuming
the electron pressure $P_{e}=10^{19}\,k_{B}$ cm-3. Table 2: The Bayesian
Information Criterion (BIC) for each pressure evaluated against the FUV line
fluxes and the X-ray spectrum used to fit the model.

Log Electron Pressure $\log_{10}P_{e}$ | BIC | Integrated EUV Flux | $\log_{10}s_{\mathcal{L}_{\mathrm{max}}}$ | $\log_{10}s$
---|---|---|---|---
$\log_{10}\left([k_{B}\ \mathrm{cm}^{-3}\ \mathrm{K}]\right)$ | — | [ergs s-1 cm-2] | — | —
$12$ | $-447.8$ | $4.9^{+4.2}_{3.4}$ | $-0.2090$ | $-0.1671^{+0.0739}_{-0.0632}$
$13$ | $-448.6$ | $4.6^{+4.3}_{-3.1}$ | $-0.1852$ | $-0.1699^{+0.0731}_{-0.0626}$
$14$ | $-450.1$ | $4.3^{+4.5}_{-2.8}$ | $-0.1929$ | $-0.1757^{+0.0735}_{-0.0619}$
$15$ | $-449.8$ | $4.5^{+4.2}_{-3.1}$ | $-0.2166$ | $-0.1741^{+0.0726}_{-0.0624}$
$16$ | $-452.3$ | $4.1^{+3.6}_{-2.8}$ | $-0.2249$ | $-0.1783^{+0.0733}_{-0.0631}$
$17$ | $-456.2$ | $3.6^{+3.5}_{-2.3}$ | $-0.2137$ | $-0.1909^{+0.0720}_{-0.0619}$
$18$ | $-460.7$ | $3.5^{+3.3}_{-2.1}$ | $-0.2311$ | $-0.2019^{+0.0697}_{-0.0615}$
$19$ | $-464.7$ | $3.5^{+3.1}_{-2.1}$ | $-0.2424$ | $-0.2178^{+0.0692}_{-0.0596}$
$20$ | $-463.3$ | $3.8^{+3.3}_{-2.4}$ | $-0.2515$ | $-0.2125^{+0.0709}_{-0.0603}$
$21$ | $-455.1$ | $6.2^{+6.1}_{-4.1}$ | $-0.2108$ | $-0.1804^{+0.0729}_{-0.0623}$
$22$ | $-426.2$ | $20^{+20}_{-17}$ | $-0.0881$ | $-0.0591^{+0.0873}_{-0.0725}$
$23$ | $-376.7$ | $59^{+94}_{-59}$ | $0.0828$ | $0.1254^{+0.1199}_{-0.0962}$
$24$ | $-367.6$ | $140^{+220}_{-140}$ | $-0.0113$ | $0.0413^{+0.1069}_{-0.0826}$
$25$ | $-364.1$ | $200^{+370}_{-200}$ | $0.0244$ | $0.0473^{+0.1067}_{-0.0809}$

Note. — The BIC penalizes model parameters by $k\ln(n)$ where $k$ is the
number of parameters being fit and $n$ is the number of data points being fit
to. An increasingly negative BIC indicates a better fit. In this case, the
most preferred models are the $P_{e}=10^{19}$ and $10^{20}$ $k_{B}$ cm-3 K
models respectively. We also show the value of $\log_{10}s$ for the maximum
likelihood sample from the posterior and the median $\pm 1\sigma$ confidence
interval for $\log_{10}s$.

###  4.2 Sensitivity to S/N

Figure 4: A comparison of fitting the DEM to the same data, a combination of
FUV line fluxes and X-ray spectra from the Woods et al. (2009) quiescent Sun
spectrum, but with different errors assigned to vary the signal-to-noise ratio
(S/N). As shown before in Figure 3, the solid lines and shaded regions
represent the different DEM models while the horizontal bars represent
constraints imposed by the flux measurements. The red, green, and blue models
represent fitting the DEM at S/N values of 1, 10, and 100 respectively. As
shown previously in Figure 3, the green dashed bars and black solid bars
represent the flux constraints from the FUV lines and X-ray spectra used to
fit these models. The dotted purple bars represent the flux constraints from
EUV lines listed in Table 1 and which were not used to fit these models.
Figure 5: The corner plot for the parameter distributions when fitting the Sun
at S/N = 100 shows the median value of $\log_{10}s$ to be
$-0.22^{+0.07}_{-0.06}$ which translates to a model uncertainty of $60\%$ for
the predicted flux. This uncertainty dominates the fitting over the
uncertainty associated with the data itself.

Figure 4 compares our DEMs fit to the Sun’s data with S/N = 1, 10, and 100.
For the higher S/N models, the variance is dominated by the uncertainty on the
predicted flux, parameterized by the $s-$factor and independent of the shape
of the DEM, so changing the S/N of the data used to fit the DEM has little
consequence. At S/N = 1, the data uncertainty dominates and the model
percentile ranges shift dramatically, with no overlap with the higher S/N
models until the higher temperature regimes constrained by many X-ray spectrum
fluxes. The low DEM values at FUV temperatures are likely a consequence of the
prior requiring a negative derivative at the start of the temperature domain.
As mentioned earlier, $s$ is a measure of our average temperature-independent
systematic uncertainty that combines the uncertainties on the emissivities
with anything else intrinsic to our method’s assumptions and approximations.
Figure 5 shows us that $s=0.63$ in the solar case with S/N$=100$. This
indicates that we should assume the systematic uncertainty on any predicted
line flux is roughly 60% of the predicted value. Some of the systematic
uncertainty may be attributed to fact that the Woods et al. (2009) spectrum
combines observations of quiescence from different instruments taken at
different times, which is a problem that will affect most stellar observations
and needs to be accounted for in modelling uncertainties.

###  4.3 Including EUV Data in the Fitting and Excluding Anomalous Ions

With the Sun, we can refer to the EUV lines observed by Woods et al. (2009) in
the SIRS, allowing us to see how much information we are losing about the DEM
in the stellar case where EUV data is not available. Including the EUV lines
gives us more temperature coverage and allows us to exclude three ions from
the Na-like and Li-like isolectronic sequences: N V, C IV, and Si IV. Del
Zanna et al. (2002) showed that a DEM informed by these ions significantly
overpredicts the flux of other transitions because of an anomaly in the
CHIANTI ionization equilibrium calculations for these isoelectronic sequences
compared to other ions for the same plasma environmental conditions. The
factor of discrepancy is not constant across all transitions and densities, so
it cannot be corrected for by a consistent known number. This discrepancy
constitutes a significant systematic uncertainty that cannot be avoided when
fitting the DEM to faint stars with few strong measurable lines that are not
from these anomalous ions. Including line fluxes from multiple transitions of
other ions can help mitigate the influence of the anomalous ions, but upper
limits can still help if the star is too faint to measure these lines.

To help characterize the magnitude of these discrepancies, we include EUV
lines from 7 ions that are currently unobservable for our M dwarf sample of
interest, listed in Table 1, and drop the anomalous ions (retaining some of
the FUV lines and the X-ray spectrum) when fitting the model labeled “Fit with
EUV Lines and without Anomalous Ions” in Figures 6, 7, 8, and 9. This new DEM
model shifts down by a factor of $\sim 5$ in between $10^{5}$ and $10^{6}$ K
(see Figure 6) to match the EUV line fluxes (see Figure 7) that are not
informing our stellar-analogous DEM. The DEMs agree with each other at the
higher temperatures constrained by the X-ray spectra resulting in nearly
identical predictions in that spectral regime (see Figures 6 and 8). The
predicted line fluxes from both models are compared to the data in Table 1 and
Figure 7, and highlight the problem of the anomalous ions. The model without
EUV lines and including the anomalous ions predicts the FUV fluxes reasonably
well, adopting a compromise position between FUV lines formed at similar
temperatures that have discrepant DEM constraints (see Figure 6), but this
compromise still overestimates the flux of the EUV lines by up to a factor of
5. When applying the DEM to M dwarfs without EUV data, we include the ions
with anomalous CHIANTI emissivities because these are the strongest lines
available and we cannot afford to simply exclude them. Measuring upper limits
for the fluxes of other transitions formed at similar temperatures can
mitigate the influence of the anomalous ions when combined with the $s-$factor
uncertainty, as demonstrated in our modeling of AU Mic in section 5.

In Figure 7, the error bars associated with the plotted data point incorporate
both the posterior distributions of the DEM shape and the $s-$factor
uncertainty of the “Fit with EUV Lines and without Anomalous Ions” model. We
do this by drawing randomly from the posterior sample of the MCMC fitting to
generate a sample $\Psi(T)$ using the Chebyshev coefficients $c_{n}$ (see
Equation 5) with an associated $s-$factor. The $\Psi(T)$ is combined with
$G(T)$ to predict the flux of an observed data point, $y_{i}$, giving a model
flux $f(x_{i})$. Multiple sample draws in this fashion would only represent
the uncertainty associated with the DEM shape and exclude the $s-$factor. To
include the model-intrinsic uncertainty, we draw randomly from the Gaussian
distribution $\mathcal{N}(\mu=f(x_{i}),\sigma={s\cdot f(x_{i})})$ and record
the prediction from the flux distribution parameterized by a single MCMC
posterior sample. This process is executed with $N=5\times 10^{4}$ draws from
the model parameters’ posterior distribution to describe the full range of the
model’s predicted flux. The errorbar for a particular datapoint represents the
width of the 16${}^{\textrm{th}}$ to 84${}^{\textrm{th}}$ percentile interval
for this distribution built up of random draws. Figures 8 and 9 show the
uncertainties of both models as errorbars on the models’ respective predicted
spectra using the same method. The $s-$factor dominating the uncertainty
results in errorbars that scale according to the magnitude of the flux
predicted by the model.

Figure 6: A comparison of the DEM fit using FUV lines and the X-ray spectrum
to a new DEM that excludes N V, C IV, Si IV from the FUV data and includes a
number of EUV lines listed in Table 1. As in Figures 3 and 4, the solid lines
and regions represent the median DEM and 1$\sigma$ confidence intervals while
the horizontal bars represent constraints imposed by the measured fluxes.

Figure 7: A comparison of the predicted line fluxes from the DEM models shown
in Figure 6 to the lines listed in Table 1. The color scheme is the same as
Figure 6, with the red pentagons marking predictions from the model including
EUV lines and the blue triangles marking predictions from the model excluding
EUV lines. The black points represent the line flux measurements with
errorbars for their measurement uncertainties while the model predictions have
errorbars for the $s$-factor model intrinsic uncertainty. The high S/N
assigned to the Woods et al. (2009) data makes the errorbars nearly invisible.
Note how the red model accurately predicts the fluxes of O III, O IV, O V
formed at roughly the same temperatures as N V, C IV, Si IV but drastically
underestimates the flux of these anomalous ions. Conversely, the blue model is
caught in a compromise that slightly underestimates the flux of these ions and
overestimates the flux of other FUV ions, but this compromise results in
significantly overestimating the flux of EUV ions.

Figure 8: The predicted X-ray spectra from both models discussed in Section
§4.3 compared to each other and the Woods et al. (2009) SIRS data used to fit
the DEM. The solid lines represent models according to the same color scheme
as Figure 6, blue for the model without EUV lines and red for the model with
EUV lines and without the anomalous FUV ions. The errorbars on each model
incorporate both the uncertainty in DEM shape and the $s-$factor uncertainty
of their respective models in the manner described in Section §4.3. The DEMs
do not differ significantly in the temperature regime associated with emitting
at X-ray wavelengths, so the spectra for both models are nearly identical. A
more complete look at the X-ray spectral data, including data not used to fit
the DEM, and our DEM prediction is shown in the top-left panel of Figure 16.

When we compare the predicted EUV flux of the Sun from both fits to the
observed spectrum itself, the DEM prediction fit without EUV lines
overestimates the data by 80% (see Figure 9 and Table 9). The $s-$factor
uncertainty for this model estimates that each line’s predicted flux has an
uncertainty $60\%$ of the predicted value, so the 1$\sigma$ confidence
interval of the model still encompasses the observed data. The DEM prediction
including EUV lines and excluding anomalous ions underestimates the integrated
flux by only $0.01\%$, but does have a significant $57\%$ $s-$factor
uncertainty. Within the EUV regime, there are 3 different recombination
continuum regions that form from excess kinetic energy emitted when an ion
captures a free electron into a bound state. Only one of them, the H I 912 Å
continuum, is a significant contributor to the total EUV flux integrated from
100 to 912 Å, accounting for 15% in the Woods et al. (2009) spectrum of the
Sun. The other two regions, He II 229 Å and He I 504 Å, contribute 3% and 2 %
respectively. This falls within our uncertainties on the predicted flux, but
it would be worth investigating if it is possible to incorporate these
recombination continua in the DEM model without added parameters. The
reconstruction of the EUV spectrum only used the emissivities of optically
thin emission lines in the $G_{\lambda}(T)$ emissivity matrices and does not
account for any contribution from continuum processes. The data required to
incorporate free-bound, free-free, and two-photon continua exist in CHIANTI,
but as of this work we have not included these sources of emissivity in the
$G_{\lambda}(T)$ emissivity matrices. In future work we hope to include these
processes for both fitting the DEM to X-ray spectra and reconstructing the EUV
spectrum.

Figure 9: The predicted EUV spectra of both DEM models discussed in Section
§4.3 with errorbars representing the uncertainties derived in the manner
described in that section. The color scheme follows that of Figures 6 and 8,
with the solid blue line representing the predicted spectrum and uncertainty
of the model fit without EUV lines and the solid red line representing the
predicted spectrum and uncertainty of the model fit with EUV lines and without
the anomalous ions. The black points are the data from Woods et al. (2009).
Both model spectra show nearly the same spectral shape as the data, excepting
the ramp feature at wavelengths between 800 and 912 Å, discussed in Section
§4.3. The blue model consistently overestimates the spectrum although the
1$\sigma$ errorbars usually include the data.

###  4.4 Abundance Sensitivity

We use the solar coronal abundances from Schmelz et al. (2012) stored in
CHIANTI and an ionization equilibrium calculated at those abundances to get
the $G_{ul}(T)$ functions for each emission line. Then for every $G_{ul}(T)$
function, we multiply by the stellar abundance if the ion’s atomic number
$Z>2$. A higher metal abundance should shift the ionization equilibrium, which
would have some effect on the self-consistency of the emissivity calculations.
To verify whether or not this effect would be significant, we tested how the
DEM fitting was sensitive to abundance in this crude framework. When we fit to
stars other than the Sun we can refer to the literature and use the best
abundance available, but we should not expect to always have an accurate and
precise abundance measurement for the stars we are fitting the DEM to. Fitting
the Sun with emissivity matrices generated using a super-solar [Fe/H] $=1$ and
a sub-solar [Fe/H] $=-1$ abundance resulted in the overall DEM shifting up or
down to compensate. The predicted EUV fluxes obtained by combining each model
DEM with their respective emissivity matrix differed by $<15\%$.

This makes sense since the $G_{ul}(T)$ of an emission line from any metal is
linearly proportional to the abundance. In fact, this harks back to the
original implementation of the emission measure distribution in Pottasch
(1963), where the author determined the relative abundances of elements by
shifting them around to minimize the scatter in the emission measure
distribution. This abundance adjustment is contingent on the assumption of the
relative abundances being constant throughout the optically thin plasma, but
everything else about this implementation of the DEM hinges on this assumption
as well. This will also let us compare the DEM shapes of stars with different
abundances by modifying the temperature independent coefficient $c_{0}$ to
shift the overall DEM shape, marginalizing over the stellar abundance for any
future study of the DEMs of a population of stars. However, if the relative
abundances of an M dwarf corona differ from the Sun, the ions of an affected
element will be consistently discrepant from the final DEM fit, and this is
likely our greatest source of systematic uncertainty for stars other than the
Sun. In the Sun, heavy elements with a first ionization potential $<9$ keV
tend to be enhanced in the corona relative to the photosphere (Drake et al.,
1995). There is evidence to suggest that this first ionization potential
effect varies with spectral type and may even be reversed in M dwarfs (Drake
et al., 1995; Wood et al., 2012; Laming, 2015).

A single line from a species constrains the DEM near its ionization
temperature, but many lines from an element across multiple ionization stages
can constrain the element’s relative abundance by fitting an element-specific
DEM and shifting it to match the DEM fit without that element. This will only
be possible for the brightest stars and the C, O, or Si species, since this
requires bright lines from at least three ionization stages formed at
conditions valid for the DEM method. Fitting DEMs for as many nearby main-
sequence stars as possible may allow us to generate this stellar coronal
abundance library. A future and intermediate improvement would be to use the
sample of stars for which _EUVE_ allowed papers like Drake et al. (1995) to
determine coronal abundances and select the most appropriate analog for a
target we are fitting the DEM to. We intend to use this approach for future
work.

###  4.5 Polynomial Degree

Higher polynomial degrees allow more flexibility between the well-constrained
temperatures. The spread of FUV and X-ray formation temperatures constrains
the slopes in their respective regimes ($2\times 10^{4}$ K $<$ FUV $<2\times
10^{5}$ K, $10^{6}$ K $<$ X-ray $<2\times 10^{7}$ K), and this rigidly
constrains lower order polynomials. The models show the most agreement with
each other around $T=10^{5}$ K, where there are a number of FUV lines with
overlapping emissivities to anchor the fit (see Figures 2 and 10), while there
are significant discrepancies at temperatures lower than the constraints
imposed by the FUV lines and higher than the constraints imposed by the X-ray
spectrum. Table 3 compares the BIC values for the different polynomial order
and since the Sun’s DEM shows very little curvature, the BIC prefers the
3${}^{\textrm{rd}}$ order model and penalizes the complexity of higher order
polynomials. Table 3 also shows that the models predict values consistent
within 1 $\sigma$ for the EUV flux integrated between 100 to 912 Å, and the
$s-$factor uncertainty decreases slightly for higher order polynomials. This
consistency is a product of both the low curvature of the Sun’s DEM and the
fact that the FUV and X-ray data constrain either end of the temperature
interval responsible for EUV emission lines. When we ran a similar test with
AU Mic, we found that orders below 5 were unable to match the complexity of
the data, so we move forward with the 5${}^{\textrm{th}}$ order as our
standard approach for fitting other stars. Most targets will have too few data
points to merit a model with many more parameters than a 5${}^{\textrm{th}}$
order polynomial.

Figure 10: This plot compares the Sun’s DEMs fit with different polynomial
orders to each other and the flux constraints imposed by the observed FUV
lines and X-ray spectrum used to fit the data according to the same scheme
described in the caption of Figure 3. The red model is a 3rd order polynomial,
the blue model is the 5${}^{\textrm{th}}$ order polynomial we adopt as our
standard approach, and the green is a 7th order polynomial. The spread in FUV
line formation temperatures sets a slope for the DEM to match between $2\times
10^{4}$ and $2\times 10^{5}$ K, while the X-ray spectrum sets the slope
between $10^{6}$ and $2\times 10^{7}$ K. Table 3 compares the different order
models to show that the different models predict consistent EUV fluxes. Table
3: The Bayesian Information Criterion (BIC) for each polynomial model.

Chebyshev Polynomial Order | BIC | Integrated EUV Flux | $\log_{10}s_{\mathcal{L}_{\mathrm{max}}}$ | $\log_{10}s$
---|---|---|---|---
— | — | [ergs s-1 cm-2] | — | —
$3$ | $-472.0$ | $3.8^{+2.6}_{-2.3}$ | $-0.2438$ | $-0.2182^{+0.0696}_{-0.0600}$
$4$ | $-468.4$ | $3.7^{+2.7}_{-2.3}$ | $-0.2516$ | $-0.2173^{+0.0690}_{-0.0610}$
$5$ | $-464.7$ | $3.5^{+3.1}_{-2.1}$ | $-0.2515$ | $-0.2177^{+0.0696}_{-0.0603}$
$6$ | $-464.7$ | $3.4^{+2.9}_{-2.1}$ | $-0.2372$ | $-0.2208^{+0.0692}_{-0.0597}$
$7$ | $-464.9$ | $3.0^{+2.4}_{-1.8}$ | $-0.2662$ | $-0.2348^{+0.0706}_{-0.0599}$
$8$ | $-461.1$ | $2.9^{+2.4}_{-1.8}$ | $-0.2852$ | $-0.2350^{+0.0698}_{-0.0600}$
$9$ | $-457.1$ | $2.9^{+2.5}_{-1.8}$ | $-0.2709$ | $-0.2337^{+0.0709}_{-0.0593}$
$10$ | $-454.0$ | $3.1^{+3.0}_{-1.9}$ | $-0.2858$ | $-0.2307^{+0.0717}_{-0.0616}$

Note. — The Sun’s DEM shows very little complexity, so the 3${}^{\textrm{rd}}$
order polynomial model is most preferred, with each subsequent order scoring a
worse BIC. This was not true for the AU Mic test for which orders below 5 were
too inflexible to match the constraints of the prior and the data.

##  5 Applying Our DEM Method to AU Mic

AU Mic is a nearby young and active M dwarf with an observable debris disk,
$9.979\pm 0.04$ pc away and $22\pm 3$ Myr old (MacGregor et al., 2013; Ibañez
Bustos et al., 2019; Plavchan et al., 2020). The system also hosts at least
one confirmed planet and possibly a second planet candidate (Plavchan et al.,
2020), making it a rare case of an M dwarf planetary system with a precisely
known young age. Despite being nearby, and having been observed multiple times
with _EUVE_ , the quiescent EUV spectrum of AU Mic is poorly constrained (see
Section 7 and Figure 16), so some method of reconstructing the EUV spectrum is
required to study this planetary system in detail. On the other hand, the
X-ray and FUV data for this star are extremely precise considering how
intrinsically faint the star is, allowing us to fit a very well-constrained
DEM and compare our implementation to earlier work published in Del Zanna et
al. (2002). Pagano et al. (2000) published a very thorough list of emission
line fluxes measured from spectra taken during quiescence, a subset of which
we use to fit the DEM for AU Mic and list in Table 4. Redfield et al. (2002)
published separate quiescent and flare FUSE measurements of FUV lines and
Redfield et al. (2003) reported the quiescent coronal line fluxes listed in
Table 4. The bandpass for _FUSE_ overlaps with that of COS, so we fit separate
DEMs including and excluding the _FUSE_ measurements to demonstrate the
usefulness of transitions observed between $900$ and $1100$ Å in constraining
the high-temperature end of the DEM, motivating future COS observations of
cool dwarfs.

The majority of the lines listed in Table 4 will be unobservable for other
fainter cool dwarfs, but having the ground truth of which lines are emitted
from the upper atmospheres of cool dwarfs is immensely useful for future DEM
fitting. If a line is observed for AU Mic, we can place upper limits on the
flux from that line for another star, constraining the DEM near that line’s
formation temperature. We fit the DEM assuming solar coronal abundances from
Schmelz et al. (2012), using the FUV lines listed in Table 4 and an X-ray
spectrum from the Reflection Grating Spectrometer (RGS, den Herder et al.
2001) on _XMM-Newton_ , observed in October 2018 (Kowalski et al., 2019). The
spectrum was resampled at 1 Å resolution before fitting. The lightcurve for
this observation showed multiple flares and we use the quiescent X-ray
spectrum from the work of Kowalski et al. (in prep) and Tristan et al. (in
prep). This observation was not concurrent with the data obtained by Pagano et
al. (2000) and Redfield et al. (2003), so it is possible that the X-ray and
FUV data are not from identical levels of quiescence. Very few targets will be
likely to have concurrent X-ray and FUV observations, so this issue will
plague most of the stars needing EUV reconstruction. We compare our fit using
FUV lines and an X-ray spectrum to another fit using the same FUV lines and
X-ray line flux measurements reported by Wood et al. (2018). We also note that
whenever possible, a line flux measurement is more useful than a spectral bin,
and if high temperature emission line strengths can be measured, these should
be favored over the use of spectral bins.

Figure 11: Comparing the DEM of AU Mic fit to a combination of FUV line fluxes
from STIS (published in Pagano et al. 2000) and _FUSE_ (published in Redfield
et al. 2002 and Redfield et al. 2003) and either a coarsely sampled X-ray
spectrum or a list of X-ray line fluxes reported in Wood et al. (2018). The
red line and region correspond to the median DEM value and 1$\sigma$
confidence intervals for AU Mic fit using only line flux measurements, while
the blue line and region are the model for AU Mic constrained by the X-ray
spectrum instead of line fluxes. The horizontal bars represent constraints
from the FUV lines listed in Table 4 and from the _XMM-Newton_ X-ray spectrum.
The dotted purple bars correspond to the STIS lines published in Pagano et al.
(2000) while the dashed red bars correspond to _FUSE_ lines published in
Redfield et al. (2003) and the black solid bars are from the quiescent _XMM-
Newton_ spectrum presented in Kowalski et al. (2019) (to be published in
Kowalski et al. in prep, and Tristan et al. in prep). There were too many
X-ray lines listed in Wood et al. (2018) to represent in this figure.

In Figure 11 we compare two fits for AU Mic to the constraints imposed by the
data, with the red model using only line flux measurements and the blue model
using the combination of FUV lines and an X-ray spectrum that will be
applicable to fainter M dwarfs. AU Mic has a higher DEM than the Sun across
the entire temperature domain, and significantly higher at temperatures
greater than $2\times 10^{6}$ K, corresponding to the corona. The presence of
detectable coronal iron lines, formed at $T>10^{6}$ K, places a strong
constraint that lifts AU Mic’s DEM far higher than the Sun’s which had many
X-ray spectrum points depressing the DEM in the vicinity of this temperature.
The iron emission lines are not solely responsible for the differences, as the
constraints from the X-ray spectra set a slope at these high temperatures
leading toward the iron lines. The shape of AU Mic’s DEM beyond $T=2\times
10^{5}$ K demonstrates the importance of including a corona in calculating the
total EUV flux. Peacock et al. (2019b) finds significant differences in the
total flux and spectral shape between the PHOENIX models without a corona and
the semi-empirical SRPM model from Fontenla et al. (2016). Peacock et al.
(2019b) also simulates the addition of a corona to their models by using the
DEM of AU Mic available in CHIANTI from Del Zanna et al. (2002), showing
potential opportunities for supplementing stellar atmosphere models with DEMs
fit to observations of specific stars.

Figure 12 compares the predicted X-ray spectra to the observed spectrum and
Figure 13 compares the predicted line fluxes to the FUV line profile
measurements listed in Table 4. Both figures incorporate the $s-$factor
uncertainties and DEM shape variation (Section §4.3) and show that the model
predictions are typically consistent with the data to 1$\sigma$. Figure 14
shows the parameter distributions for the DEM fit to line fluxes, and we find
that the $s-$factor for AU Mic is 0.4, comparable to that of the Sun. At lower
wavelengths the red model, which was not fit to the spectrum itself,
significantly underestimates the flux. This may be a consequence of the higher
energy emission including flux from free-free or free-bound continuum sources,
creating the discrepancies between both these DEMs. To reproduce the flux in
these bins without accounting for this extra emissivity, the DEM fit to the
spectrum must enhance the amount of material at these high temperatures.
Including free-free and free-bound continuum emissivities in the
$G_{\lambda}(T)$ matrices should mitigate or eliminate these discrepancies.

Figure 12: The X-ray data used to fit the DEM compared to the predicted
spectra for the DEM models of AU Mic shown in Figure 11. The black points
represent the downsampled RGS X-ray spectrum at a 1 Å wavelength resolution
while the red and blue lines show the DEM prediction for the flux density in
the same wavelength bins, incorporating the $s-$factor uncertainty in their
error bars, for the models fit with X-ray lines and the X-ray spectrum
respectively. They are consistent with each other, although the median
prediction of the DEM fit to the spectrum is consistently above the DEM fit to
X-ray line fluxes.

Figure 13: The FUV line fluxes predicted by the AU Mic DEM models shown in
Figure 11 compared to the data from Pagano et al. (2000) and Redfield et al.
(2003). Following the color scheme of Figures 11 and 12, the red triangles
represent the predictions from the DEM model fit using only line fluxes while
the blue pentagons represent the predictions from the DEM model using the
X-ray spectrum. The black points represent the combined fluxes of lines
emitted by the ion annotating the point, with errorbars for the measurement
uncertainty while the model predictions have errorbars representing the
$s-$factor model intrinsic uncertainty.

Table 4: Ion fluxes of AU Mic compared to the predictions from the red DEM
model shown in Figure 11.

Ion | Wavelengths | $\log_{10}T_{f}$ | Observed Flux | DEM Prediction
---|---|---|---|---
| [Å] | $\log_{10}(\textrm{[K]})$ | [$10^{-15}$ erg s-1 cm-2] | [$10^{-15}$ erg s-1 cm-2]
C II | 1324, 1336aamultiplet | 4.6 | $214\pm 16.0$ | $149^{+71.7}_{-96.2}$
C III | 977, 1176aamultiplet | 4.9 | $280\pm 16.0$ | $681^{+309}_{344}$
C IV | 1548, 1551 | 5.0 | $327\pm 24.2$ | $208^{+164}_{-174}$
N II | 1086 | 4.6 | $10.0\pm 2.00$ | $6.94^{+3.33}_{-4.21}$
N III | 990, 992 | 4.9 | $29.0\pm 2.24$ | $68.6^{+30.9}_{-33.3}$
N IV | 1487 | 5.1 | $0.90\pm 0.09$ | $0.747^{+0.337}_{-0.366}$
N V | 1238.8, 1242.8 | 5.3 | $66.3\pm 4.99$ | $41.6^{+19.1}_{-22.5}$
O III | 1666 | 4.9 | $1.30\pm 0.13$ | $1.00^{+0.448}_{-0.489}$
O IV | 1400, 1401, 1407 | 5.1 | $5.90\pm 0.423$ | $12.3^{+5.55}_{-6.17}$
O V | 1218, 1371 | 5.3 | $21.8\pm 1.86$ | $21.4^{+9.84}_{-12.2}$
O VI | 1032, 1038 | 5.5 | $315\pm 26.4$ | $183^{+89.7}_{-121}$
Ne V | 1146 | 5.4 | $2.00\pm 1.00$ | $1.30^{+0.599}_{-0.784}$
Si II | 1260, 1265, 1304, 1527, 1533 | 4.3 | $22.0\pm 1.06$ | $23.3^{+12.3}_{-20.7}$
Si III | 1108, 1110, 1113, 1206, 1295, 1297, 1299, 1301, 1303 | 4.8 | $101\pm 8.11$ | $624^{+289}_{-340}$
Si IV | 1394, 1403 | 4.9 | $78.5\pm 5.73$ | $126^{+57.3}_{-62.8}$
S II | 1253.8, 1259.5 | 4.4 | $2.60\pm 0.18$ | $2.78^{+1.45}_{-2.16}$
S III | 1016, 1021, | 4.8 | $7.00\pm 2.24$ | $7.11^{+3.29}_{-3.75}$
S IV | 1063, 1073 | 5.0 | $6.00\pm 1.41$ | $28.2^{+12.6}_{-13.7}$
Fe XVIII | 975 | 6.9 | $6.00\pm 2.00$ | $8.75^{+4.10}_{-5.08}$
Fe XIX | 1118 | 7.0 | $7.00\pm 1.00$ | $4.51^{+2.05}_{-2.42}$
Fe XXI | 1354 | 7.0 | $7.70\pm 0.77$ | $3.73^{+1.91}_{-2.93}$

Note. — All entries are from STIS line fluxes reported in Pagano et al. (2000)
and/or _FUSE_ measurements reported in Redfield et al. (2002) and Redfield et
al. (2003).

Figure 14: The corner plot showing the parameter distributions when fitting AU
Mic with a 5${}^{\textrm{th}}$ order polynomial to FUV lines and the X-ray
spectrum. The model uncertainty for AU Mic is slightly better than that of the
Sun, with $s=0.4$ instead of 0.6.

##  6 Our DEMs Compared to Published Literature

Differential emission measure techniques are primarily applied in the solar
context to resolved regions, highlighting individual structures like coronal
holes or flare loops. Vernazza & Reeves (1978) fit DEMs to many such
structures, including “quiet regions” with minimal observed activity during
the 9-month interval from 1973 May to 1974 August, shortly before the solar
minimum of Cycle 21. This is not perfectly analogous to our quiescent Sun DEM
fit to data integrated over the entire solar disk, but it is the best
comparison for which we had access to a published DEM via CHIANTI (Dere et
al., 1997; Del Zanna et al., 2015). The left panel of Figure 15 compares the
Vernazza & Reeves (1978) quiet region DEM to our disk-integrated quiescent Sun
DEM, fit under the assumptions applicable to observing the Sun as a star
(albeit with much higher signal-to-noise): the few FUV lines listed in Table 1
and the X-ray spectrum at a low resolution of $R\leq 50$. We also include the
solar DEM fit excluding anomalous ions and including EUV lines.

Similarly, the right panel of Figure 15 compares our DEM for AU Mic to the DEM
published by Del Zanna et al. (2002), which combined _FUSE_ , STIS, and _EUVE_
observations. CHIANTI has a volume differential emission measure version of
the Del Zanna et al. (2002) DEM which needed to be divided by $\frac{4\pi
R_{\star}^{2}}{d^{2}}$ to match our formulation of a column differential
emission measure, and we use the stellar radius and distance assumed by Del
Zanna et al. (2002), $R_{\star}=0.68R_{\odot}$ and $d=9.94$ pc, for
consistency in the scaling factor. Like our dataset, they did not have access
to simultaneous observations from these different instruments. Unlike our
dataset, they included the integrated fluxes of lines observed with _EUVE_ ,
but we believe the _EUVE_ observations they used were flare contaminated for
reasons discussed later in Section §7. We also compare our AU Mic DEM to the
active solar region DEM of Vernazza & Reeves (1978), demonstrating a small
resemblance between the active region and the active star AU Mic. The active
region DEM seems to shift the shape of the quiet sun DEM to a higher
temperature and dramatically enhance the DEM near $10^{6}$ K. AU Mic, a star
more active than the Sun, has more material at nearly all temperatures.

Figure 15: The left panel compares both of our DEMs for the quiescent disk-
integrated Sun described in Section §4.3 to the “quiet Sun” DEM derived by
Vernazza & Reeves (1978) from observations of quiet regions of the Sun, while
the right panel compares our DEM for AU Mic to the AU Mic DEM published by Del
Zanna et al. (2002) and the DEM of an active region published by Vernazza &
Reeves (1978). These DEMs used for comparison were available in the CHIANTI
database (Dere et al., 1997; Del Zanna et al., 2015).

Both panels show that our DEMs are significantly higher than their literature
counterparts at temperatures near $10^{5}$ K, roughly corresponding to the
transition region. For the AU Mic panel on the right, the Del Zanna et al.
(2002) DEM is significantly lower than ours between $3\times 10^{4}$ and
$3\times 10^{5}$ K. Monsignori Fossi et al. (1996) published the time-
evolution of the AU Mic DEM during a flare observed by _EUVE_ in July 1992,
and while they use a slightly different formulation of the DEM
$=n_{e}^{2}\frac{dV}{dT}$, the shape of our DEM strongly resembles theirs
published in panels a and g of their Figure 9, which correspond to quiescent
phases, while multiplying our DEM by $4\pi R_{\star}^{2}\approx 3\times
10^{22}$ cm2 scales to the approximately the same order of magnitude as their
DEM. Unfortunately, their DEM was not available in CHIANTI for direct
comparison.

The discrepancy at temperatures below $10^{6}$ K is largely driven by the line
list we are using to fit the DEM, where our strongest lines are from Na-like
and Li-like species N V, Si IV, and C IV, and the solar DEM excluding these
ions is much closer to the Vernazza & Reeves (1978) quiet region DEM. Del
Zanna et al. (2002) fits the DEM of AU Mic without these lines and has enough
individual lines to constrain the DEM in this temperature regime without them.
We will not be able to afford this luxury for nearly every other M dwarf
unless we are observing a strong flare. Instead we fit our DEM including these
lines and use the $s-$parameter boost to our variance to account for the
systematic uncertainties involved. As long as some other ions formed in the
same region are included in the line list, even with just upper limits on
their fluxes, the DEM will shift down to accommodate these lines. In the
future, we plan to test dividing the observed fluxes of these lines by a
corrective factor $\sim 5$ before fitting the DEM to see how this improves the
fit and affects the predicted EUV spectrum.

Working in a data-limited regime is also why we fit for the DEM using an
assumed functional form instead of interpolating between the emission measure
loci estimated from individual lines (Pagano et al., 2000), or fitting for the
value of the DEM in discrete temperature bins (Del Zanna et al., 2002, 2015).
Without measured lines in the temperature regime corresponding to most EUV
lines, we must use our assumption of a continuous function anchored on both
ends of the inaccessible temperature/wavelength regime. The DEM for AU Mic
derived by Del Zanna et al. (2002) is poorly constrained between
$10^{6}<T<10^{6.7}K$ because they do not use X-ray line fluxes or a coarse
spectrum, but our DEMs agree on the position and magnitude of the coronal peak
DEM ($T=10^{7}$ K), if not the shape of the decline. Other contributions to
the discrepancies between DEMs are differences in the atomic data and
calculations for abundances and the ionization equilibrium and level
populations. While the DEMs created by Sanz-Forcada et al. (2011) were not
available in CHIANTI, we expect significant discrepancies at the lower
temperature end of the DEM because the majority of their cool dwarfs did not
have UV data available.

##  7 Comparing Model EUV Spectra to Data

With our implementation of the DEM well-characterized, we can move on to the
main objective of this project: generating EUV spectra in a format useful to
the astronomical community with errorbars that self-consistently account for
both statistical and systematic uncertainties. Our model spectra range from 1
to 2000 Å at a constant $R=\frac{\Delta\lambda}{\lambda}=500$, but we advise
using data instead of our model in the regimes where that is possible. The
model spectrum files include uncertainties derived according to the method
described in Section §4.3. Figure 16 compares our model spectrum of the Sun to
the Woods et al. (2009) data and our model spectrum of AU Mic to an _EUVE_
observation from July 1992 and quiescent _FUSE_ data from a different time
described in Redfield et al. (2002, 2003). This _EUVE_ observation was during
the calibration phase of the mission and happened to catch a flare, first
reported by Cully et al. (1993) and later studied in more detail in Monsignori
Fossi et al. (1996). We compare our DEM-generated spectrum to data from the
quiescent time segment before the flare, extracted using standard _EUVE_ Guest
Observer Center IRAF procedures. Useful data were obtained with the short
wavelength (70 to 190 Å) and medium wavelength (140 to 380 Å) detectors, but
no useful signal was present in the long wavelength (280 to 760 Å) region. The
spectral resolution of the 3 spectrometers is $\sim$ 0.5, 1.0, and 2 Å , which
corresponds to 7 pixels per resolution element on the detectors. The photon
event data were screened to eliminate high background times and times when the
detectors were switched off. The IRAF routine “apall” was used for the
spectral extraction with a 14 pixel wide spectral region and two 85 pixel wide
background regions measured above and below the stellar spectrum. The count
rate spectra were converted to flux densities using the effective areas
established by the _EUVE_ mission and the oversampled spectra were smoothed to
the intrinsic spectrometer resolution. The large background area sampled
allowed precise monitoring of the time dependent background. Del Zanna et al.
(2002) used a 1993 _EUVE_ observation which has no published lightcurve to
verify the absence of flares and the time-averaged 1993 spectrum shows more
flux than the quiescent 1992 data at all wavelengths indicating possible flare
contamination.

Figure 16: The individual rows of this figure show 250 Å chunks of the X-ray,
extreme ultraviolet, and a portion of the far ultraviolet regimes. The left
side compares to the Woods et al. (2009) spectrum shown as black points to the
DEM-generated model spectrum for the Sun in blue, where the DEM was fit to EUV
lines, FUV lines excluding anomalous ions, and X-ray data. The right side
shows the DEM-generated model spectrum for AU Mic in red, compared to both
_EUVE_ data shown as black points and _FUSE_ data shown as purple squares. The
model spectrum for the Sun is largely consistent with the data, except for a
recombination continuum associated with the H I 912 Å line spanning $\sim 800$
to $912$ Å. The AU Mic _EUVE_ data is almost entirely consistent with zero,
but our model does not predict flux significantly higher than these upper
limits, while roughly matching some of the clearly detected emission lines.

Our DEM prediction for the EUV luminosity of AU Mic is
$L_{\mathrm{EUV}}=9.0^{+6.1}_{-4.6}\times 10^{28}$ erg s-1, while the Linsky
et al. (2014) relations give $1.2\times 10^{29}$ erg s-1, Sanz-Forcada et al.
(2011) relation gives $1.4\times 10^{30}$ erg s-1, and Chadney et al. (2015)
finds $L_{\mathrm{EUV}}=8.4\times 10^{28}$ erg s-1 using a coronal emission
measure distribution model. We use the Lyman-$\alpha$ flux reported by Wood et
al. (2005) for the Linsky et al. (2014) relations, and the X-ray luminosity
obtained by multiplying the integrated flux of our EPIC-MOS spectrum with
$4\pi d^{2}$, $L_{\textrm{X-ray}}=3.0\times 10^{29}$ erg s-1, for the Sanz-
Forcada et al. (2011) relation. To convert the surface EUV flux reported by
Chadney et al. (2015) to luminosity we use their assumed radius for AU Mic
$=0.68R_{\odot}$. The major advantages of our DEM approach are a well-
characterized uncertainty and a balance between ease of implementation and
specificity to each star. Using Equation 3 of Sanz-Forcada et al. (2011) has a
minimum uncertainty of 1.99 dex in predicting the EUV flux, but our method can
do significantly better than this empirical relation even for very faint
targets like TRAPPIST-1, with only a few FUV line measurements and an
integrated X-ray flux or coarse spectrum. The DEM method also provides a low-
resolution spectral shape in addition to a total flux, which may be useful for
those who wish to model more detailed effects of high energy stellar radiation
on a planet, disk, or the local interstellar medium.

##  8 Case Studies

Thus far, we have demonstrated our method on targets with extremely good data,
which are not representative of the majority of stars for which the
astronomical community needs reconstructed EUV spectra to enable other
science. In this section we apply our method to three M dwarfs of interest to
the exoplanet community, all fainter and less active than AU Mic, and show
that their coronae are hotter than the Sun and will seriously affect the
atmospheric evolution of planets orbiting in their respective habitable zones.

###  8.1 GJ 832

GJ 832 was included in the original MUSCLES survey (France et al., 2016) with
its EUV flux estimated by the Linsky et al. (2014) correlations and semi-
empiricaly modeled by both Fontenla et al. (2016) and Peacock et al. (2019b).
We fit our DEM model to the line fluxes published in Youngblood et al. (2016)
and X-ray data hosted on _MAST_ as part of the MUSCLES data products. The line
fluxes were measured using STIS data while the X-ray spectrum is from _XMM-
Newton_ EPIC (France et al., 2016). Comparing our model spectra to the semi-
empirical models and the EUV fluxes predicted by the Linsky relations in
Figure 18 shows that the different models agree with each other in different
wavelength regimes. The data we used to fit GJ 832 was not perfectly
quiescent, so some of the model discrepancy may be due to flare contamination
(Youngblood et al., 2016). The DEM inferred from fitting the data is shown in
Figure 17.

Figure 17: Comparing the DEM of GJ 832 fit to a quiescent X-ray spectrum from
_XMM-Newton_ EPIC and line fluxes from STIS, with $\overline{\Psi}(T)$
constraints shown as horizontal dotted blue and solid black bars respectively.

The DEM model tends to predict higher fluxes than the other three models at
wavelengths shorter than 600 Å. Below 400 Å all models except for the Peacock
et al. (2019b) PHOENIX model coincide quite closely, with the outlier lacking
a coronal contribution. Between 400 to 600 Å both the DEM and semi-empirical
models predict higher fluxes than the Lyman-$\alpha$ correlations. The biggest
discrepancies between models are between 800 to 1100 Å where the semi-
empirical models include the recombination continuum at the H I 912 Å line and
the blue wing of the Lyman-$\alpha$ 1216 Å line, both of which are unaccounted
for by the DEM model and which the Lyman-$\alpha$ correlation fluxes seem to
underestimate. Tilipman et al. (2020) updates the Fontenla et al. (2016) SSRPM
model of GJ 832, and when all these EUV reconstruction methods have been
applied to a larger sample of stars we may have more insight into the
conditions under which each is more likely to be accurate.

Figure 18: Comparing four methods of reconstructing the EUV spectrum of GJ
832. The green line and errorbars represent the spectrum predicted by our DEM
model with the errors enclosing the $16^{\textrm{th}}$ to $84^{\textrm{th}}$
percentile intervals for the predicted flux according to the method described
in Section §4.3. The red lines show the predicted EUV flux in 100 Å bandpasses
according to the Lyman-$\alpha$ correlations published in Linsky et al.
(2014), available in the MUSCLES dataset. The blue points represent the semi-
empirical model published by Fontenla et al. (2016), also available in the
MUSCLES dataset. The purple triangles show the PHOENIX model published in
Peacock et al. (2019b) and hosted on _MAST_. These four methods of predicting
the EUV spectrum appear to have different regions of mutual agreement,
discussed in Section 8.1.

###  8.2 Barnard’s Star

Barnard’s Star is old and inactive compared to most M dwarfs (Ribas et al.,
2018), but it still flares occasionally (Paulson et al., 2006). France et al.
(2020) obtained X-ray (_Chandra_ ACIS-S) and FUV (_HST_ STIS/COS) data of this
star both during quiescence and during a flare. We used these data to prepare
both a quiescent DEM and a flare DEM, under the assumptions that the system is
still in a collisionally dominated equilibrium and with solar coronal
abundances scaled by the stellar metallicity [Fe/H] $=-0.32$. This is a very
low S/N regime in the quiescent data but not nearly as low as TRAPPIST-1,
discussed in Section 8.3. A more physically accurate flare DEM would require
adjusting the emissivity matrix to account for the magnetic reconnection’s
influence on the level populations. In the current framework, we see that the
flaring state has more material than quiescence at temperatures between $\sim
10^{5}$ to $10^{7}$ K, but both the flare and quiescent DEMs agree beyond
$3\times 10^{7}$ K (see Figure 19). France et al. (2020) uses our quiescent
and flaring model EUV spectra to investigate the influence of EUV variability
on a hypothetical planet orbiting in the habitable zone of Barnard’s star,
demonstrating the applicability of our low-resolution spectra to models of
atmospheric escape more complicated than simple energy-limited
photoevaporation. We note that this DEM shape is very unusual compared to past
published DEMs, but emphasize that the wide error intervals are likely to
encompass the true DEM shape. If more data become available, perhaps the model
constraints will narrow to something more familiar.

Figure 19: The left panel compares the DEM of Barnard’s Star fit to a coarse
X-ray spectrum from _Chandra_ ACIS-S and line fluxes measured from _HST_ STIS.
The right panel compares the DEM of Barnard’s Star fit to an integrated X-ray
flux from _Chandra_ ACIS-S and line fluxes from STIS/COS during a flare. For
both panels, the dotted blue and solid black bars represent the
$\overline{\Psi}(T)$ constraints for the X-ray and FUV data respectively.

###  8.3 TRAPPIST-1

TRAPPIST-1 is an ultracool dwarf with a gaggle of seven planets discovered
through transits (Gillon et al., 2017). It is included in the Mega-MUSCLES
survey, an extension of the original MUSCLES survey from France et al. (2016),
and the measured FUV line fluxes (_HST_ STIS) and X-ray spectrum (_XMM-Newton_
EPIC) will be published shortly in Wilson et al. (submitted), where we discuss
the implementation of our DEM method for this specific target. Wilson et al.
(submitted) also compares the MegaMUSCLES spectral energy distribution of
TRAPPIST-1 to the PHOENIX model published in Peacock et al. (2019a).
TRAPPIST-1 tests the fitting in a very low S/N regime constrained by a few FUV
lines and a very faint X-ray spectrum, but we still get meaningful fits and
constraints on the EUV flux (see Figures 20 and 22).

Figure 20: Comparing the DEM of TRAPPIST-1 fit to a quiescent X-ray spectrum
from _XMM-Newton_ EPIC and line fluxes from STIS, with $\overline{\Psi}(T)$
constraints shown as horizontal dotted blue and solid black bars respectively.

###  8.4 Comparing the Entire Sample

Figure 21 shows the DEMs of all the stars considered within this work and
shows some preliminary trends with activity and spectral type: more active
stars have a higher mean DEM while the decline in the DEM associated with the
move from the transition region to the corona appears to shift to higher
temperatures at cooler spectral types. With a more complete sample, we could
go a step further to interpolate the EUV flux of cool dwarfs that lack
observed FUV and X-ray data by relating the DEM to the more accessible stellar
parameters $T_{\textrm{eff}}$ and age, as traced by rotation and/or activity
indicators from optical spectra. Determining the best method for this
interpolation and testing its accuracy and precision is left to future work.

Figure 21: The DEMs of all stars considered in this work are plotted here with
the solid line representing the median predicted value of $\Psi(T)$ and shaded
intervals enclosing the 16${}^{\textrm{th}}$ to 84${}^{\textrm{th}}$
percentile intervals. Filling out this plot with more stars across a grid of
effective temperature and stellar activity will allow us to investigate trends
in the DEM and EUV spectra of cool dwarfs across the entire population.

We summarize the data products of this paper by plotting all of our DEM-
generated EUV spectra scaled to the flux density observed at a distance of 1
AU from the host star in Figure 22. The slope of the spectra across the EUV
seems to vary as a function of both spectral type and activity. The hotter
stars seem to have more EUV flux between 800 to 900 Å while the more active
stars have more EUV flux between 100 to 600 Å. These wavelength regions
correspond to lines formed roughly at temperatures $2\times 10^{5}$ and
$3\times 10^{6}$ K respectively (see Figure 1). AU Mic, which is one of the
more active stars and also one of the hotter stars in this sample, has a
roughly flat EUV spectrum. TRAPPIST-1, which is also very active but much
cooler, shows a strong negative slope from 100 to 1000 Å. The spectral shape
of the EUV is controlled by the relative strengths of the corona and
chromosphere, and these preliminary observations of our EUV spectra
conceptually agree with the findings of Linsky et al. (2020). Linsky et al.
(2020) measured the relationship between X-ray and Lyman-$\alpha$ flux for a
large sample of FGKM dwarfs and found that for older and relatively inactive
stars, the inverse relationship between coronal emission and effective
temperature is much stronger than the inverse relationship between
chromospheric emission and effective temperature. Trends in the shape of the
EUV spectrum generated by the DEM should be investigated along with trends in
the DEM and the data used to inform the fitting, and our current sample is
simply too small to make stronger claims than these extremely tentative
observations. Furthermore, the overall shape of the EUV spectrum may be
significantly altered by the inclusion of continuum processes and an analysis
of trends without this source of emissivity is premature. Table 9 lists the
integrated EUV flux at 1 AU, median $s-$factor, effective temperature, stellar
radius, and distance for each star considered in this work.

Figure 22: A comparison of all of our model spectra generated for the stars
considered in this work, downsampled from a resolution of
$R=\frac{\Delta\lambda}{\lambda}=500$ to 100 Å wavelength bins, scaled to the
flux received at 1 AU. The distance and radius assumed for each star is listed
in Table 9, and can be scaled to surface flux or habitable zone distances by
multiplying with the appropriate factor. On the left panel, we have the hotter
stars of the sample: the Sun (spectrum generated from the DEM model without
anomalous ions discussed in Section §4.3) in blue, AU Mic in red, and GJ 832
in green. On the right panel we have Barnard’s Star in quiescence represented
by the brown line, Barnard’s Star while flaring shown in pink, and TRAPPIST-1
in dark purple. All the spectra shown have errorbars determined by the method
described in Section §4.3.

##  9 Conclusions and Future Work

Our tests with the Sun show that with a low-resolution X-ray spectrum and the
strongest FUV emission lines, we can predict the EUV flux to within a factor
of 2 across the entire EUV region. Furthermore, our characterization of the
uncertainties in our method show that our predicted EUV spectra are consistent
with the data for the Sun, with meaningful error bars to support that
statement. We have demonstrated that our choice of functional form is able to
describe the DEM within the temperature region relevant to predicting EUV
flux, and that uncertainties in the abundance and average electron pressure
can be accounted for and propagated to the final output spectra. While our
approach to fitting the DEM has serious limitations discussed in Sections §4.3
and §6, our tradeoff exchanging precision for simplicity allows us to fit DEMs
and estimate EUV spectra for many more stars than more sophisticated methods
which require EUV data or laborious iterations, both manual and computational.
The method can be refined in the future to handle different stellar abundances
more carefully, to use temperature-pressure profiles from stellar atmosphere
models incorporating more physics, and to find a better way to handle the
discrepant CHIANTI ionization equilibrium of the Na-like and Li-like
isoelectronic sequences. Including the free-free and free-bound emission of
hydrogen and helium species is possible with existing CHIANTI data and these
emissivity sources will be accounted for in future DEM fits and EUV spectral
reconstructions. The greatest problem with our current approach is the lack of
elemental abundances tailored to individual stars based on their effective
temperature and age. Updating the abundances we use and refitting the DEM may
dramatically improve our precision, but this is left to future work.

_Hubble_ is the only observatory with FUV spectroscopy and the capability to
observe a large sample of M dwarfs, but its lifetime is limited and there will
be no replacement in the next few decades. Given the interest in M dwarf
planetary systems, a number of survey programs have proposed using _Hubble_ to
build up a spectral atlas of M dwarfs across the broad range of effective
temperature and activity level represented within the spectral type. MUSCLES
(France et al., 2016), HAZMAT (Shkolnik & Barman, 2014), FUMES (Pineda et al.,
in prep), and Mega-MUSCLES (Froning et al., 2019) are completed programs with
available data while observations for MEATS, a _Hubble_ survey targeting cool
dwarf exoplanet hosts scheduled to be observed by _Webb_ (HST-GO-16166, PI-
France), are forthcoming. Between all these surveys and archival data, we will
have enough FUV and X-ray data of cool dwarfs to start characterizing them as
a population: fitting DEMs to all nearby cool dwarf stars with sufficiently
available data and estimating their EUV flux. With a sufficiently
comprehensive DEM library, it may be possible to interpolate DEMs for stars
too faint for FUV or X-ray observations, and build on the work of Sanz-Forcada
et al. (2011) to calculate the EUV luminosities of all known planet-hosting
main-sequence stars. Until stellar atmosphere models become sophisticated
enough to have a grid of models varying both effective temperature and stellar
activity across the entire cool dwarf regime, or an EUV observatory (France et
al., 2019, e.g.) is able to provide directly observed EUV spectra of nearby
stars, differential emission measure techniques can satisfy the need for
stellar EUV spectra which are physically informed and empirically calibrated.

Star | Integrated 100 - 912 Å EUV Flux at 1 AU | Median $s$ | Star | Effective Temperature | Stellar Radius | Distance
---|---|---|---|---|---|---
— | [erg s-1 cm-2] | — | — | [K] | [$R_{\odot}$] | [pc]
Sun (Woods et al. 2009) | $1.99$ | — |  | 5772 [1][1]footnotemark:  | 1 | $4.848\times 10^{-6}$
Sun (DEM fit to FUV + X-ray)aaThis integrated flux is from the DEM fit to the Sun using the X-ray spectrum and FUV lines. | $3.54^{+3.12}_{-2.12}$ | $0.60$ |  | 5772[1][1]footnotemark:  | 1 | $4.848\times 10^{-6}$
Sun (DEM fit without anomalous ions)bbThis integrated flux is from the DEM fit to the Sun using the X-ray spectrum, EUV lines, and FUV lines excluding the anomalous ions N V, C IV, and Si IV. | $1.99^{+1.28}_{-1.15}$ | $0.57$ |  | 5772[1][1]footnotemark:  | 1 | $4.848\times 10^{-6}$
AU Mic | $31.9^{+21.8}_{-16.2}$ | $0.42$ |  | 3700[2][2]footnotemark:  | 0.75[2][2]footnotemark:  | 9.979[2][2]footnotemark:
GJ 832 | $1.66^{+1.30}_{-1.06}$ | $0.62$ |  | 3657[3][3]footnotemark:  | 0.499[4][4]footnotemark:  | 4.965[5][5]footnotemark:
TRAPPIST-1 | $0.762^{+1.30}_{-0.744}$ | $0.74$ |  | 2516[6][6]footnotemark:  | 0.121[6][6]footnotemark:  | 12.43[6][6]footnotemark:
Barnard’s Star (Quiescent) | $0.0183^{+0.109}_{-0.00860}$ | $0.43$ |  | 3278[7][7]footnotemark:  | 0.178[7][7]footnotemark:  | 1.83[5][5]footnotemark:
Barnard’s Star (Flaring) | $0.146^{+0.112}_{-0.0965}$ | $0.45$ |  | 3278[7][7]footnotemark:  | 0.178[7][7]footnotemark:  | 1.83[5][5]footnotemark:

Table 5: The DEM predictions for the integrated EUV flux at 1 AU and median
$s-$factor uncertainty for each star considered in this work.

Note. — To enable scaling the spectrum to other quantities, we also list the
stellar effective temperature, radius, and distance assumed in this work.

References. — [1] Mamajek et al. (2015), [2] Plavchan et al. (2020), [3]
Bailey et al. (2009), [4] Houdebine (2010), [5] Gaia Collaboration (2018), [6]
Van Grootel et al. (2018), [7] Ribas et al. (2018)

The authors thank the the referee for their constructive comments which
improved the quality of this paper. G.M.D. acknowledges the contributions of
conversations with Steve Cranmer in the initial stages of this project. This
work was funded in part by grants HST-GO-14640, HST-GO-15071, and HST-
GO-15264, which were provided by NASA through a grant from the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Incorporated, under NASA contract NAS5-26555. ZKBT
acknowledges funding for this work by the National Science Foundation under
Grant No. 1945633. This research is based on observations made with the _Far
Ultraviolet Spectroscopic Explorer_ and _Extreme Ultraviolet Explorer_ ,
obtained from the MAST data archive at the Space Telescope Science Institute,
which is operated by the Association of Universities for Research in
Astronomy, Inc., under NASA contract NAS 5–26555. We also used observations
obtained with _XMM-Newton_ , an ESA science mission with instruments and
contributions directly funded by ESA Member States and NASA. The analysis of
these observations was supported by NASA/_XMM-Newton_ AO-17 grant number
82274, ”A Unified Understanding of Flare Heating”.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Ayres (2014) Ayres, T. R. 2014, AJ, 147, 59, doi: 10.1088/0004-6256/147/3/59
  * Bailey et al. (2009) Bailey, J., Butler, R. P., Tinney, C. G., et al. 2009, ApJ, 690, 743, doi: 10.1088/0004-637X/690/1/743
  * Ballard (2019) Ballard, S. 2019, AJ, 157, 113, doi: 10.3847/1538-3881/aaf477
  * Barclay et al. (2018) Barclay, T., Pepper, J., & Quintana, E. V. 2018, ApJS, 239, 2, doi: 10.3847/1538-4365/aae3e9
  * Baron & Hauschildt (2007) Baron, E., & Hauschildt, P. H. 2007, A&A, 468, 255, doi: 10.1051/0004-6361:20066755
  * Chadney et al. (2015) Chadney, J. M., Galand, M., Unruh, Y. C., Koskinen, T. T., & Sanz-Forcada, J. 2015, Icarus, 250, 357, doi: 10.1016/j.icarus.2014.12.012
  * Claire et al. (2012) Claire, M. W., Sheets, J., Cohen, M., et al. 2012, ApJ, 757, 95, doi: 10.1088/0004-637X/757/1/95
  * Cox & Reynolds (1987) Cox, D. P., & Reynolds, R. J. 1987, ARA&A, 25, 303, doi: 10.1146/annurev.aa.25.090187.001511
  * Craig et al. (1997) Craig, N., Abbott, M., Finley, D., et al. 1997, ApJS, 113, 131, doi: 10.1086/313052
  * Cully et al. (1993) Cully, S. L., Siegmund, O. H. W., Vedder, P. W., & Vallerga, J. V. 1993, ApJ, 414, L49, doi: 10.1086/186993
  * Del Zanna et al. (2015) Del Zanna, G., Dere, K. P., Young, P. R., Landi, E., & Mason, H. E. 2015, A&A, 582, A56, doi: 10.1051/0004-6361/201526827
  * Del Zanna et al. (2002) Del Zanna, G., Landini, M., & Mason, H. E. 2002, A&A, 385, 968, doi: 10.1051/0004-6361:20020164
  * den Herder et al. (2001) den Herder, J. W., Brinkman, A. C., Kahn, S. M., et al. 2001, A&A, 365, L7, doi: 10.1051/0004-6361:20000058
  * Dere et al. (1997) Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149, doi: 10.1051/aas:1997368
  * Drake et al. (1995) Drake, J. J., Laming, J. M., & Widing, K. G. 1995, ApJ, 443, 393, doi: 10.1086/175533
  * Dressing & Charbonneau (2015) Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45, doi: 10.1088/0004-637X/807/1/45
  * Fontenla et al. (2016) Fontenla, J. M., Linsky, J. L., Witbrod, J., et al. 2016, ApJ, 830, 154, doi: 10.3847/0004-637X/830/2/154
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * France et al. (2018) France, K., Arulanantham, N., Fossati, L., et al. 2018, ApJS, 239, 16, doi: 10.3847/1538-4365/aae1a3
  * France et al. (2016) France, K., Loyd, R. O. P., Youngblood, A., et al. 2016, ApJ, 820, 89, doi: 10.3847/0004-637X/820/2/89
  * France et al. (2019) France, K., Fleming, B. T., Drake, J. J., et al. 2019, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11118, Proc. SPIE, 1111808, doi: 10.1117/12.2526859
  * France et al. (2020) France, K., Duvvuri, G., Egan, H., et al. 2020, arXiv e-prints, arXiv:2009.01259.  https://arxiv.org/abs/2009.01259
  * Froning et al. (2019) Froning, C. S., Kowalski, A., France, K., et al. 2019, ApJ, 871, L26, doi: 10.3847/2041-8213/aaffcd
  * Gaia Collaboration (2018) Gaia Collaboration. 2018, VizieR Online Data Catalog, I/345
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456, doi: 10.1038/nature21360
  * Hauschildt (1993) Hauschildt, P. H. 1993, J. Quant. Spec. Radiat. Transf., 50, 301, doi: 10.1016/0022-4073(93)90080-2
  * Hauschildt & Baron (2006) Hauschildt, P. H., & Baron, E. 2006, A&A, 451, 273, doi: 10.1051/0004-6361:20053846
  * Hawley (1993) Hawley, S. L. 1993, PASP, 105, 955, doi: 10.1086/133262
  * Henry et al. (2006) Henry, T. J., Jao, W.-C., Subasavage, J. P., et al. 2006, AJ, 132, 2360, doi: 10.1086/508233
  * Houdebine (2010) Houdebine, E. R. 2010, MNRAS, 407, 1657, doi: 10.1111/j.1365-2966.2010.16827.x
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ibañez Bustos et al. (2019) Ibañez Bustos, R. V., Buccino, A. P., Flores, M., et al. 2019, MNRAS, 483, 1159, doi: 10.1093/mnras/sty3147
  * Jones & West (2016) Jones, D. O., & West, A. A. 2016, ApJ, 817, 1, doi: 10.3847/0004-637X/817/1/1
  * Kass & Raftery (1995) Kass, R. E., & Raftery, A. E. 1995, Journal of the American Statistical Association, 90, 773.  http://www.jstor.org/stable/2291091
  * Kowalski et al. (2009) Kowalski, A. F., Hawley, S. L., Hilton, E. J., et al. 2009, AJ, 138, 633, doi: 10.1088/0004-6256/138/2/633
  * Kowalski et al. (2019) Kowalski, A. F., Allred, J., Axelson, R., et al. 2019, in American Astronomical Society Meeting Abstracts, Vol. 233, American Astronomical Society Meeting Abstracts #233, 360.14
  * Laming (2015) Laming, J. M. 2015, Living Reviews in Solar Physics, 12, 2, doi: 10.1007/lrsp-2015-2
  * Linsky et al. (2014) Linsky, J. L., Fontenla, J., & France, K. 2014, ApJ, 780, 61, doi: 10.1088/0004-637X/780/1/61
  * Linsky et al. (2020) Linsky, J. L., Wood, B. E., Youngblood, A., et al. 2020, ApJ, 902, 3, doi: 10.3847/1538-4357/abb36f
  * Louden et al. (2017) Louden, T., Wheatley, P. J., & Briggs, K. 2017, MNRAS, 464, 2396, doi: 10.1093/mnras/stw2421
  * Loyd et al. (2018a) Loyd, R. O. P., Shkolnik, E. L., Schneider, A. C., et al. 2018a, ApJ, 867, 70, doi: 10.3847/1538-4357/aae2ae
  * Loyd et al. (2018b) Loyd, R. O. P., France, K., Youngblood, A., et al. 2018b, ApJ, 867, 71, doi: 10.3847/1538-4357/aae2bd
  * MacGregor et al. (2013) MacGregor, M. A., Wilner, D. J., Rosenfeld, K. A., et al. 2013, ApJ, 762, L21, doi: 10.1088/2041-8205/762/2/L21
  * Mamajek et al. (2015) Mamajek, E. E., Prsa, A., Torres, G., et al. 2015, arXiv e-prints, arXiv:1510.07674.  https://arxiv.org/abs/1510.07674
  * Mariska (1992) Mariska, J. T. 1992, The Solar Transition Region (Cambridge University Press)
  * Monsignori Fossi et al. (1996) Monsignori Fossi, B. C., Landini, M., Del Zanna, G., & Bowyer, S. 1996, ApJ, 466, 427, doi: 10.1086/177522
  * Oliphant (2006) Oliphant, T. 2006, NumPy: A guide to NumPy, USA: Trelgol Publishing.  http://www.numpy.org/
  * Owen & Jackson (2012) Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931, doi: 10.1111/j.1365-2966.2012.21481.x
  * Pagano et al. (2000) Pagano, I., Linsky, J. L., Carkner, L., et al. 2000, ApJ, 532, 497, doi: 10.1086/308559
  * Paulson et al. (2006) Paulson, D. B., Allred, J. C., Anderson, R. B., et al. 2006, PASP, 118, 227, doi: 10.1086/499497
  * Peacock et al. (2019a) Peacock, S., Barman, T., Shkolnik, E. L., Hauschildt, P. H., & Baron, E. 2019a, ApJ, 871, 235, doi: 10.3847/1538-4357/aaf891
  * Peacock et al. (2019b) Peacock, S., Barman, T., Shkolnik, E. L., et al. 2019b, ApJ, 886, 77, doi: 10.3847/1538-4357/ab4f6f
  * Peacock et al. (2020) —. 2020, ApJ, 895, 5, doi: 10.3847/1538-4357/ab893a
  * Penz & Micela (2008) Penz, T., & Micela, G. 2008, A&A, 479, 579, doi: 10.1051/0004-6361:20078873
  * Plavchan et al. (2020) Plavchan, P., Barclay, T., Gagné, J., et al. 2020, Nature, 582, 497, doi: 10.1038/s41586-020-2400-z
  * Pottasch (1963) Pottasch, S. R. 1963, ApJ, 137, 945, doi: 10.1086/147569
  * Redfield et al. (2003) Redfield, S., Ayres, T. R., Linsky, J. L., et al. 2003, ApJ, 585, 993, doi: 10.1086/346129
  * Redfield et al. (2002) Redfield, S., Linsky, J. L., Ake, T. B., et al. 2002, ApJ, 581, 626, doi: 10.1086/344153
  * Ribas et al. (2018) Ribas, I., Tuomi, M., Reiners, A., et al. 2018, Nature, 563, 365, doi: 10.1038/s41586-018-0677-y
  * Sanz-Forcada et al. (2011) Sanz-Forcada, J., Micela, G., Ribas, I., et al. 2011, A&A, 532, A6, doi: 10.1051/0004-6361/201116594
  * Sanz-Forcada et al. (2010) Sanz-Forcada, J., Ribas, I., Micela, G., et al. 2010, A&A, 511, L8, doi: 10.1051/0004-6361/200913670
  * Scalo et al. (2007) Scalo, J., Kaltenegger, L., Segura, A. G., et al. 2007, Astrobiology, 7, 85, doi: 10.1089/ast.2006.0125
  * Schmelz et al. (2012) Schmelz, J. T., Reames, D. V., von Steiger, R., & Basu, S. 2012, ApJ, 755, 33, doi: 10.1088/0004-637X/755/1/33
  * Schwarz (1978) Schwarz, G. 1978, Annals of Statistics, 6, 461
  * Sekiya et al. (1980) Sekiya, M., Nakazawa, K., & Hayashi, C. 1980, Progress of Theoretical Physics, 64, 1968, doi: 10.1143/PTP.64.1968
  * Shields et al. (2016) Shields, A. L., Ballard, S., & Johnson, J. A. 2016, Phys. Rep., 663, 1, doi: 10.1016/j.physrep.2016.10.003
  * Shkolnik & Barman (2014) Shkolnik, E. L., & Barman, T. S. 2014, AJ, 148, 64, doi: 10.1088/0004-6256/148/4/64
  * Tian & Ida (2015) Tian, F., & Ida, S. 2015, Nature Geoscience, 8, 177, doi: 10.1038/ngeo2372
  * Tilipman et al. (2020) Tilipman, D., Vieytes, M., Linsky, J. L., Buccino, A. P., & France, K. 2020, arXiv e-prints, arXiv:2012.11738.  https://arxiv.org/abs/2012.11738
  * Van Grootel et al. (2018) Van Grootel, V., Fernandes, C. S., Gillon, M., et al. 2018, ApJ, 853, 30, doi: 10.3847/1538-4357/aaa023
  * Vernazza & Reeves (1978) Vernazza, J. E., & Reeves, E. M. 1978, ApJS, 37, 485, doi: 10.1086/190539
  * Warren et al. (1998) Warren, H. P., Mariska, J. T., & Lean, J. 1998, J. Geophys. Res., 103, 12077, doi: 10.1029/98JA00810
  * West et al. (2008) West, A. A., Hawley, S. L., Bochanski, J. J., et al. 2008, AJ, 135, 785, doi: 10.1088/0004-6256/135/3/785
  * West et al. (2004) West, A. A., Hawley, S. L., Walkowicz, L. M., et al. 2004, AJ, 128, 426, doi: 10.1086/421364
  * Winters et al. (2015) Winters, J. G., Henry, T. J., Lurie, J. C., et al. 2015, AJ, 149, 5, doi: 10.1088/0004-6256/149/1/5
  * Wood et al. (2012) Wood, B. E., Laming, J. M., & Karovska, M. 2012, ApJ, 753, 76, doi: 10.1088/0004-637X/753/1/76
  * Wood et al. (2018) Wood, B. E., Laming, J. M., Warren, H. P., & Poppenhaeger, K. 2018, ApJ, 862, 66, doi: 10.3847/1538-4357/aaccf6
  * Wood et al. (2005) Wood, B. E., Redfield, S., Linsky, J. L., Müller, H.-R., & Zank, G. P. 2005, ApJS, 159, 118, doi: 10.1086/430523
  * Woods et al. (2009) Woods, T. N., Chamberlin, P. C., Harder, J. W., et al. 2009, Geophys. Res. Lett., 36, L01101, doi: 10.1029/2008GL036373
  * Youngblood et al. (2016) Youngblood, A., France, K., Loyd, R. O. P., et al. 2016, ApJ, 824, 101, doi: 10.3847/0004-637X/824/2/101
  * Youngblood et al. (2017) —. 2017, ApJ, 843, 31, doi: 10.3847/1538-4357/aa76dd
  * Zahnle & Catling (2017) Zahnle, K. J., & Catling, D. C. 2017, ApJ, 843, 122, doi: 10.3847/1538-4357/aa7846