[a,b]Matteo Bresciani

# Progresses on high-temperature QCD: Equation of State and energy-momentum
tensor

Mattia Dalla Brida  Leonardo Giusti  Michele Pepe

###### Abstract

We present first non-perturbative results for the renormalization constants of
the QCD energy-momentum tensor, based on the framework of thermal QCD with
shifted and twisted (for quarks only) boundary conditions in the compact
direction. We also show preliminary results for the entropy density obtained
with the very same numerical strategy. This opens the way to the determination
of the QCD Equation of State up to very high temperatures.

CERN-TH-2023-219

##  1 Introduction

Quantum Chromodynamics (QCD) is the fundamental field theory describing the
interaction among quarks and gluons. Its thermal state from the deconfinement
temperature $150$-$200$ MeV up to the electroweak scale $\sim 100$ GeV, the
so-called Quark-Gluon Plasma (QGP), is of great interest for instance in
cosmology and high energy physics. Indeed the Big Bang model describes the
Universe in some of its early stages as a plasma of quarks and gluons.
Moreover, the QGP is currently under investigation in heavy-ion collision
experiments. On the theoretical side, many physical properties of the QGP are
encoded by the correlation functions of the energy-momentum tensor (EMT) of
QCD, see Ref. [1] for a recent review. They can be computed from first
principles on the lattice once the discretized EMT has been properly
renormalized. We propose a numerical strategy to determine the renormalization
constants of the EMT non-perturbatively, based on considering QCD with shifted
and twisted boundary conditions [7, 8]. The very same framework allows for the
determination of the Equation of State (EoS) of QCD up to very high
temperatures. At present the EoS is known non-perturbatively for temperatures
below $1$ GeV [2, 3] because of some technical limitations of the state of the
art techniques, which are overcome by our alternative strategy. An analogous
approach has been exploited for the SU$(3)$ Yang-Mills theory in Ref. [6],
where perturbation theory was shown to be unreliable up to very high
temperatures.

##  2 The energy-momentum tensor in the continuum

We consider the symmetric and gauge-invariant definition of the QCD EMT [4,
5]. In the continuum the latter splits into a singlet and a 9-dimensional
representation of the SO$(4)$ group. Under renormalization, the singlet
component mixes with the identity operator, while there are no other fields of
mass dimension $\leq 4$ which may mix with the non-singlet representation.
Moreover, the bare and renormalized non-singlet EMT coincide because the
latter is a conserved current associated to the invariance of the theory under
the group of space-time translations. This property becomes manifest if we
formulate QCD at finite temperature $T$, in a moving reference frame at
Euclidean speed $\boldsymbol{\xi}$ and with an imaginary chemical potential
$\mu_{I}$ [7, 8]. This setup is realized with the following boundary
conditions in the (compact) temporal direction for the gauge and quark fields
[8]:

$A_{\mu}(x_{0}+L_{0},\boldsymbol{x})=A_{\mu}(x_{0},\boldsymbol{x}-L_{0}\boldsymbol{\xi})\;,$
$\psi(x_{0}+L_{0},\boldsymbol{x})=-e^{i\theta_{0}}\,\psi(x_{0},\boldsymbol{x}-L_{0}\boldsymbol{\xi})\;,$
(2.1)
${\bar{\psi}}(x_{0}+L_{0},\boldsymbol{x})=-e^{-i\theta_{0}}\,{\bar{\psi}}(x_{0},\boldsymbol{x}-L_{0}\boldsymbol{\xi})\;.$

The Euclidean boost is encoded by the $L_{0}\boldsymbol{\xi}$ shift ($L_{0}$
is the temporal extension) in the spatial directions at the temporal boundary.
The usual relation between the temperature and $L_{0}$ is modified by the
Euclidean Lorentz factor in $T=L_{0}^{-1}(1+\boldsymbol{\xi}^{2})^{-1/2}$, and
the twist phase $\theta_{0}=-L_{0}\mu_{I}$ for the fermions is related to the
imaginary chemical potential. The space-time components of the EMT satisfy the
following identity [8],

$\langle
T_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}=-\dfrac{\partial
f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}}{\partial\xi_{k}}\;,\quad
f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}=-\frac{1}{L_{0}L^{3}}\ln{\cal
Z}(L_{0},L,\boldsymbol{\xi},\theta_{0})\;,$ (2.2)

where
$f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}$
is the free-energy density of thermal QCD and $\mathcal{Z}$ its partition
function. The free-energy density is a spectral quantity with an additive
power-divergent term. Its derivative with respect to the shift
$\boldsymbol{\xi}$ is thus finite, once the bare parameters of the theory have
been renormalized. Therefore $T_{0k}$ and all the other non-singlet components
of the EMT are finite too. The traceless diagonal components of the EMT are
related to the space-time components by means of Lorentz transformations [7,
8]:

$\langle
T_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}=\xi_{k}\left(\langle
T_{00}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}-\langle
T_{jj}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}\right)\quad(j\neq
k,\xi_{j}=0)\;.$ (2.3)

Finally, the dependence of $\langle
T_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}$
on the twist phase $\theta_{0}$ can be written as [8]

$\langle
T_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{A}}}-\langle
T_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{B}}}=\dfrac{i}{L_{0}}\int_{\theta_{0}^{A}}^{\theta_{0}^{B}}d\theta_{0}\dfrac{\partial}{\partial\xi_{k}}\langle
V_{0}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}\equiv-{{\cal
V}_{0,k}^{AB}}\;,$ (2.4)

where $V_{\mu}$ is the flavour-singlet vector current. We will take advantage
of these relations in our lattice calculations.

##  3 Renormalization of the energy-momentum tensor on the lattice

The lattice regularization breaks the SO$(4)$ symmetry down to its discrete
hypercubic subgroup SW4. The non-singlet EMT components split into the sextet
and triplet representations of this discrete group [5]. These two
representations renormalize separately, and no other lattice field is involved
in their renormalization. Since only discrete space-time translations are
allowed, the bare non-singlet EMT is not a conserved current. As a
consequence, the renormalization constants of the operators in each
representation are finite functions of the bare coupling $g_{0}$ and
independent from the external states, up to discretization effects. We focus
on the non-perturbative definition of the sextet and triplet representations
of the lattice EMT, which follow the renormalization pattern

$T_{\mu\nu}^{R,\\{i\\}}=Z_{G}^{\\{i\\}}T_{\mu\nu}^{G,\\{i\\}}+Z_{F}^{\\{i\\}}T_{\mu\nu}^{F,\\{i\\}},\;\;\;\;i=3,6\;.$
(3.1)

$T_{\mu\nu}^{G,\\{i\\}}$ and $T_{\mu\nu}^{F,\\{i\\}}$ are the dimension-4
gauge-invariant lattice operators in the $i-$dimensional representation of the
hypercubic group [5, 8].

The lattice shifted and twisted boundary conditions are analogous to the
continuum ones (2.1), but with the gauge potential $A_{\mu}(x)$ replaced by
the link field $U_{\mu}(x)$. In addition, we consider the usual periodic
boundary conditions for the three compact spatial directions of the lattice
box, each of size $L$.

Using the definition (3.1) together with equation (2.2) at two values of the
twist phase, $\theta_{0}^{A}$ and $\theta_{0}^{B}$, and exploiting equation
(2.4), we can write our master formula for the computation of the
renormalization constants of the sextet representation:

$\displaystyle\begin{pmatrix}\langle
T^{G,\\{6\\}}_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{A}}}&&\langle
T^{F,\\{6\\}}_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{A}}}\\\
\\\ \\\ \langle
T^{G,\\{6\\}}_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{B}}}&&\langle
T^{F,\\{6\\}}_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{B}}}\end{pmatrix}\begin{pmatrix}Z_{G}^{\\{6\\}}\\\
\\\ \\\ Z_{F}^{\\{6\\}}\end{pmatrix}=\begin{pmatrix}-\dfrac{\Delta
f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{A}}}}{\Delta\xi_{k}}\\\
\\\ -\dfrac{\Delta
f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{A}}}}{\Delta\xi_{k}}+{{\cal
V}_{0,k}^{AB}}\end{pmatrix}+{\rm O}(a^{2})\;.$ (3.2)

The partial derivative has been replaced by finite differences up to
discretization effects. Then, the lattice counterpart of the Ward Identity
(2.3)

$\langle
T^{R,\\{6\\}}_{0k}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}=\xi_{k}\langle
T^{R,\\{3\\}}_{0j}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}\;\;\;\;(j\neq
k,\xi_{j}=0)\;$ (3.3)

gives access to the renormalization constants of the triplet representation.
By dimensional arguments, the leading lattice artifacts of the renormalization
constants when $a\to 0$ are

$Z\left(g_{0}^{2},aT\right)=Z(g_{0}^{2})+C_{1}\cdot\left(aT\right)^{2}+C_{2}\cdot(a\Lambda_{\text{QCD}})aT+C_{3}\cdot(a\Lambda_{\text{QCD}})^{2}+\;...\;,$
(3.4)

where $Z$ is a shorthand for any of the renormalization constants in equation
(3.1). We include a $T\to 0$ limit at fixed $g_{0}^{2}$ in the definition of
the renormalization constants, so that the leading discretization effects
$C_{1}$, $C_{2}$ are suppressed. All the residual lattice artifacts will
disappear when a renormalized correlator of the lattice EMT will be
extrapolated to the continuum limit.

##  4 Renormalization strategy at work

This procedure proved to be effective in the non-perturbative renormalization
of the energy-momentum tensor in the SU$(3)$ Yang-Mills theory on the lattice
[9], and of the flavour-singlet local vector current [10] in lattice QCD. In
the following we give some details on the lattice determination of the
quantities appearing in the master equation (3.2). Although the strategy is
general, we choose the Wilson formulation of lattice QCD with ${N_{f}}=3$
flavours of massless O($a$)-improved Wilson fermions.

###  4.1 Shift derivative of the free-energy density

In order to estimate the shift derivative of the free-energy from a lattice
calculation we add and subtract inside the discrete differential operator, at
fixed $g_{0}^{2}$ and $L_{0}/a$, the free-energy density
$f^{\infty}_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}$
of QCD with quarks at infinite mass:

$\dfrac{\Delta
f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}}{\Delta\xi_{k}}=\dfrac{\Delta}{\Delta\xi_{k}}\left(f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}^{\infty}-\int_{0}^{\infty}dm_{q}\;\dfrac{\partial
f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}}{\partial
m_{q}}\right)=\dfrac{\Delta}{\Delta\xi_{k}}\left(f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}^{\infty}-\int_{0}^{\infty}dm_{q}\;\langle{\bar{\psi}}\psi\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}\right)\;.$
(4.1)

We also replace
$f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}-f^{\infty}_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}$
with an integral in the bare subtracted quark mass $m_{q}=m_{0}-m_{c}$
($m_{c}$ is the critical mass) of the derivative in $m_{q}$ of
$f_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}$.
This derivative is the chiral condensate, as we write in the last step. The
integral in the mass is a well defined quantity because the additive
divergences of the chiral condensate are canceled by the overall
differentiation in the shift, while its multiplicative renormalization drops
together with the renormalization of the quark mass in the integration
measure.

Quarks at infinite mass decouple from QCD, therefore the shift derivative of
$f^{\infty}_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}$
can be computed in the ${N_{f}}=0$ theory as described in [9]. We estimate the
integral in the quark mass in lattice QCD, and the left panel of Figure 1
shows a preliminary plot of the integrand function coming from lattice
simulations on a $6\times 96^{3}$ volume and at $\beta=8.8727$. The integral
is estimated with two Gauss quadratures, one for the peak (5 points) and one
for the tail (10 points). The integration interval is cut when the remainder
of the integral is expected to be negligible compared with the target
accuracy. Using this strategy we obtain the area below the curve with a
relative precision of a few permille.

Figure 1: Left: integrand function of eq. (4.1) sampled at the mass values
prescribed by the Gauss quadratures for peak and tail. Right: integrand
function of eq. (4.2) sampled at the values of the twist phase $\theta_{0}$
prescribed by the Gauss quadrature.

###  4.2 Integral in the twist phase

We move to the determination of the integral ${{\cal V}_{0,k}^{AB}}$ in
equation (3.2), defined in equation (2.4). On the lattice we have

${{\cal
V}_{0,k}^{AB}}=-\dfrac{i}{L_{0}}\int_{\theta_{0}^{A}}^{\theta_{0}^{B}}d\theta_{0}\;\dfrac{\Delta}{\Delta\xi_{k}}\langle
V_{0}^{R}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}=-\dfrac{i}{L_{0}}\int_{(\theta_{0}^{A})^{2}}^{(\theta_{0}^{B})^{2}}d\theta_{0}^{2}\;\frac{1}{2\theta_{0}}\dfrac{\Delta}{\Delta\xi_{k}}\langle
V_{0}^{R}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}}}\;,$
(4.2)

and in place of $V_{0}^{R}$ we may use the conserved current $V_{0}^{c}$, or
the renormalized local current $Z_{V}V_{0}^{l}$. For this study, we choose the
former. As a last technical step, we perform a change of variable from
$\theta_{0}$ to $\theta_{0}^{2}$, and we estimate the resulting integral with
a $n=3$ Gauss quadrature. The change of variable is convenient because the
integrand function is odd in $\theta_{0}$, meaning that a polynomial
representation of it, like the Gauss quadrature, would contain odd powers of
$\theta_{0}$ only. With the change to $\theta_{0}^{2}$ we effectively halve
the degree of the polynomial, so that with a 3-point quadrature we obtain the
same accuracy that we would get with a 6-point quadrature of the original
integral. The interplay of the twist phase with the $\mathbb{Z}_{3}$ center
symmetry of the gauge group SU$(3)$ makes the free-energy density effectively
periodic in $\theta_{0}$, with period $2\pi/3$ [11]. The choice
$\theta_{0}^{A}=0$ and $\theta_{0}^{B}=3\pi/10$ covers the most of the allowed
interval and thus maximizes the signal. The right panel of Figure 1 shows the
integrand function, compared with its 1-loop result at finite $L_{0}/a$. The
lattice setup is the same as for the integral in the mass, and as in that case
we obtain a relative error of a few permille on the integral.

###  4.3 1-point functions of the bare EMT

The last missing ingredients in equation (3.2) are the 1-point functions of
the bare EMT appearing in the matrix on the left-hand side. Table 1 shows some
preliminary numerical results at the usual bare parameters. Data at
$\theta_{0}=0$ come from ensembles with $L/a=288$ while data at
$\theta_{0}=3\pi/10$ from ensembles with $L/a=96$. Finite volume effects are
expected to be negligible because exponentially suppressed with the mass of
the lightest screening state [7], which is proportional to the temperature.
The 1-point functions are averaged over the volume, and the number of
measurements is tuned so as to reach a comparable statistical precision in the
two ensembles. The relative error on the $\langle T^{F}\rangle$ components is
a few permille, while on the $\langle T^{G}\rangle$ components is about $1\%$.

$\theta_{0}^{A}=0$
---
$\langle T^{F,\\{6\\}}_{01}\rangle/T^{4}$ | -6.343(9) | 0.15$\%$
$\langle T^{G,\\{6\\}}_{01}\rangle/T^{4}$ | -2.822(24) | 0.85$\%$
$\langle T^{F,\\{3\\}}_{02}\rangle/T^{4}$ | -6.846(11) | 0.17$\%$
$\langle T^{G,\\{3\\}}_{02}\rangle/T^{4}$ | -3.13(4) | 1.20$\%$
$L/a=288$, $n_{\rm trj}=100$

$\theta_{0}^{B}=3\pi/10$
---
$\langle T^{F,\\{6\\}}_{01}\rangle/T^{4}$ | -4.053(9) | 0.21$\%$
$\langle T^{G,\\{6\\}}_{01}\rangle/T^{4}$ | -2.677(27) | 1.00$\%$
$\langle T^{F,\\{3\\}}_{02}\rangle/T^{4}$ | -4.375(15) | 0.34$\%$
$\langle T^{G,\\{3\\}}_{02}\rangle/T^{4}$ | -2.86(4) | 1.53$\%$
$L/a=96$, $n_{\rm trj}=2000$

Table 1: Preliminary results for the 1-point functions of the sextet and
triplet bare EMT. The bare parameters are $L_{0}/a=6$, $\beta=8.8727$.

###  4.4 Renormalization constants

Using the results described so far in equation (3.2) we obtain a first
estimation of the renormalization constants of the sextet energy-momentum
tensor. Then, equation (3.3) leads to the renormalization constants of the
triplet representation. Figure 2 shows the relative errors on the four
renormalization constants, and the pie charts give the breakdown of the
contributions to the variances. The error is largely dominated by the one of
the 1-point functions $\langle T^{G}\rangle$. The reason can be traced back to
the fact that the solution of the master equation (3.2) depends on differences
like $\langle
T^{G}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{B}}}-\langle
T^{G}\rangle_{{\scriptscriptstyle\boldsymbol{\xi},}{\scriptscriptstyle\theta_{\scalebox{0.5}{0}}^{A}}}$.
These 1-point functions are weakly dependent on $\theta_{0}$, therefore the
signal mostly cancels while the errors, which are basically independent from
the value of $\theta_{0}$, sum in quadrature.

Figure 2: Breakdown of the contributions to the variance of the
renormalization constants.

##  5 Entropy density

The numerical determination of the shift derivative of the free-energy density
is closely related to the computation of the entropy density of QCD [7, 8]:

$\dfrac{s}{T^{3}}=\dfrac{1+\boldsymbol{\xi}^{2}}{\xi_{k}}\dfrac{1}{T^{4}}\dfrac{\Delta
f_{\boldsymbol{\xi}}}{\Delta\xi_{k}}=\dfrac{1+\boldsymbol{\xi}^{2}}{\xi_{k}}\dfrac{1}{T^{4}}\left(\dfrac{\Delta
f^{\infty}_{\boldsymbol{\xi}}}{\Delta\xi_{k}}-\int_{0}^{\infty}dm_{q}\dfrac{\Delta\langle{\bar{\psi}}\psi\rangle_{\boldsymbol{\xi}}}{\Delta\xi_{k}}\right)\;.$
(5.1)

The free-energy is expressed in terms of the infinite quark mass case and of
an integral in the quark mass of the shift derivative of the chiral
condensate, similarly to what described in Subsection 4.1. Using these
formulas, and the lattice data at $L_{0}/a=6$, $\beta=8.8727$ (corresponding
to a temperature of about $140$ GeV [12]) we estimate the entropy density with
a relative error of half a percent. After the tree-level improvement of the
lattice data (which is expected to remove the leading discretization effects
at small $g_{0}^{2}$) we observe a deviation of a few percent with respect to
the Stefan-Boltzmann limit.

##  6 Conclusions

In these proceedings we described our ongoing work on the non-perturbative
renormalization of the QCD energy-momentum tensor, and on the first principle
determination of the Equation of State of the Quark-Gluon Plasma at high
temperature. We carry out our computations in the framework of QCD in a moving
reference frame and with an imaginary chemical potential, where some
convenient Ward Identities can be employed to constrain the renormalization
constants of the energy-momentum tensor, and to compute the entropy density.
We showed some preliminary non-perturbative results from a lattice setup of
temporal extension $L_{0}/a=6$ and $\beta=8.8727$, where we obtained a final
accuracy of $\sim 3\%$ on the renormalization constants $Z_{F}^{\\{3,6\\}}$,
and of $\sim 5\%$ on $Z_{G}^{\\{3,6\\}}$. On the same ensembles, the entropy
density is determined with a relative error of $\sim 0.5\%$. These results are
encouraging for the extension of this computation to the temperature range $1$
GeV - $100$ GeV considered in Ref. [12].

##  7 Acknowledgements

We acknowledge PRACE for awarding us access to the HPC system MareNostrum4 at
the Barcelona Supercomputing Center (Proposals n. 2018194651 and 2021240051)
where most of the numerical results presented in this paper have been
obtained. We also thank CINECA for providing us with computer time on Marconi
(CINECA-INFN, CINECA-Bicocca agreements). The R&D has been carried out on the
PC clusters Wilson and Knuth at Milano-Bicocca. We thank all these
institutions for the technical support. This work is (partially) supported by
ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and
Quantum Computing, funded by European Union – NextGenerationEU. MB would like
to thank CERN for the kind hospitality and support extended to him during the
write-up of this contribution.

## References

  * [1] H. B. Meyer, Eur. Phys. J. A 47 (2011), 86 doi:10.1140/epja/i2011-11086-3 [arXiv:1104.3708 [hep-lat]].
  * [2] A. Bazavov et al. [HotQCD], Phys. Rev. D 90 (2014), 094503 doi:10.1103/PhysRevD.90.094503 [arXiv:1407.6387 [hep-lat]].
  * [3] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabo, Phys. Lett. B 730 (2014), 99-104 doi:10.1016/j.physletb.2014.01.007 [arXiv:1309.5258 [hep-lat]].
  * [4] L. D. Landau and E. M. Lifschits, Pergamon Press, 1975
  * [5] S. Caracciolo, G. Curci, P. Menotti and A. Pelissetto, Annals Phys. 197 (1990), 119 doi:10.1016/0003-4916(90)90203-Z
  * [6] L. Giusti and M. Pepe, Phys. Lett. B 769 (2017), 385-390 doi:10.1016/j.physletb.2017.04.001 [arXiv:1612.00265 [hep-lat]].
  * [7] L. Giusti and H. B. Meyer, JHEP 01 (2013), 140 doi:10.1007/JHEP01(2013)140 [arXiv:1211.6669 [hep-lat]].
  * [8] M. Dalla Brida, L. Giusti and M. Pepe, JHEP 04 (2020), 043 doi:10.1007/JHEP04(2020)043 [arXiv:2002.06897 [hep-lat]].
  * [9] L. Giusti and M. Pepe, Phys. Rev. D 91 (2015), 114504 doi:10.1103/PhysRevD.91.114504 [arXiv:1503.07042 [hep-lat]].
  * [10] M. Bresciani, M. Dalla Brida, L. Giusti, M. Pepe and F. Rapuano, Phys. Lett. B 835 (2022), 137579 doi:10.1016/j.physletb.2022.137579 [arXiv:2203.14754 [hep-lat]]
  * [11] A. Roberge and N. Weiss, Nucl. Phys. B 275 (1986), 734-745 doi:10.1016/0550-3213(86)90582-1
  * [12] M. Dalla Brida, L. Giusti, T. Harris, D. Laudicina and M. Pepe, JHEP 04 (2022), 034 doi:10.1007/JHEP04(2022)034 [arXiv:2112.05427 [hep-lat]].


# The JWST Resolved Stellar Populations Early Release Science Program VI.
Identifying Evolved Stars in Nearby Galaxies

Martha L. Boyer Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Giada Pastorelli Padova Astronomical Observatory,
Vicolo dell’Osservatorio 5, Padova, Italy  Léo Girardi Padova Astronomical
Observatory, Vicolo dell’Osservatorio 5, Padova, Italy  Paola Marigo
Department of Physics and Astronomy G. Galilei, University of Padova, Vicolo
dell’Osservatorio 3, I-35122, Padova, Italy  Andrew E. Dolphin RTX, 1151 E.
Hermans Road, Tucson, AZ 85756, USA  Steward Observatory, University of
Arizona, 933 N. Cherry Avenue, Tucson, AZ 85719, USA  Kristen. B. W. McQuinn
Rutgers University, Department of Physics and Astronomy, 136 Frelinghuysen
Road, Piscataway, NJ 08854, USA  Max J. B. Newman Rutgers University,
Department of Physics and Astronomy, 136 Frelinghuysen Road, Piscataway, NJ
08854, USA  Alessandro Savino Department of Astronomy, University of
California, Berkeley, CA 94720, USA  Daniel R. Weisz Department of Astronomy,
University of California, Berkeley, CA 94720, USA  Benjamin F. Williams
Department of Astronomy, University of Washington, Box 351580, U.W., Seattle,
WA 98195-1580, USA  Jay Anderson Space Telescope Science Institute, 3700 San
Martin Drive, Baltimore, MD 21218, USA  Roger E. Cohen Rutgers University,
Department of Physics and Astronomy, 136 Frelinghuysen Road, Piscataway, NJ
08854, USA  Matteo Correnti INAF Osservatorio Astronomico di Roma, Via
Frascati 33, 00078, Monteporzio Catone, Rome, Italy  ASI-Space Science Data
Center, Via del Politecnico, I-00133, Rome, Italy  Andrew A. Cole School of
Natural Sciences, University of Tasmania, Private Bag 37, Hobart, Tasmania
7001, Australia  Marla C. Geha Department of Astronomy, Yale University, New
Haven, CT 06520, USA  Mario Gennaro Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD 21218, USA  The William H. Miller III
Department of Physics & Astronomy, Bloomberg Center for Physics and Astronomy,
Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA
Nitya Kallivayalil Department of Astronomy, The University of Virginia, 530
McCormick Road, Charlottesville, VA 22904, USA  Evan N. Kirby Department of
Physics and Astronomy, University of Notre Dame, 225 Nieuwland Science Hall,
Notre Dame, IN 46556, USA  Karin M. Sandstrom Department of Astronomy &
Astrophysics, University of California San Diego, 9500 Gilman Drive, La Jolla,
CA 92093, USA  Evan D. Skillman University of Minnesota, Minnesota Institute
for Astrophysics, School of Physics and Astronomy, 116 Church Street, S.E.,
Minneapolis, MN 55455, USA  Christopher T. Garling Department of Astronomy,
The University of Virginia, 530 McCormick Road, Charlottesville, VA 22904, USA
Hannah Richstein Department of Astronomy, The University of Virginia, 530
McCormick Road, Charlottesville, VA 22904, USA  Jack T. Warfield Department of
Astronomy, The University of Virginia, 530 McCormick Road, Charlottesville, VA
22904, USA

###### Abstract

We present an investigation of evolved stars in the nearby star-forming galaxy
WLM, using NIRCam imaging from the JWST resolved stellar populations early-
release science (ERS) program. We find that various combinations of the F090W,
F150W, F250M, and F430M filters can effectively isolate red supergiants (RSGs)
and thermally-pulsing asymptotic giant branch (TP-AGB) stars from one another,
while also providing a reasonable separation of the primary TP-AGB subtypes:
carbon-rich C-type stars and oxygen-rich M-type stars. The classification
scheme we present here agrees very well with the well-established Hubble Space
Telescope (HST) medium-band filter technique. The ratio of C to M-type stars
(C/M) is 0.8$\pm$0.1 for both the new JWST and the HST classifications, which
is within one sigma of empirical predictions from optical narrow-band CN and
TiO filters. The evolved star colors show good agreement with the predictions
from the PARSEC$+$COLIBRI stellar evolutionary models, and the models indicate
a strong metallicity dependence that makes stellar identification even more
effective at higher metallicity. However, the models also indicate that
evolved star identification with NIRCam may be more difficult at lower
metallicies. We test every combination of NIRCam filters using the models and
present additional filters that are also useful for evolved star studies. We
also find that $\approx$90% of the dusty evolved stars are carbon-rich,
suggesting that carbonaceous dust dominates the present-day dust production in
WLM, similar to the findings in the Magellanic Clouds. These results
demonstrate the usefulness of NIRCam in identifying and classifying dust-
producing stars without the need for mid-infrared data.

JWST (2291), Hubble Space Telescope (761), Asymptotic giant branch stars
(2100), Red supergiant stars (1375), Local Group (929)

††facilities: HST(WFC3), JWST(NIRCam)††software: Astropy (Astropy
Collaboration et al., 2013, 2018, 2022), DOLPHOT (Dolphin, 2000, 2016), Scipy
(Virtanen et al., 2020), TRILEGAL (Girardi et al., 2005)

##  1 Introduction

The thermally-pulsing Asymptotic Giant Branch (TP-AGB) stage is the end-phase
of stellar evolution for $\sim$1–8 $M_{\odot}$ stars, marked by complex
physical processes such as pulsation, mass loss, hot bottom burning, and
dredge up (Habing & Olofsson, 2004; Karakas, 2017). TP-AGB stars are
relatively rare due to the short lifetimes ($\sim 0.5$–$20$ Myr) of this
phase, yet they have an outsized impact on the evolution of their host
galaxies owing to their production of metals and dust. Their luminosities are
also extremely high, affecting interpretations of galaxy observations.

Stellar models have historically struggled to reproduce TP-AGB observables
owing to both the complexity of this stage of evolution and to a scarcity of
high-quality calibration data. Recent progress has been made in calibrating
the TP-AGB phase (e.g., Rosenfield et al., 2014, 2016) using data of nearby
galaxies, but many of these efforts were unable to separate different TP-AGB
subtypes in the data, which limits the scope of the calibration. The only
galaxies with high quality, complete TP-AGB calibration data (including
stellar subtypes) are the Magellanic Clouds, for which there is an abundance
of near- to mid-infrared (IR) data. As a result, most models are anchored to
the Magellanic Clouds (e.g., Pastorelli et al., 2019, 2020), leaving the
metallicity dependencies poorly constrained. These model uncertainties
propagate when TP-AGB stars are present in observations that are used to
derive various galaxy properties.

The presence of TP-AGB stars can make precise measurements of the tip of the
red giant branch (TRGB) difficult since TP-AGB stars can overlap with RGB
stars. In addition, fitting a stellar population to derive model-dependent
parameters such as age, star-formation history (SFH), or metallicity can
result in highly inaccurate estimates if the TP-AGB stars are not excluded or
fully accounted for (Conroy et al., 2009; Baldwin et al., 2018). Moreover,
failure to separate oxygen-rich M-type AGB stars from carbon-rich C-type AGB
stars can also potentially affect distance estimates measured from period-
luminosity (P-L) diagrams and the ‘JAGB’ method (Freedman & Madore, 2020; Lee,
2023).

Classifying TP-AGB stars can also provide useful galaxy diagnostics. For
example, the ratio of C- to M-type stars (C/M) can provide information about a
galaxy’s metallicity and/or age (e.g., Groenewegen, 2006a; Boyer et al.,
2019). Also, since each subtype produces a different type of dust
(carbonaceous or siliceous), separating them allows for a better understanding
of the source of a galaxy’s interstellar dust: AGB stars, supernovae, or grain
growth in the interstellar medium itself.

Yet, identifying TP-AGB stars, and especially classifying them into their main
subtypes (M- and C-type) and separating them from other evolved stars like red
supergiants (RSGs), is challenging without spectroscopy and pulsation
information. While a fraction of TP-AGB stars is always found below the TRGB
because of the luminosity dips of their thermal pulse cycles, in the optical
wavelengths $<$1 µm this fraction increases due to molecular absorption in
their atmospheres and/or to variability and extinction from circumstellar dust
(Boyer et al., 2019). This is less of an issue in the near-IR, where molecular
absorption is not as strong, pulsation amplitudes are smaller, and extinction
is weaker. However, clean and reliable separation of M- and C-type AGB stars
in the near-IR has so far proven impossible with broadband filters (e.g.,
Girardi et al., 2010; Boyer et al., 2013, 2015; Rosenfield et al., 2014).

Recent studies have shown that medium-band filters on HST’s Wide Field Camera
3 (WFC3/IR) cleanly sample key molecular features and can thus effectively and
efficiently identify TP-AGB stars and classify them into their subtypes.
However, given the limitations in HST’s resolution and sensitivity (and its
finite lifetime), finding effective filter combinations with JWST’s Near-
Infrared Camera (NIRCam; Rieke et al., 2005; Gardner et al., 2023; Rigby et
al., 2023; Rieke et al., 2023) is prudent.

To investigate JWST’s ability to photometrically classify evolved stars, we
target the nearby dwarf irregular galaxy Wolf–Lundmark–Melotte (WLM) as part
of the JWST Stellar Populations Early Release Science (ERS) program (Boyer et
al., 2022; Warfield et al., 2023; Weisz et al., 2023; McQuinn et al., 2023).
Wolf (1909) discovered the galaxy on photographic plates, and Melotte (1926,
along with K. Lundmark) re-discovered the galaxy and noted its similar
appearance to NGC 6822, the nearest dwarf irregular galaxy. WLM is seen nearly
edge-on to its rotating disk of gas and stars (Leaman et al., 2012). Like most
dwarf irregular galaxies, WLM has ongoing star formation
($0.008~{}M_{\sun}~{}{\rm yr}^{-1}$ Mondal et al., 2018). It also has a star
formation history weighted toward recent ages, with the majority of star
formation occurring in the last 5 Gyr (Weisz et al., 2014; Albers et al.,
2019; McQuinn et al., 2023). WLM’s stellar metallicity distribution is typical
for dwarf galaxies of its mass, with a mean iron abundance of
${\rm[Fe/H]}=-1.27$ and a shallow radial metallicity gradient (Leaman et al.,
2009, 2013).

One goal of our ERS program is to demonstrate whether a filter combination
that is popular for many galaxy studies can enable accurate AGB and RSG
identification and classification. In this paper, we explore NIRCam data for
WLM and compare the evolved star population to Hubble observations and to that
predicted by stellar models. In §2 we describe the observations and data
reduction. In §3 we outline our selection of TP-AGB stars and other objects.
In §4 we discuss our final selection recommendation and galaxy properties that
might affect it, such as the metallicity.

##  2 Data

The NIRCam images of WLM111https://doi.org/10.17909/71kb-ga31 (catalog DOI:
10.17909/71kb-ga31) were observed as part of ERS program 1334 (PI Weisz) in
July 2022 using the F090W and F150W filters in the short wavelength (SW)
channel and the F250M and F430M filters in the long wavelength (LW) channel.
See Weisz et al. (2023) for full details about the observations. The exposure
setup included 8 groups with the MEDIUM8 readout pattern with 7–9 integrations
and 4 subpixel dithers (4-POINT-MEDIUM-WITH-NIRISS). The total exposure time
was $\approx$24–30 ks. NIRCam’s Module A covers a central field in WLM, while
Module B covers a field to the south.

The data were processed with JWST Pipeline version 1.11.4, with CRDS context
jwst_p1147.pmap, using the default parameters in the pipeline through Stage-2
processing. Photometry was performed using the new NIRCam module for DOLPHOT
(Dolphin, 2000, 2016; Weisz et al., 2023), using parameters described by Weisz
et al. (2024). The 50% photometric completeness limit is $\gtrsim$28 mag in
the broadband filters and $\sim$25 mag in the medium band filters, both of
which are several magnitudes fainter than the TRGB. The photometric
uncertainties are within 0.01 mag. The photometry was culled following
Warfield et al. (2023), see §3.3.4. A handful of stars in our images are
saturated, so it is possible that a few of the most extreme examples of TP-AGB
stars are not included in our catalog, but we expect the TP-AGB population is
very near fully complete in the NIRCam FOV.

Our analysis also includes Hubble Space Telescope (HST) data from the Local
UltraViolet and Infrared Treasury (LUVIT) program (K. Gilbert et al., in prep,
and M. Boyer et al., in prep). That program includes imaging with the Advanced
Camera for Surveys (ACS), Wide Field Camera 3 (WFC3), and Wide-Field and
Planetary Camera 2 (WFPC2) instruments in at least nine filters from F275W to
F160W, obtained as part of programs GO-15275 (PI: K. Gilbert) and GO-16162
(PI: M. Boyer). Photometry on those data was performed using DOLPHOT and the
photometry pipeline initially developed for the Panchromatic Hubble Andromeda
Treasury (PHAT) program (Dalcanton et al., 2012; Williams et al., 2014). Here,
we focus mainly on the IR medium-band filters from GO-16162: F127M, F139M, and
F153M, which are particularly useful for separating M-type from C-type AGB
stars (§3.2). Figure 1 shows the overlap between the JWST and HST fields of
view.

Figure 1: Map of the SW NIRCam point sources (black dots), with the medium-
band HST WFC3/IR footprints from GO-16162 overlaid in black. North is to the
left and east is down. Each NIRCam module (4 SW detectors with
$\approx$4″small gaps in between) is 2$\farcm$2$\times$2$\farcm$2\.

Throughout our analysis, we assume a foreground extinction of $A_{V}=0.10$
based on Galactic dust maps (Schlegel et al., 1998; Schlafly & Finkbeiner,
2011), and adopt the $A_{\lambda}$ values using the Cardelli et al. (1989) $+$
O’Donnell (1994) extinction curve with $R_{V}=3.1$ and the post-launch JWST
filters transmissions updated in Nov 2022.222Derived using the CMD v3.2
website

##  3 Analysis

###  3.1 TRGB Estimates

At most wavelengths, TP-AGB stars are usually brighter than the TRGB, making
the TRGB a key constraint for identifying AGB stars. However, multiple factors
can push TP-AGB stars fainter than this limit. Thermal pulses, variability
and/or strong molecular absorption in the atmosphere can cause these stars to
periodically dip below the TRGB at some wavelengths, especially in optical and
near-IR bands where these effects cause the largest brightness changes.
Circumstellar dust can also cause a highly evolved TP-AGB star to appear
fainter than the TRGB in bluer bands. For all of these reasons, including the
TRGB at longer wavelengths in the selection criteria results in a more
complete TP-AGB sample.

To measure the TRGB in each band, we follow the procedure in Méndez et al.
(2002). We create a Gaussian-smoothed luminosity function for each filter and
apply an edge-detection Sobel kernel of the form $[-2,0,+2]$. We repeat this
with 1000 Monte Carlo bootstrap resampling trials, each time adding random
1-$\sigma$ noise to the photometry and randomly varying the starting magnitude
of the luminosity function bins. The final TRGB magnitude and its uncertainty
are the peak and standard deviation of a Gaussian function fit to the
distribution of TRGBs measured in the trials; these are listed in Table 1.
With our adopted distance modulus of 24.93 mag (Albers et al., 2019), the
absolute magnitude of the F090W TRGB is $M_{\rm TRGB}^{\rm F090W}=-4.31$ mag,
within 0.05 mag of the F090W estimate recently reported by Anand et al.
(2024).

Table 1: Measured TRGBs  Filter | TRGB |  $\sigma$TRGB
---|---|---
| (mag) | (mag)
—JWST NIRCam—
F090W | 20.62 | 0.03
F150W | 19.27 | 0.05
F250M | 18.98 | 0.07
F430M | 18.80 | 0.05
—HST WFC3/IR—
F127M | 19.79 | 0.04
F139M | 19.57 | 0.04
F153M | 19.26 | 0.04

###  3.2 Selecting TP-AGB stars with HST filters

The medium-band HST filters overlap with the H2O and CN+C2 features in M-type
and C-type TP-AGB spectra, respectively. A combination of the F127M, F139M,
and F153M filters are therefore well-suited to isolating these types of TP-AGB
stars (Boyer et al., 2013, 2017, 2019; Goldman et al., 2019; Jones et al.,
2023). The LUVIT HST fields overlap with the ERS JWST fields (Fig. 1), so we
can use the HST classifications to guide TP-AGB classification in the ERS
filters.

Figure 2: HST medium-band color-color diagram showing the separation between
C-type and M-type TP-AGB stars in the overlapping HST$+$JWST fields (Fig. 1).
All points in this diagram are brighter than the TRGB in at least one of the
HST filters. Grey points are dominated by foreground, main sequence, and RSGs
(Boyer et al., 2013, 2017). The solid line shows the boundaries derived by
Boyer et al. (2019), which we adjust slightly blueward (dotted line) to
include a handful of additional stars that appear to be an extension of the
branch of C-type stars. The red arrow shows the direction and magnitude of
circumstellar extinction for $E(J-K_{\rm s})=1$ mag for a mix of 70% amorphous
carbon and 30% SiC dust (Groenewegen, 2006b).  Figure 3: JWST CMDs with HST-
selected TP-AGB stars overplotted. Light grey points are all stars, and dark
points are stars in the region where the WFC3/IR and NIRCam footprints
overlap. Blue squares and red diamonds are stars classified using the HST
medium-band photometry, which covers a smaller FOV inside the NIRCam footprint
(Fig. 1). While the separation of the TP-AGB subtypes is poor when only
considering the ERS SW (left) or LW (middle) bands in isolation, we find good
separation by combining both (right).

In Figure 2, we show the Boyer et al. (2019, B19) selection boundaries, which
recover more M-type stars than the Boyer et al. (2017, B17) boundaries. The
original B17 boundaries are based on stellar atmosphere models from Aringer et
al. (2009, 2016), which B19 adjusted slightly to better fit HST data of the
M31 galaxy. The exact boundaries between stellar types may be slightly
dependent on stellar metallicity and the age of the population. Here, we see
that a handful of stars around $\rm{F127M}-\rm{F139M}=0.25$ mag are excluded
from the B19 C-type selection, despite appearing to be an extension of the C
star branch. We therefore adjust the boundary slightly blueward to include
these stars (dotted line in Fig. 2). Inspection of Aringer’s library of C star
spectra indicates that the F127M$-$F139M interval between 0.2 and 0.25 mag can
be filled with C stars of relatively high effective temperature ($\sim$4000
K), which probably form only at low metallicities.

Figure 3 shows three NIRCam color-magnitude diagrams (CMDs) with the HST-
selected TP-AGB stars plotted in red and blue. It is clear that, with our
filters, a combination of SW and LW filters is needed for a robust separation
of C- and M-type stars since the SW-only (left) and LW-only (middle) CMDs show
the highest degree of overlap between the two types.

###  3.3 Selecting TP-AGB stars with JWST filters

The ERS SW filters (F090W and F150W) were selected to provide a good recovery
of the WLM SFH (McQuinn et al., 2023). The LW filters were observed
simultaneously at no additional cost; we selected F250M and F430M for the LW
channel, in part because the COLIBRI TP-AGB models (Marigo et al., 2013, 2017)
predicted that they would provide a reasonable separation between TP-AGB
stellar types when used together with the SW filters. The two LW filters
sample CO in the atmospheres of C-type stars and H2O in M-type stars (Fig. 4).
The F430M filter is also sensitive to emission from circumstellar dust. On the
SW side, the F150W filter samples H2O in M-type stars, while both filters
sample CN$+$C2 in C-type stars. See Figure 4 for example model atmosphere
spectra of C- and M-type TP-AGB stars from Aringer et al. (2009, 2016), which
are used to calculate the bolometric corrections applied to TP-AGB stars in
the COLIBRI tracks. In this section, we use the HST classifications and the
COLIBRI models to guide the JWST TP-AGB classification.

Figure 4: Model stellar atmospheres (bottom panel) of C- and M-type stars from
Aringer et al. (2009, 2016). The ERS NIRCam filters are shown in blue in the
top panel. Medium-band HST WFC3/IR filters are shown in gray and hatched
(F127M, F139M, F153M). Additional NIRCam filters discussed in §4.2.1 are
plotted in yellow. The dashed and dotted lines in the bottom panel indicate
the center wavelengths of the filters. The most prominent molecular features
are marked in the bottom panel for C-type stars (red) and M-type stars (blue).
Figure 5: CMDs showing which stars are brighter/fainter than the TRGB at
F090W, F150W, and F250M. Green and Yellow points are those that are brighter
than the TRGB in only F150W or F250M. Green points are those with
F090W$-$F150W $<0.8$ mag and include contaminating main sequence and blue
helium-burning stars. Yellow points are dominated by AGB stars and RSGs. Black
points are remaining stars that are brighter than the TRGB in any of the four
ERS filters, and include mostly contaminating sources (main sequence and blue
helium-burning stars). Horizontal lines mark the TRGB in the F090W (left),
F150W (middle), and F250M (right) filters. To maximize the evolved star sample
while minimizing contaminating sources, our classification scheme starts with
stars that are brighter than either F150W or F250M.

####  3.3.1 The effect of the JWST TRGB

The TRGB limit provides a magnitude cut in each band to eliminate RGB stars
and other faint objects. Here, we investigate which TRGB limits provide the
most complete sample of TP-AGB stars with minimal contamination. As mentioned
in §3.1, TP-AGB stars can sometimes fall below the TRGB due to variations
during the pulsation cycle, molecular absorption, and/or circumstellar dust
attenuation, with the strength of the effect depending on the wavelength.
Figure 5 shows that some dusty AGB stars fall below the F090W and F150W TRGBs
(the yellow points with the reddest colors in the CMDs). In addition, a large
fraction of blue contaminating stars, such as main sequence and blue He-
burning stars, are brighter than the F090W TRGB.

Both of these issues are mitigated by instead using a TRGB measured in longer
wavelengths. The right panel in Figure 5 shows that the F250M TRGB minimizes
contamination from blue sources while also including dusty AGB stars. We
therefore compile our starting TP-AGB sample by including stars that are
either brighter than the F150W or F250M TRGBs. We exclude the F430M TRGB
because very red contaminating sources such as unresolved background galaxies
and embedded young stellar objects can be brighter than the TRGB at
wavelengths redder than about 3 µm (e.g., Jones et al., 2017).

Figure 6: Color-color diagrams used to classify C- and M-type TP-AGB stars.
The upper panels (A and B) show the JWST data, while the lower panels (C and
D) show simulations based on the COLIBRI models and the SFH from McQuinn et
al. (2023), scaled to match the number of AGB stars in the HST$+$NIRCam
footprint. In all panels, only stars that are brighter than either the F150W
or F250M TRGBs are included. Red diamonds are C-type stars and blue squares
are M-type stars, classified via HST colors (panels A and B) or identified via
their model parameters (panels C and D). Light pink and cyan dots in the lower
panels show the same simulation for 10$\times$ WLM’s stellar mass, to
illustrate more sparsely populated regions of these diagrams. The black lines
mark the JWST color selections (Eqs. 1 and 2). Stars of interest are outlined
in black circles, squares, and triangles, see text. In panel A, possible
contaminants 1 and 2 are marked (see §3.3.4).

####  3.3.2 JWST TP-AGB Classification

Based on the HST TP-AGB classifications for WLM, we find that the widest color
separation between stellar types is in the F090W$-$F150W vs. F150W$-$F250M
color-color diagram (CCD), shown in Figure 6, panel A. The CCD with
F150W$-$F430M on the y-axis (panel B) is also promising, but there is more
overlap between M-type stars (marked by large black squares or circles) and
the C-type stars. This overlap is concerning for systems that may have a
larger population of late-type M stars with deep H2O absorption. In Panel A,
these stars are more separated from the C star population. The black lines in
these panels mark our best separation between AGB subtypes based on the HST
classification, chosen to minimize cross-contamination between C- and M-type
stars in these JWST colors:

${\rm F250M:}\\\ (F090W-F150W,F150W-F250M)_{\rm HST{\text{-}}driven}\\\
=[(1.4,2.0),(1.4,0.36),(1.6,0.36),(3.0,1.0)],$ (1)

${\rm F430M:}\\\ (F090W-F150W,F150W-F430M)_{\rm HST{\text{-}}driven}=\\\
[(1.4,2.0),(1.4,0.6),(1.55,0.60),(2.75,2.0)].$ (2)

These AGB star selections are supported by comparing to model predictions.
Using the star formation history (SFH) from McQuinn et al. (2023), we simulate
the population of WLM using the TRILEGAL code (Girardi et al., 2005) and the
COLIBRI TP-AGB tracks from Pastorelli et al. (2019, 2020). The bolometric
corrections for AGB stars are based on the stellar spectral libraries by
Aringer et al. (2009) for C-rich stars, and Aringer et al. (2016) for O-rich
stars. The effect of circumstellar dust is taken into account following the
approach of Marigo et al. (2008) and adopting the following dust mixtures:
amorphous carbon (85%) and SiC (15%) for C-rich stars, and silicates for
O-rich stars (Groenewegen, 2006b). For non-TP-AGB stars, the simulation uses
tracks from the PARSEC v1.2S database (Bressan et al., 2012; Tang et al.,
2014; Chen et al., 2015).

The COLIBRI tracks are calibrated with TP-AGB stars from the Small Magellanic
Cloud (SMC) for $0.0005\leq Z_{\rm i}\leq 0.006$ and from the Large Magellanic
Cloud (LMC) for $0.008\leq Z_{\rm i}\leq 0.02$. The SMC is slightly more
metal-rich than WLM ($12+\log({\rm O/H})\sim$ 8–8.4, compared to 7.8 for WLM;
Lee et al., 2005; Toribio San Cipriano et al., 2017), and it is as yet unclear
how strongly this metallicity difference affects the AGB colors (see §4.2).
The simulation is scaled to match the total number of AGB stars in the
combined HST$+$JWST footprint (Fig. 1), and is shown in the lower panels of
Figure 6 with large symbols. The boundary between C and M stars remains
similar to the boundary derived using the HST classification (Eqs. 1 and 2).

The lower panels in Figure 6 also show the same simulation for 10$\times$
WLM’s stellar mass to highlight the more sparsely populated regions of the CCD
(C- and M-type stars are marked with light red and blue colors). This more
massive simulation shows that the JWST color cuts devised using the HST
classifications may miss a non-negligible number of C stars in a more
populated system. In fact, the massive simulation suggests that the F430M
color (panels B and D) may be the better choice for classifying AGB stars
overall. See §4 for more discussion.

Given WLM’s low stellar mass and sparsely populated M star branch, it is
difficult to know exactly how much overlap occurs between C-type and M-type
stars at metallicities lower than the SMC. If we look at the two HST-confirmed
M-type stars with the deepest H2O absorption (marked with large black squares
in Fig. 6), it appears as though the M star branch might curve upwards in the
F430M CCD, creating almost total overlap with the C star branch. We discuss
this point further in §4.2. The two stars marked with circles that also follow
this upwards trajectory are also HST-confirmed M-type stars, though they
appear to have less extreme H2O absorption, and may be red in F150W$-$F430M
due to circumstellar dust (see § 4.3). There is one star that is not in the
HST footprint that appears to follow the simulated M star branches well
(marked with a large black triangle), which may indicate that a more populated
galaxy would have more stars along this branch. On the other hand, this star
is a possible contaminant (§3.3.4). A larger JWST footprint in WLM and/or
observations in similarly metal-poor systems would help to verify the
behaviour of metal-poor M type stars.

####  3.3.3 Isolating Massive Red Supergiants

Figure 6 also shows several stars brighter than the TRGB that are not
classified as TP-AGB stars in either the HST data or the simulations
(black/grey points with F090W–F150W $\lesssim$ 1.3 mag). This population
includes a combination of massive RSGs (also referred to as red helium-burning
– or RHeB – stars), main sequence (MS) stars, blue helium burning (BHeB)
stars, and foreground stars.

We find that RSGs are most easily isolated in the F090W-F430M vs. F430M color-
magnitude diagram (Fig. 7), where they form a distinct branch blueward of the
M-type AGB stars, similar to what is seen in other near-IR surveys of nearby
star-forming galaxies (e.g., Nikolaev & Weinberg, 2000; Blum et al., 2006;
Boyer et al., 2011; Rosenfield et al., 2016). We follow Rosenfield et al.
(2016) to define a boundary for the RSGs that depends on the color of the
TRGB, which is itself dependent on the metallicity of the galaxy:

$\begin{split}M_{\rm F430M}=-7.4-4.5[&(M_{\rm F090W}-M_{\rm F430M})-\\\
&(M_{\rm F090W}^{\rm TRGB}-M_{\rm F430M}^{\rm TRGB})];\end{split}$ (3)

where we assume a true distance modulus of $\mu=24.93$ mag (Albers et al.,
2019). Note that this boundary is conservative, in that it minimizes the loss
of M-type AGB stars by keeping the dividing line closer to the RSG branch. It
is, however, likely that a small handful of RSGs fall in our M-type AGB
sample. Using Equation 3, we find 33 RSGs that are brighter than either the
F150W or F250M TRGB.

Points blueward of F090W-F430M = 1 mag are classified as MS and BHeB stars.
Foreground stars can also populate this color-magnitude space, though NIRCam’s
small FOV and WLM’s Galactic latitude ($b=-73.6^{\circ}$) together result in
minimal foreground in our sample.

Figure 7: F090W$-$F430M vs. F430M CMD showing the best separation of RSGs and
AGB stars in the ERS filters (dashed line; Eq. 3). Points to the left of the
dotted line are mostly likely MS stars or BHeB stars. AGB stars are color-
coded according to Equation 1 (C-type stars are red, M-type are dark blue).
Figure 8: JWST CMDs showing the extended objects that were cut from the final
point-source catalog in black. Most of these are likely background objects or
very red objects belonging to WLM, such as YSOs. Without the quality cuts, a
handful of these would be in our AGB sample since they are brighter than the
F250M TRGB (black line). Even with the quality cuts, two objects remain in our
sample, circled in black and marked 1 and 2.

####  3.3.4 Contamination

To produce our final point-source catalog, we applied several quality cuts to
eliminate objects that are slightly extended based on the DOLPHOT sharp,
crowd, flag, and object type parameters (Warfield et al., 2023; Weisz et al.,
2024). Figure 8 shows that a large fraction of the eliminated extended objects
are redward of the RGB. The Spitzer surveys of the SMC and LMC (Meixner et
al., 2006; Gordon et al., 2011) showed that these objects are dominated by
unresolved background galaxies, young stellar objects (YSOs), and rare objects
such as planetary nebulae (Boyer et al., 2011; Ruffle et al., 2015; Jones et
al., 2017). While our quality cuts eliminated the vast majority of these
objects, two remain in this region of the CMD that are classified as AGB stars
based on their colors/brightnesses, marked as 1 and 2 in Figure 8. Given their
positions below the extension of the C-type AGB branch in Figure 8, their
status as AGB stars is suspect. Figure 9 shows their JWST spectral energy
distributions (SEDs). Source 2 peaks near 1 $\mu$m, and Source 1 peaks
somewhere past 4.3 $\mu$m. Based on their SED shapes and positions in the CMD,
these could be AGB stars, or perhaps post-AGB stars, given that their SEDs are
similar to those in Fig. 6 from Ruffle et al. (2015). Spitzer observations did
not detect any variability (Boyer et al., 2015; Goldman et al., 2019). Longer
wavelength data and/or spectroscopy is needed to confirm their nature.

By requiring our TP-AGB and RSG candidates to be brighter than the F250M TRGB,
most interlopers are eliminated even without the DOLPHOT quality cuts.
However, we note that several of these slightly extended objects are brighter
than the TRGB at F430M (Fig. 8), and would be classified primarily as C-rich
AGB stars if they remained in our sample. We show their SEDs in Figure 9,
compared to some example dusty C-type AGB stars (F090W$-$F150W $>$ 2) from
this dataset. The SEDs of the two populations are mostly indistinguishable at
these wavelengths, and we cannot eliminate the possibility that objects in
this region of the CMD are evolved stars. However, since they fall below the
primary C-type AGB branch at F430M = 14–16 mag, it is likely that if they
belong to WLM, they are YSOs, planetary nebulae, or post-AGB stars. We exclude
these extended objects from our sample, but include sources 1 and 2 identified
in Figure 8.

Figure 9: Left: SEDs of the two sources highlighted in Figure 8. Middle: SEDs
of objects brighter than the F430M TRGB that were identified as extended based
on DOLPHOT quality cuts. Right: SEDs of the brightest and reddest C-type AGB
stars in our sample.  Table 2: Number of Evolved Stars

Classification | Field | $N_{\rm M}$ | $N_{\rm C}$ | $N_{\rm RSG}$
---|---|---|---|---
JWST | North NIRCam | 61 | 40 | 24
JWST | South NIRCam | 26 | 30 | 9
HST | HST FOV | 83 | 66 |  $\cdots$
HST | NIRCam FOV | 35 | 32 |  $\cdots$

Note. — For the JWST classifications, the AGB stars (C- and M-type) are based
on Equation 1 and the RSGs are based on Equation 3. In both cases, all stars
are brighter than either the F150W or F250M TRGB.

##  4 Discussion

###  4.1 Final Evolved Star Selection

Based on the HST classifications and the PARSEC-COLIBRI models, our preferred
selections are Equation 1 for TP-AGB stars and Equation 3 for RSGs, resulting
in the statistics listed in Table 2. This selection may minimize overlap
between C stars and late-type M stars. Out of a total of 157 TP-AGB stars in
WLM, only two are outliers and possibly contaminants based on their JWST
colors (stars 1 and 2 in §3.3.4). Only one star conflicts with the HST
classification, a carbon star marked with a black circle in Figures 6 and 10.
This star sits right on the edge of our JWST-based classification, but is
likely an M-type star based on its HST colors. The classification based on
F430M (Eq. 2) mis-identifies 4 stars: 3 M-type and 1 C star.

Our JWST classification scheme has also picked up 16 more M-type stars in the
HST footprints with $\rm{F139M-F153M}<0.3$ (Fig. 10). Two of these are clearly
outliers and are likely contaminants, based on their positions in the HST CCD
(see Boyer et al., 2017). According to the PARSEC-COLIBRI models, the rest of
the stars in this region (about 30% of the total M star sample) are a mix of
upper RGB stars and early-AGB stars. Based on their effective temperatures,
they may be K-type stars, but they cannot be excluded from our M-star sample
with our set of JWST filters.

Figure 10: JWST-classified TP-AGB stars (Eq. 1) shown on the HST CCD. The
large squares and circles indicate the stars marked with the same symbols in
Figure 6. Only one star conflicts with the HST classification: the red diamond
in the M star region of the CCD.

We can use the ratio of C to M stars (C/M) to compare the HST and JWST
classifications, which differ mostly by the additional M/K-type stars picked
up by the JWST classification. The C/M ratio is sometimes used as a proxy for
metallicity (e.g., Groenewegen, 2006a). Here, we use the relation between C/M
and the metallicity derived by Boyer et al. (2019) and assume a metallicity of
${\rm[M/H]}=-1.27\pm 0.04$, measured using medium resolution spectroscopy of
the Calcium II Triplet equivalent widths in RGB stars (Leaman et al., 2009).
Within the combined HST+JWST footprints (Fig. 1), the HST and JWST-derived C
to M ratios are ${\rm(C/M)_{HST}}=0.91\pm 0.16$ and ${\rm(C/M)_{JWST}}=0.68\pm
0.12$, yielding [M/H] values of $-1.23\pm 0.07$ for HST and $-1.11\pm 0.07$
for JWST. The HST value is clearly closer to the metallicity measured by
Leaman et al. (2009). However, low number statistics are playing a role here,
since the C/M value in the entire HST footprint, derived using the HST
classifications is $0.80\pm 0.10$, while the overall JWST C/M in the entire
NIRCam footprint is also $0.80\pm 0.10$. These values yield a metallicity of
${\rm[M/H]}=-1.18\pm 0.05$, within one sigma of the Leaman value.

Moreover, all of these C/M values are within the uncertainties predicted from
C/M values from several nearby dwarf galaxies, which were derived using the
HST classification scheme and an optical classification technique that uses
narrow band CN and TiO filters (see Fig. 9 from Boyer et al., 2019). We
conclude that while the HST and JWST classification schemes give slightly
different results, the differences result primarily from low number
statistics.

Figure 11: Simulations using a constant star-formation rate at two
metallicities, showing only stars brighter than the TRGB. The insets in the
left hand panels are zoomed into the region outlined by a black box. At low
metallicity (left), it is difficult to isolate the bluest C stars with F250M,
and difficult to isolate dusty M-type AGB stars with F430M. At high
metallicity (right), either filter works across the entire F090W$-$F150W color
range.

###  4.2 Metallicity Effect on TP-AGB colors

As noted in §3.3.2, the COLIBRI models that we compare to here are calibrated
to the Magellanic Clouds, which are more metal-rich than WLM. In this section,
we investigate the effect metallicity has on the models and on the JWST
colors. Figure 11 shows simulations from COLIBRI+TRILEGAL, this time with a
constant SFH to ensure all TP-AGB ages are well populated. We simulated two
metallicities, one at SMC metallcity ($Z=0.004)$, and one at high metallicity
($Z=0.014$, about 70% solar). The higher-metallicity COLIBRI models are
currently calibrated to the LMC (Pastorelli et al., 2020), at $Z=0.008$. The
most dramatic differences between these simulations is the color of the M-type
stars on the x-axes and the bluest extent of the C stars. These differences
are presumably caused by the cooler temperatures and deeper molecular
absorption at high metallicity, as well as the narrower age range of C stars
at high metallicity.

As the metallicity increases, the separation between C- and M-type stars
improves significantly. At SMC metallicity, the different stellar types have
quite a narrow separation, made even narrower by the presence of many dusty
M-type stars, which appear to be rare in WLM. This metallicity dependence
suggests that in lower metallicity systems like WLM, we should expect even
more overlap between the C- and M-type stars. Indeed, this does appear to be
the case, as illustrated in the upper panels of Figure 6. However, it is
difficult to discern how severe the overlap is, given the small number of
late-type M stars in WLM. In a more populated metal-poor system, it is
possible that the JWST colors used by this ERS program would not provide a
good separation between AGB subtypes. On the other hand, late-type M stars are
quite rare in metal-poor systems in general, so this overlap may not have a
significant effect on the overall statistics. Nevertheless, we advise caution
when using the ERS filters (F090W, F150W, F250M, and F430M) to classify AGB
stars in very metal-poor systems.

####  4.2.1 Other JWST Filters for Evolved Star Classification

Other JWST filters sample interesting features in TP-AGB and RSG stars. For
example, F140M, F182M, and F210M sample H2O, F300M, F335M, F410M, and F430M
sample HCN$+$C2H2. Despite this, we were unable to find filter combinations
that effectively isolate AGB stellar types from one another, independent of
metallicity. Using the simulations with constant star formation rate described
above, we tested every possible combination, and found that most combinations
are either worse or similar to the ERS filters in their classification
effectiveness. Only a handful of options show some marginal improvement, shown
in Figure 12. All of these options include F210M and/or F430M (see Fig. 4).
The F210M$-$F430M color specifically has the widest separation between the
bulk of the M-type stars and C-type stars, though circumstellar dust can
complicate matters. The upper and lower panels in Figure 12 show combinations
that are also quite good at separating RSGs (in grey) from the TP-AGB stars.

As with the ERS filters discussed in §4.2, it is possible that the C and M
stars will show more overlap at lower metallicities. However, the HST medium-
band filters appear to work well over a wide range of metallicity (Boyer et
al., 2017, 2019). If a clean sample of these stellar types in a metal-poor
galaxy is essential, the HST medium-band WFC3/IR filters may be the better
option, despite HST’s decreased sensitivity and angular resolution in the
near-IR. We are providing the simulations shown in Figure 12 along with this
manuscript, in every NIRCam filter, to assist in future JWST program
designs.333 We plan to post the simulations to CDS and will link to them in
this footnote.

Figure 12: Additional JWST filters that show good separation between AGB
subtypes, according to the same constant star-formation rate simulation shown
in Figure 11. For reference, these filters are also plotted in Figure 4. The
F210M$-$F430M shows the widest separation between the bulk of the C- and
M-type AGB stars. The combinations shown in the upper and lower panels also
show good isolation of the RSGs (grey). Only stars brighter than the TRGB are
plotted here.

###  4.3 Dust

With JWST, we find 10 very dusty AGB star candidates in WLM (F150W$-$F430M $>$
2 mag), and many more slightly dusty candidates (F150W$-$F430M = 0.7–2 mag;
Fig. 3). With Spitzer, contamination from unresolved background galaxies made
it impossible to identify dusty objects belonging to WLM, so Spitzer surveys
relied on variability information with 2–8 epochs of data (Boyer et al., 2015;
Goldman et al., 2019). With just two epochs of data, Boyer et al. (2015) were
able to identify 9/10 of the JWST very dusty AGB stars. However, even with 8
epochs, Goldman et al. (2019) could not identify any of the less-dusty AGB
stars, presumably because their pulsation amplitudes are too small. JWST data
are thus clearly essential to obtaining a full picture of the dust budget in
galaxies more distant than the Magellanic Clouds.

All of the dusty stars with F150W$-$F430M $>$ 2 mag are identified as carbon
stars using the classification described in §4.1. Of the less dusty stars, 11%
are M-type. While we cannot measure the dust mass produced by these stars
without data beyond 5 µm, it is clear that the present-day dust input by
evolved stars in WLM is dominated by carbonaceous dust, and thus by stars with
masses between about 1.5–3 $M_{\odot}$ (Marigo et al., 2013).

Dusty massive M-type stars such as those seen in the extremely metal-poor
galaxies Sextans A and Sag DIG (Boyer et al., 2017) are of particular interest
to the evolution of dust in the early Universe since they can start producing
dust as early as 30 Myr after they form. However, based on their brightnesses,
the most massive dusty M-type AGB stars that we see in WLM are likely
$\sim$2.5–3.5 $M_{\odot}$, which corresponds to a dust input timescale of
about 500 Myr – too late to be proxies of dusty AGB stars in the early
Universe. Based on their colors and magnitudes, it is also possible that two
of the dusty M-type AGB candidates are actually RSGs, which is possible
because WLM has ongoing star formation and evidence for very recent star
formation (Mondal et al., 2018). If so, that would put their masses between
about 8 and 23 $M_{\odot}$, and currently existing dust would be at least
partially destroyed by their future supernova explosions.

##  5 Conclusions

We imaged WLM in four NIRCam filters as part of the JWST Resolved Stellar
Populations ERS program. These filters were selected, in part, to investigate
how effectively different subtypes of evolved stars can be identified with
JWST, specifically TP-AGB stars and RSGs. In a single NIRCam footprint, we
find 70 carbon TP-AGB stars, 87 M-type TP-AGB stars, and 33 RSGs. Of these, 72
(or about 46%) show evidence of circumstellar dust.

We find that the ERS filters (F090W, F150W, F250M, and F430M) are excellent
for isolating RSGs, but the efficacy of the AGB selection is likely highly
dependent on metallicity. A comparison to the well-established HST medium-band
AGB classification strategy shows very good agreement to the ERS
classification defined here for WLM, resulting in an identical C/M ratio
derived from both techniques. This C/M ratio (0.8$\pm$0.1) follows the
relationship of C/M versus metallicity derived using HST medium-band data of
other nearby dwarf galaxies.

However, WLM is a low-mass system and short-lived stellar types such as late-
type M AGB stars and dusty AGB stars are rare. Metal-poor systems with more
massive populations may suffer from significant overlap between these rarer
objects and carbon-rich AGB stars using these filters. The COLIBRI stellar
evolution models suggest that the ERS filters will be more effective in more
metal-rich systems, where cooler temperatures and a narrowing mass range for
carbon stars leads to a much wider separation of AGB subtypes. We also
investigated other NIRCam filters and found that combinations using F182M,
F210M, and F460M may provide some improvements at low metallicity.

We also compared WLM’s dusty AGB population to that identified with Spitzer by
Boyer et al. (2015). We find that Spitzer was only able to identify the 9 most
extreme examples out of 71 total dusty AGB stars in the NIRCam footprint. Of
these dusty stars, 89% are carbon-rich, confirming that carbonaceous dust
dominates the dust input into WLM’s interstellar medium. While dusty objects
are typically studied in the mid-IR, these results demonstrate NIRCam’s
ability to efficiently identify dust-producing stars photometrically, while
also providing information on the basic dust species (carbonaceous or
siliceous). NIRCam is thus a valuable tool for the identification of dust-
producing stars in much more distant systems than can be reached by JWST’s
Mid-InfraRed Instrument (MIRI), which has lower spatial resolution and
sensitivity.

This work is based in part on observations made with the NASA/ESA/CSA JWST.
The data were obtained from the Mikulski Archive for Space Telescopes at the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127
for JWST. These observations are associated with program ERS-1334. This
research is also based in part on observations made with the NASA/ESA Hubble
Space Telescope obtained from the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS 5–26555. These observations are associated with
program HST-GO-16162. Support for this work was provided by NASA through
grants ERS-1334 from the Space Telescope Science Institute, which is operated
by AURA, Inc., under NASA contract NAS5-26555. This research has made use of
NASA Astrophysics Data System Bibliographic Services and the NASA/IPAC
Extragalactic Database (NED), which is operated by the Jet Propulsion
Laboratory, California Institute of Technology, under contract with the
National Aeronautics and Space Administration. E.N.K. acknowledges support
from NSF CAREER grant AST-2233781. We thank the JWST teams at STScI for a
smooth execution of this ERS program and for assistance in its design and
calibration.

## References

  * Albers et al. (2019) Albers, S. M., Weisz, D. R., Cole, A. A., et al. 2019, MNRAS, 490, 5538, doi: 10.1093/mnras/stz2903
  * Anand et al. (2024) Anand, G. S., Riess, A. G., Yuan, W., et al. 2024, arXiv e-prints, arXiv:2401.04776, doi: 10.48550/arXiv.2401.04776
  * Aringer et al. (2016) Aringer, B., Girardi, L., Nowotny, W., Marigo, P., & Bressan, A. 2016, MNRAS, 457, 3611, doi: 10.1093/mnras/stw222
  * Aringer et al. (2009) Aringer, B., Girardi, L., Nowotny, W., Marigo, P., & Lederer, M. T. 2009, A&A, 503, 913, doi: 10.1051/0004-6361/200911703
  * 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
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Baldwin et al. (2018) Baldwin, C., McDermid, R. M., Kuntschner, H., Maraston, C., & Conroy, C. 2018, MNRAS, 473, 4698, doi: 10.1093/mnras/stx2502
  * Blum et al. (2006) Blum, R. D., Mould, J. R., Olsen, K. A., et al. 2006, AJ, 132, 2034, doi: 10.1086/508227
  * Boyer et al. (2011) Boyer, M. L., Srinivasan, S., van Loon, J. T., et al. 2011, AJ, 142, 103, doi: 10.1088/0004-6256/142/4/103
  * Boyer et al. (2013) Boyer, M. L., Girardi, L., Marigo, P., et al. 2013, ApJ, 774, 83, doi: 10.1088/0004-637X/774/1/83
  * Boyer et al. (2015) Boyer, M. L., McQuinn, K. B. W., Barmby, P., et al. 2015, ApJ, 800, 51, doi: 10.1088/0004-637X/800/1/51
  * Boyer et al. (2017) Boyer, M. L., McQuinn, K. B. W., Groenewegen, M. A. T., et al. 2017, ApJ, 851, 152, doi: 10.3847/1538-4357/aa9892
  * Boyer et al. (2019) Boyer, M. L., Williams, B. F., Aringer, B., et al. 2019, ApJ, 879, 109, doi: 10.3847/1538-4357/ab24e2
  * Boyer et al. (2022) Boyer, M. L., Anderson, J., Gennaro, M., et al. 2022, Research Notes of the American Astronomical Society, 6, 191, doi: 10.3847/2515-5172/ac923a
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Chen et al. (2015) Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068, doi: 10.1093/mnras/stv1281
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Dalcanton et al. (2012) Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18, doi: 10.1088/0067-0049/200/2/18
  * Dolphin (2016) Dolphin, A. 2016, DOLPHOT: Stellar photometry, Astrophysics Source Code Library, record ascl:1608.013.  http://ascl.net/1608.013
  * Dolphin (2000) Dolphin, A. E. 2000, PASP, 112, 1383, doi: 10.1086/316630
  * Freedman & Madore (2020) Freedman, W. L., & Madore, B. F. 2020, ApJ, 899, 67, doi: 10.3847/1538-4357/aba9d8
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001, doi: 10.1088/1538-3873/acd1b5
  * Girardi et al. (2005) Girardi, L., Groenewegen, M. A. T., Hatziminaoglou, E., & da Costa, L. 2005, A&A, 436, 895, doi: 10.1051/0004-6361:20042352
  * Girardi et al. (2010) Girardi, L., Williams, B. F., Gilbert, K. M., et al. 2010, ApJ, 724, 1030, doi: 10.1088/0004-637X/724/2/1030
  * Goldman et al. (2019) Goldman, S. R., Boyer, M. L., McQuinn, K. B. W., et al. 2019, ApJ, 884, 152, doi: 10.3847/1538-4357/ab418a
  * Gordon et al. (2011) Gordon, K. D., Meixner, M., Meade, M. R., et al. 2011, AJ, 142, 102, doi: 10.1088/0004-6256/142/4/102
  * Groenewegen (2006a) Groenewegen, M. A. T. 2006a, in Planetary Nebulae Beyond the Milky Way, ed. L. Stanghellini, J. R. Walsh, & N. G. Douglas, 108, doi: 10.1007/3-540-34270-2_13
  * Groenewegen (2006b) Groenewegen, M. A. T. 2006b, A&A, 448, 181, doi: 10.1051/0004-6361:20054163
  * Habing & Olofsson (2004) Habing, H. J., & Olofsson, H. 2004, Asymptotic Giant Branch Stars, doi: 10.1007/978-1-4757-3876-6
  * Jones et al. (2023) Jones, O. C., Boyer, M. L., McDonald, I., Meixner, M., & van Loon, J. T. 2023, arXiv e-prints, arXiv:2301.10613, doi: 10.48550/arXiv.2301.10613
  * Jones et al. (2017) Jones, O. C., Woods, P. M., Kemper, F., et al. 2017, MNRAS, 470, 3250, doi: 10.1093/mnras/stx1101
  * Karakas (2017) Karakas, A. I. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 461, doi: 10.1007/978-3-319-21846-5_117
  * Leaman et al. (2009) Leaman, R., Cole, A. A., Venn, K. A., et al. 2009, ApJ, 699, 1, doi: 10.1088/0004-637X/699/1/1
  * Leaman et al. (2012) Leaman, R., Venn, K. A., Brooks, A. M., et al. 2012, ApJ, 750, 33, doi: 10.1088/0004-637X/750/1/33
  * Leaman et al. (2013) —. 2013, ApJ, 767, 131, doi: 10.1088/0004-637X/767/2/131
  * Lee (2023) Lee, A. J. 2023, ApJ, 956, 15, doi: 10.3847/1538-4357/acee69
  * Lee et al. (2005) Lee, H., Skillman, E. D., & Venn, K. A. 2005, ApJ, 620, 223, doi: 10.1086/427019
  * Marigo et al. (2013) Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488, doi: 10.1093/mnras/stt1034
  * Marigo et al. (2008) Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883, doi: 10.1051/0004-6361:20078467
  * Marigo et al. (2017) —. 2017, ApJ, 835, 77, doi: 10.3847/1538-4357/835/1/77
  * McQuinn et al. (2023) McQuinn, K. B. W., Newman, M. J. B., Savino, A., et al. 2023, arXiv e-prints, arXiv:2312.03060, doi: 10.48550/arXiv.2312.03060
  * Meixner et al. (2006) Meixner, M., Gordon, K. D., Indebetouw, R., et al. 2006, AJ, 132, 2268, doi: 10.1086/508185
  * Melotte (1926) Melotte, P. J. 1926, MNRAS, 86, 636, doi: 10.1093/mnras/86.8.636
  * Méndez et al. (2002) Méndez, B., Davis, M., Moustakas, J., et al. 2002, AJ, 124, 213, doi: 10.1086/341168
  * Mondal et al. (2018) Mondal, C., Subramaniam, A., & George, K. 2018, AJ, 156, 109, doi: 10.3847/1538-3881/aad4f6
  * Nikolaev & Weinberg (2000) Nikolaev, S., & Weinberg, M. D. 2000, ApJ, 542, 804, doi: 10.1086/317048
  * O’Donnell (1994) O’Donnell, J. E. 1994, ApJ, 422, 158, doi: 10.1086/173713
  * Pastorelli et al. (2019) Pastorelli, G., Marigo, P., Girardi, L., et al. 2019, MNRAS, 485, 5666, doi: 10.1093/mnras/stz725
  * Pastorelli et al. (2020) —. 2020, MNRAS, 498, 3283, doi: 10.1093/mnras/staa2565
  * Rieke et al. (2005) Rieke, M. J., Kelly, D., & Horner, S. 2005, in SPIE Conference Series, Vol. 5904, Cryogenic Optical Systems and Instruments XI, ed. J. B. Heaney & L. G. Burriesci, 1–8, doi: 10.1117/12.615554
  * Rieke et al. (2023) Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001, doi: 10.1088/1538-3873/acac53
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001, doi: 10.1088/1538-3873/acb293
  * Rosenfield et al. (2016) Rosenfield, P., Marigo, P., Girardi, L., et al. 2016, ApJ, 822, 73, doi: 10.3847/0004-637X/822/2/73
  * Rosenfield et al. (2014) —. 2014, ApJ, 790, 22, doi: 10.1088/0004-637X/790/1/22
  * Ruffle et al. (2015) Ruffle, P. M. E., Kemper, F., Jones, O. C., et al. 2015, MNRAS, 451, 3504, doi: 10.1093/mnras/stv1106
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Tang et al. (2014) Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287, doi: 10.1093/mnras/stu2029
  * Toribio San Cipriano et al. (2017) Toribio San Cipriano, L., Domínguez-Guzmán, G., Esteban, C., et al. 2017, MNRAS, 467, 3759, doi: 10.1093/mnras/stx328
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Warfield et al. (2023) Warfield, J. T., Richstein, H., Kallivayalil, N., et al. 2023, Research Notes of the American Astronomical Society, 7, 23, doi: 10.3847/2515-5172/acbb72
  * Weisz et al. (2014) Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014, ApJ, 789, 147, doi: 10.1088/0004-637X/789/2/147
  * Weisz et al. (2023) Weisz, D. R., McQuinn, K. B. W., Savino, A., et al. 2023, ApJS, 268, 15, doi: 10.3847/1538-4365/acdcfd
  * Weisz et al. (2024) Weisz et al., D. R. 2024, ApJ, in prep
  * Williams et al. (2014) Williams, B. F., Lang, D., Dalcanton, J. J., et al. 2014, ApJS, 215, 9, doi: 10.1088/0067-0049/215/1/9
  * Wolf (1909) Wolf, M. 1909, Astronomische Nachrichten, 183, 187, doi: 10.1002/asna.19091831204