# Coexisting Magnetism, Ferroelectric, and Ferrovalley Multiferroic in
Stacking-Dependent Two-Dimensional Materials

Wei Xun  State Key Laboratory for Mechanical Behavior of Materials, Center for
Spintronics and Quantum System, School of Materials Science and Engineering,
Xi’an Jiaotong University, Xi’an, Shaanxi, 710049, China  Chao Wu  State Key
Laboratory for Mechanical Behavior of Materials, Center for Spintronics and
Quantum System, School of Materials Science and Engineering, Xi’an Jiaotong
University, Xi’an, Shaanxi, 710049, China  Hanbo Sun  State Key Laboratory for
Mechanical Behavior of Materials, Center for Spintronics and Quantum System,
School of Materials Science and Engineering, Xi’an Jiaotong University, Xi’an,
Shaanxi, 710049, China  Weixi Zhang  Department of Physics and Electronic
Engineering, Tongren University, Tongren 554300, China  Yin-Zhong Wu  School
of Physical Science and Technology, Suzhou University of Science and
Technology, Suzhou 215009, China  Ping Li<EMAIL_ADDRESS>State Key
Laboratory for Mechanical Behavior of Materials, Center for Spintronics and
Quantum System, School of Materials Science and Engineering, Xi’an Jiaotong
University, Xi’an, Shaanxi, 710049, China

###### Abstract

The two-dimensional (2D) multiferroic materials have widespread of application
prospects in facilitating the integration and miniaturization of nanodevices.
However, it is rarely coupling between the magnetic, ferroelectric, and
ferrovalley in one 2D material. Here, we propose a mechanism for manipulating
magnetism, ferroelectric, and valley polarization by interlayer sliding in 2D
bilayer material. Monolayer GdI2 exhibits a ferromagnetic semiconductor with
the valley polarization up to 155.5 meV. More interestingly, the magnetism and
valley polarization of bilayer GdI2 can be strongly coupled by sliding
ferroelectricity, appearing these tunable and reversible. In addition, we
uncover the microscopic mechanism of magnetic phase transition by spin
Hamiltonian and electron hopping between layers. Our findings offer a new
direction for investigating 2D multiferroic in the implication for next-
generation electronic, valleytronic, and spintronic devices.

###### keywords:

Two-dimensional multiferroic, Sliding ferroelectricity, Ferrovalley, Magnetic
phase transition, d-orbital hopping

Faculty of Electronic Information Engineering, Huaiyin Institute of
Technology, Huaian 223003, China State Key Laboratory for Surface Physics and
Department of Physics, Fudan University, Shanghai, 200433, China

Two-dimensional (2D) van der Waals (vdW) materials are burgeoning as one of
the top candidates due to their various structural, miniaturized
dimensionality, electronic engineering degrees of freedom, and they have
highly tunable magnetic, electronic, and optical properties1, 2, 3, 4, 5, 6.
The interlayer weak vdW interaction also unleashes the flexibility to affect
physical properties by vdW stacking. Recently, sliding ferroelectricity
originating from polar stacking of nonpolar monolayers transition metal
dichalcogenides (TMDs)7, 8, 9, and interlayer sliding in AlN and BN
bilayers10, 11. Besides, stacking orders effectively regulate the magnetic
ground state, resulting in the magnetic phase transition, such as the CrI312
and CrBr313 bilayers. More interestingly, the topological states can be tuned
by stacking order in MnBr3 bilayer14. Through interlayer sliding of 2D
materials, Ma et al. proved in theory that MnBi2Te415, Co2CF216, FeCl217,
etc., can realize layer-polarized anomalous Hall effect (LPAHE). Additionally,
stacking order can effectively tune the intrinsic valley polarization of
ferrovalley materials, such as the YI218 and VSiGeP419. These findings suggest
that the stacking order has a profound effect on determining the crystal
symmetry, reallocating spin and charge between neighboring layers, and thus
regulating the strength of ferroelectric, magnetic, topological, and valley
polarization properties.

Multiferroics, materials exhibiting coupled two or more ferroic orders (i.e.,
ferromagnetic, ferroelectric, ferroelastic, ferrovalley), are especially
interesting for the next generation electronic devices. At present, 2D
multiferroic systems are often multiferroic heterostructures, such as
Cr2Ge2Te6/In2Se320, LaCl/In2Se321, and Cr2COOH/Sc2CO222. However, the
heterostructure system will increase the difficulty of fabrication device in
an experiment, which is disadvantage to wider application. In contrast, if the
2D multiferroic can be realized by vdW stacking order in the same kind of
material, it becomes more attractive. It will realize highly flexible and
controllable coupling of the multiple ferroic orders. Therefore, the 2D
multiferroic materials will further promote the application prospect in
nanodevices.

In this work, we propose a new design of 2D multiferroics including magnetism,
ferroelectric and ferrovalley. Here we focus on the coupling and phase
transition of magnetic, ferroelectric, and ferrovalley by the effect of
stacking orders in bilayer GdI2. Firstly, we theoretically proposed that
monolayer GdI2 shows the spontaneously sizable valley polarization, which is
manipulated by switching magnetization. Then, we found that bilayer GdI2 endow
the coexistence magnetism, ferroelectric and ferrovalley, supporting the
designed target for manipulating magnetism and valley polarization via
ferroelectric switching by interlayer sliding. Moreover, we reveal the
microscopic mechanism of magnetic phase transition by spin Hamiltonian and
electron hopping between layers. The highly flexible tunable multiferroic in
bilayer GdI2 offers a practical way for designing advanced valleytronic and
spintronic devices on account of the couplings between multiferroic orders.

The 2D multiferroic systems with the coexistence of magnetism, ferroelectric,
and ferrovalley are established. Practically, for a monolayer system with
intrinsically spontaneous valley polarization, it can realize multiferroic by
bilayer stacking and interlay sliding. As shown in Figure 1, it shows
realizing the mechanism of 2D multiferroic with the coexistence of magnetism,
ferroelectric, and ferrovalley. In AA stacking, the magnetic ground state is
antiferromagnetic (AFM) coupling. Simultaneously, there are no spontaneous
valley polarization and ferroelectric polarization in the system. When the AA
stacking slide to AB or BA stacking, it will occur the magnetic phase
transition from the AFM state to ferromagnetic (FM) state, and inducing the
spontaneous valley polarization and ferroelectric polarization (signed by dark
gray arrows). Moreover, the AB and BA stacking can be transformed into each
other by interlayer sliding, leading to switch the ferroelectric polarization
between +z and -z axis direction.

Figure 1: Illustration of the mechanism of coexisting magnetism,
ferroelectric, and ferrovalley multiferroic. (a) The magnetic ground state of
AA stacking bilayer lattice is AFM state, and without spontaneous valley
polarization. The AA stacking slide to AB stacking (b) or BA stacking (c), the
magnetic ground state will switch to FM state. Simultaneously, spontaneous
valley polarization will be realized. The valley polarization can also be
manipulated via ferroelectric switching in AB stacking (b) and BA stacking
(c). Orange and light green cones represent spin up and spin down bands,
respectively. Dark gray arrows denote ferroelectric polarization P.  Figure 2:
(a) Crystal structures of monolayer GdI2 from the top view and side view. (b)
Electron localization function of monolayer GdI2. (c) The calculated phonon
spectrum along the high-symmetry. (d) Spin-polarized band structure of
monolayer GdI2. (e) The orbital-resolved band structure with considering SOC.
(f, g) The band structures with considering SOC of magnetic moments along +z
(f) and -z (g) axis directions, respectively.

As shown in Figure 2(a), it exhibits the crystal structure of monolayer GdI2.
The inversion symmetry is broken in monolayer GdI2. The in-plane lattice
constant for the fully optimized is 4.17 ${\rm\AA}$. Besides, in order to
study the bonding characteristics of monolayer GdI2, we calculate the electron
localization function (ELF). As shown in Figure 2(b), the electrons are mainly
localized around the Gd and I atoms, while those are negligible between the
atoms, showing a typical ionic bonding for all the bonds. In addition, to
evaluate the stability of monolayer GdI2, we calculate the phonon dispersion
spectrum. As shown in Figure 2(c), the absence of imaginary modes along the
high-symmetry lines confirms the dynamical stability of monolayer GdI2.
Moreover, we calculate the formation energy that is -11.140 eV. It indicates
that GdI2 is easily prepared.

In order to make sure that the magnetic ground state of monolayer GdI2, the
rectangle supercell is constructed (see Figure S1). The energy difference
$\Delta E$ = EFM \- EAFM between FM and AFM state is found to be -266 meV,
meaning that the magnetic ground state of monolayer GdI2 is FM state. However,
the successful preparation of 2D magnetic Cr2Ge2Te6 3 and CrI3 4 reveal that
the magnetic anisotropy energy (MAE) plays an important role in the stability
of magnetic order. Therefore, we investigate the MAE, it is defined as the
energy difference MAE = E001 \- E100 between the magnetic moment along [001]
and [100]. The MAE is 0.777 meV per unit cell, standing for the magnetization
along the x axis. For the monolayer GdI2, the magnetic field (H =
-$\frac{4}{3}\frac{K}{\mu_{0}M}$) of 3.58 T is needed to tune the direction of
magnetization from the x to z axis.

Next, we concern on the band structures and associated valley properties of
monolayer GdI2. As shown in Figure 2(d), it is an indirect semiconductor in
the absence of spin-orbit coupling (SOC), the valence band maximum (VBM)
belong to the spin up band, located at the high symmetry K and K’ points.
While the conduction band minimum (CBM) is contributed by the spin down band,
located at the high symmetry $\Gamma$ point. Note that the K and K’ valleys of
VBM are energetic degeneracy. When the SOC is switched on, as shown in Figure
2(e, f), the degeneracy between K and K’ valleys are broken. The valley
splitting is 155.5 meV in monolayer GdI2, which is larger than VSe2 (90 meV)
23, 24, VSiXN4 ($\thicksim$ 70 meV) 25, MoTe2/EuO ($\thicksim$ 20 meV) 26, and
other ferrovalley materials. More fascinatingly, when the magnetization
direction is switched from the +z axis to -z axis, the valley polarization can
be effectively tuned [see Figure 2(g)]. From the orbital-resolved band
structure [see Figure 2(e) and Figure S2], the VBM bands are mainly dominated
by Gd dxy and d${}_{x^{2}-y^{2}}$ orbitals, while the CBM bands are
contributed by Gd dxy, d${}_{x^{2}-y^{2}}$, and d${}_{z^{2}}$ orbitals.
Moreover, we calculate the valley polarization of U values from 2.5 eV to 5
eV. As shown in Figure S3, We can clearly observe that the valley polarization
increases first and then decreases. However, the range of variation is very
small only 0.5 meV.

Figure 3:  Sliding energy barrier between upper and lower layers of (a) AA
stacking and (d) AA’ stacking bilayer GdI2. The gray and blue ball indicates
that the magnetic ground state is AFM and FM state, respectively. The orange
ball represents the energy minima. The energy difference between interlayer
AFM and FM states as a function of interlayer translation along [1$\bar{1}$0]
(solid red line) and [100] (dotted blue line) direction for AA stacking (b)
and AA’ stacking (e) bilayer GdI2. The full space of lateral shifts for AA
stacking (c) and AA’ stacking (f) bilayer GdI2. The positive or negative
regions represent FM or AFM state, respectively.

To understand the microscopic mechanism in the valley splitting of valence
bands, we use an effective Hamiltonian model to describe the physical nature
of valley polarization induced by the SOC effect. Taking the SOC effect as the
perturbation term,

$\hat{H}_{SOC}=\lambda\hat{L}\cdot\hat{S}=\hat{H}_{SOC}^{0}+\hat{H}_{SOC}^{1},$
(1)

where $\hat{L}$ and $\hat{S}$ are orbital angular and spin angular operators,
respectively. $\hat{H}_{SOC}^{0}$ denotes the interaction between the same
spin states, while $\hat{H}_{SOC}^{1}$ represents the interaction between the
opposite spin states. Since the VBM is contributed by spin up band, therefore,
the term $\hat{H}_{SOC}^{1}$ can be neglected. For the $\hat{H}_{SOC}^{0}$, it
can be written by the polar angles 27

$\hat{H}_{SOC}^{0}=\lambda\hat{S}_{z^{\prime}}(\hat{L}_{z}cos\theta+\frac{1}{2}\hat{L}_{+}e^{-i\phi}sin\theta+\frac{1}{2}\hat{L}_{-}e^{+i\phi}sin\theta),$
(2)

When the magnetization direction of monolayer GdI2 along the out-of-plane,
$\theta$ = $\phi$ = 0∘, therefore, the $\hat{H}_{SOC}^{0}$ term can be
simplified as

$\hat{H}_{SOC}^{0}=\lambda\hat{S}_{z}\hat{L}_{z},$ (3)

Considering the C3 symmetry and the contribution of orbital, we chose
$|$$\psi$vτ$\rangle$=$\frac{1}{\sqrt{2}}$($|$dxy$\rangle$+i$\tau$$|$d${}_{x^{2}-y^{2}}$$\rangle$)$\otimes$$|$$\uparrow$$\rangle$,
as the orbital basis functions for the VBM, where $\tau$ = $\pm$1 represent
the valley index corresponding to $\rm K/\rm K^{\prime}$. The energy of the K
and K’ valleys for the VBM can be written as Evτ = $\langle$ $\psi$vτ $|$
$\hat{H}$${}_{SOC}^{0}$ $|$ $\psi$vτ $\rangle$. Then, the valley splitting in
the valence band can be written as

$E_{v}^{K}-E_{v}^{K^{\prime}}=i\langle
d_{xy}|\hat{H}_{SOC}^{0}|d_{x^{2}-y^{2}}\rangle-i\langle
d_{x^{2}-y^{2}}|\hat{H}_{SOC}^{0}|d_{xy}\rangle\approx 4\beta,$ (4)

where the $\hat{L}_{z}|d_{xy}\rangle$ = -2i$\hbar$$|d_{x^{2}-y^{2}}\rangle$,
$\hat{L}_{z}|d_{x^{2}-y^{2}}\rangle$ = 2i$\hbar$$|d_{xy}\rangle$, and
$\beta=\lambda\langle
d_{x^{2}-y^{2}}|\hat{S}_{z^{\prime}}|d_{x^{2}-y^{2}}\rangle$. The analytical
result proves that the valley splitting of valence band is consistent with our
DFT calculations ($E_{v}^{K}$ \- $E_{v}^{K^{\prime}}$ = 155.5 meV).

Layer stacking of 2D magnetic materials likely induces intriguing physical
phenomena, which differ from the monolayer materials 12, 13, 14, 15, 16, 17,
28. Here, we investigate novel physical behaviors of bilayer GdI2 by the
stacking orders. Six typical stacking structures, AA, AB, BA, AA’, AB’, and
BA’, are shown in Figure S4. The AA stacking bilayer GdI2 is constructed by
primitively placing one layer on top of the other layer. Therefore, AA
stacking layer GdI2 has mirror symmetry between the two layers about the xy
plane. AB and BA stacking are obtained by sliding t1//[$\frac{2}{3}$,
$\frac{1}{3}$, 0] and t1//[$\frac{1}{3}$, $\frac{2}{3}$, 0] of AA stacking,
respectively. However, the AA’ stacking can be realized by rotating the upper
layer 180∘ of AA stacking about the xy plane. AA’ stacking could be
transformed into the AB’ and BA’ stacking by by sliding t1//[$\frac{2}{3}$,
$\frac{1}{3}$, 0] and t1//[$\frac{1}{3}$, $\frac{2}{3}$, 0], respectively.

Table 1:  The interlayer distance of nearest-neighboring and second-neighboring between Gd atoms, the number of magnetic exchange interactions per unit cell in square brackets, the calculated Heisenberg exchange parameters, and the energy difference between interlayer AFM and FM states for the different stacking configurations.  Stacking | Gd-Gd NN (${\rm\AA}$) | Gd-Gd second-NN (${\rm\AA}$) | J1⊥ (meV) | J2⊥ (meV) | J∥ (meV) |  $\Delta E$ (meV)
---|---|---|---|---|---|---
AA | 8.139 [1] | 9.145 [6] | 0.243 | -0.003 | -8.253 | -1.794
AB | 7.888 [3] | 8.922 [3] | -0.036 | -0.007 | -8.242 | 1.046
AB’ | 7.505 [1] | 8.586 [6] | 0.002 | -0.004 | -8.235 | 0.206
BA’ | 7.915 [3] | 8.946 [3] | 0.066 | -0.008 | -8.256 | -1.403

The sliding energy barriers of two kinds of stacking patterns are investigated
to acquire the ground state stacking configurations. As shown in Figure 3(a),
for the one kind stackings, AB and BA stackings exhibit the lowest energy,
while the AA stacking has the highest energy. The transition AB stacking into
BA stacking requires overcoming an energy barrier of 16.42 meV. For the other
kinds of stacking [see Figure 3(d)], BA’ stacking shows the lowest energy,
while the energy of AB’ stacking is 0.75 meV larger than that of BA’ stacking.
The above analysis shows that the AB, BA, and BA’ are stable for bilayer GdI2.
To understand the effect of stacking on magnetism, we calculate the energy
difference between the interlayer AFM and FM states, as shown in Figures 3(b)
and 3(e) for the kinds of AA and AA’ stackings, respectively. Seeing is
believing, the magnitude of interlayer magnetic interaction for the kind of AA
stacking is stronger than that of the kind of AA’ stacking. Firstly, for the
kind of AA stacking, we find that the AA stacking strongly prefers AFM, and
the corresponding interlayer exchange energy is about -1.794 meV (see Figure
S5). On the contrary, the interlayer exchange energy of AB and BA stackings is
1.046 meV (see Figure S5), it indicates that the magnetic phase transition
occurs from AA stacking sliding to AB or BA stacking. Furthermore, for the
other (the kind of AA’ stacking), the vicinity of AB’ stacking prefers to FM
state, while the other regions favor AFM coupling. In a word, it indicates a
strong coupling between the magnetism and stacking order.

Are there other intriguing stacking orders? In order to answer the question,
we calculated the interlayer exchange energy of the kinds of AA and AA’
stackings for the entire 2D space of lateral shifts. We chose a $6\times 6$
grid to employ the lateral shifts. Figure 3(c, f) shows the energy difference
between interlayer AFM and FM states for the full 2D space of lateral shifts.
One can see that the magnetic ground state switches between AFM and FM states
as the interlayer stacking order changes. Moreover, the D3h symmetric AA
stacking prefers the AFM state. However, the magnetic interlayer interaction
become FM coupling for AB and BA stackings, when the D3h symmetric broken the
xy mirror symmetry (C3v). It is well known that the broken xy mirror symmetry
for D3h point group induces the out-of-plane electric polarization P 10. As
shown in Figure S6, the electrostatic potential difference with the positive
(negative) discontinuity $\Delta\Phi$ = 22.86 (-22.86) meV is generated
between the vacuum levels of upper and lower layers, which clearly proves that
the spontaneous P produced along +z (-z) axis in the AB (BA) stacking bilayer.
The electric polarization value of AB stacking bilayer GdI2 is up to
3.68$\times$10-12 C/m. The AB and BA stacking bilayers can realize reciprocal
transformation by interlayer sliding, therefore, we came to the conclusion
that the sliding ferroelectric exists in bilayer GdI2, which can also be
confirmed from the general theory of bilayer sliding ferroelectricity 29. It
indicates that the magnetic ground state can effectively be tuned by the
sliding ferroelectric. Note that the kind of AA’ stacking inexist spontaneous
polarization due to the 180∘ rotation between two layers induced the out-of-
plane polarization to cancel (see Figure S7).

To understand the microscopic mechanism of magnetic coupling, we consider the
simple spin Hamiltonian

$H=E_{0}+\sum_{i,j}J_{\|}S_{i}\cdot S_{j}+\sum_{i,k}J_{1\bot}S_{i}\cdot
S_{k}+\sum_{i,l}J_{2\bot}S_{i}\cdot S_{l},$ (5)

where $E_{0}$ is the ground state energy independent of the spin
configurations. _S_ i, _S_ j, _S_ k, and _S_ l denote the magnetic moments at
sites i, j, k, and l, respectively. _J_ ∥, _J_ 1⊥, and _J_ 2⊥ represent the
intralayer, nearest-neighbor (NN) interlayer, and second-neighbor (second-NN)
interlayer Gd-Gd exchange interactions, respectively (see Figure 4). The
details for the calculation of Heisenberg exchange parameters based on the
above spin Hamiltonian and total energy calculations employing DFT are shown
in the Supporting Information and Table I. The intralayer exchange is strongly
FM coupling ($\thicksim$ -8.25 meV), while the interlayer exchange interaction
depends on the stacking order.

Figure 4: (a) Schematic of the orbital dependent interlayer exchange
interactions. The hopping from d${}_{z^{2}}$-d${}_{z^{2}}$ is allowed in the
AFM exchange, whereas it is prohibited in the FM exchange. (b-e) The
interlayer Gd nearest neighbor and second neighbor corresponds to J1⊥ (red)
and J2⊥ (blue) for AA (b), AB (c), AB’ (d), and BA’ (e) stacking,
respectively. (f-i) The schematics of super-superexchange for AA (f), AB (g),
AB’ (h), and BA’ (i) stacking.

Since the Gd atoms in GdI2 are in a 4f75d1 electronic configuration (see
Figure S2), according to the Pauli exclusion principle and Hund’s rule, the
electron configuration will half fill the f and d${}_{z^{2}}$. In the trigonal
prismatic crystal field, the schematic illustrations of Gd-5d and Gd-4f
orbital splitting is shown in Figure S2(b). According to the energy minimizing
principle of hopping, the magnetic exchange between the two layers is
determined by the d orbital hopping. As shown in Figure 4(a), it exhibits the
schematic of the FM and AFM exchange interactions between Gd atoms in the
interlayer. The hopping between d${}_{z^{2}}$ and d${}_{z^{2}}$ is prohibited
for FM spin configuration (red), while it is allowed for AFM spin
configuration (blue). Hence, d${}_{z^{2}}$-d${}_{z^{2}}$ hybridization results
in AFM. On the other hand, hopping of the from
d${}_{z^{2}}$-dxy/d${}_{x^{2}-y^{2}}$ and d${}_{z^{2}}$-dxz/dyz gives rise to
an exchange coupling that is overwhelmingly FM from the viewpoint of the local
Hund coupling 30, 31. The interlayer Gd-Gd exchange interactions of all
stacking are mediated via the hybridization between the I-p orbitals.
Therefore, the nature of interlayer exchange interaction is super-
superexchange. The magnetism of stacking order dependence derives from the
competition between the orbital hybridizations of different interlayer.

In the following, we take four typical stacking as examples for detailed
analysis. For the AA stacking, Figure 4(b) shows the interlayer NN exchange
interaction of Gd-Gd J1⊥ (red) and second-NN exchange interaction of Gd-Gd J2⊥
(blue) for AA stacking. As shown in Figure 4(f), J1⊥ is predominated via
virtual excitations between Gd half-filled dz2 orbitals, causing an AFM
coupling 32. By contrast, J2⊥ is dominated by virtual excitations between Gd
half-filled d${}_{z^{2}}$ and the empty dxy/d${}_{x^{2}-y^{2}}$ orbitals,
leading to the FM coupling. As listed in Table I, although AA stacking has
only one J1⊥ bond and six J2⊥ bonds per unit cell, J1⊥ is much larger than
J2⊥. Therefore, the NN interlayer AFM coupling dominates the second-NN
interlayer FM coupling, resulting in the AFM state for AA stacking GdI2. For
the AB stacking, the sliding one layer t1//[$\frac{2}{3}$, $\frac{1}{3}$, 0]
with relative to the other breaks the interlayer hybridization between I-p
electrons and produces new ones. As shown in Figure 4(c), the NN ligancy
number is increased compared with the AA stacking the second-NN ligancy number
is decreased (see Table I). Simultaneously, the virtual excitations of J1⊥
become between Gd half-filled d${}_{z^{2}}$ and the empty
dxy/d${}_{x^{2}-y^{2}}$ orbitals, indicating the FM coupling. While J2⊥ is
dominated by virtual excitations between Gd half-filled d${}_{z^{2}}$ and the
empty dxz/dyz orbitals, contributing to also the FM coupling. As a result, the
magnetic ground state becomes FM state. Although the coordination number of
AB’ (BA’) stacking is consistent with AA (AB) stacking, they exhibit opposite
magnetic ground states. It originates from the symmetry difference [see Figure
S4, Figure 4(b-i), and Table I].

Stacking orders not only induces the changes of magnetic ground state and
ferroelectric polarization, the valley polarization is also transformed due to
the different symmetries. As shown in Figure 5(a) and Figure S8(a-d), the band
structures of AA stacking bilayer GdI2 are investigated based on the AFM
(M$\uparrow\downarrow$) magnetic ground state. In the absence of SOC, spin up
and spin down bands are degenerate [see Figure S8(a)]. When the SOC is
switched on, as shown in Figure 5(a), the degeneracy of spin up and spin down
bands disappear. In detail, we focus on the valence band. The spin up band
shifts above the spin down band in K valley, while the spin up band shifts
below the spin down band in K’ valley. More interestingly, the spin up and
spin down bands of K and K’ valleys are degenerate in energy. From the
orbital-resolved band structure [see Figure S8(b)], the VBM bands are also
dominated by Gd dxy and d${}_{x^{2}-y^{2}}$ orbitals, it is consistent with
monolayer GdI2. Moreover, as shown in Figure 5(b), the Berry curvatures of K
and K’ valleys have equal magnitudes and opposite signs (26.2 ${\rm\AA}^{2}$
for K point and -26.2 ${\rm\AA}^{2}$ for K’ point).

Figure 5:  Spin-resolved band structures and Berry curvatures of AA (a, b), AB
(c, d), and BA (e, f) stacking bilayer GdI2 with considering SOC effect.

When the AA stacking slides to AB stacking, the degeneracy of spin up and spin
down bands disappear in the absence of SOC [see Figure S8(e)]. It originates
from the magnetic ground state from AFM state transition to FM state. When the
SOC is included, as shown in Figure 5(c), it produces a valley polarization of
158 meV. It indicates that the interlayer sliding induces the out-of-plane
ferroelectric polarization, simultaneously, leading to the magnetic phase
transition and valley polarization. Due to the valley polarization generating,
as shown in Figure 5(d), the Berry curvature magnitudes of K and K’ valleys
are no longer equal, and the signs are still opposite (25.0 ${\rm\AA}^{2}$ for
K point and -28.8 ${\rm\AA}^{2}$ for K’ point). When further slide to BA
stacking, the valley polarization is reversed by the ferroelectric switching,
as shown in Figure 5(e). Meanwhile, the Berry curvatures of K and K’ valleys
change to 28.8 ${\rm\AA}^{2}$ and -25.0 ${\rm\AA}^{2}$ [see Figure 5(f)]. In
addition, the other stacking also produces intriguing coupling of magnetic,
ferroelectric, and ferrovalley orders (see Figure S9).

In conclusion, we propose a mechanism to realize coexisting magnetism,
ferroelectric, and ferrovalley multiferroic in 2D materials. The mechanism is
proved in bilayer GdI2. By transforming the interlayer stacking order, one can
regulate the magnetic ground state, ferroelectric polarization, and valley
polarization for bilayer GdI2. The magnetic ground state of AA stacking is AFM
state with the degeneracy at the VBM of K and K’ valleys. As the AA stacking
slides to AB (BA) stacking, its spontaneous valley polarization occurs due to
the magnetic ground state transforming from AFM to FM state. The AB stacking
changes to the BA stacking state by interlayer sliding, and the ferroelectric
polarization and valley polarization are switched. Moreover, we reveal the
microscopic mechanism of magnetic phase transition by spin Hamiltonian and
electron hopping between layers. Our work offers not only a novel 2D
multiferroic material but also an efficient means to tune magnetic,
ferroelectric, and ferrovalley properties.

###  0.1 AUTHOR INFORMATION

Corresponding Authors
<EMAIL_ADDRESS>
Notes
The authors declare no competing financial interest.

This work is supported by the National Natural Science Foundation of China
(Grant No. 12004295). P. Li thanks China’s Postdoctoral Science Foundation
funded project (Grant No. 2022M722547), the Fundamental Research Funds for the
Central Universities (xxj03202205), and the Open Project of State Key
Laboratory of Surface Physics (No. KF2022$\\_$09).

Details of the calculation methods, spin Hamiltonian model, spin charge
density, orbital-resolved band structures, the effect of U value on valley
polarization, stacking structures of bilayer GdI2, the energy difference
between interlayer AFM and FM for bilayer, the plane averaged electrostatic
potential, and band structures for bilayer GdI2.

## References

  * 1 Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. _Science_ 2004, 306, 666.
  * 2 Cao, T.; Wang, G.; Han, W.; Ye, H.; Zhu, C.; Shi, J.; Niu, Q.; Tan, P.; Wang, E.; Liu, B.; Feng, J.; Valley-selective circular dichroism of monolayer molybdenum disulphide. _Nat. Commun._ 2012, 3, 887.
  * 3 Gong, C.; Li, L.; Li, Z.; Ji, H.; Stern, A.; Xia, Y.; Cao, T.; Bao, W.; Wang, C.; Wang, Y.; Qiu, Z. Q.; Louie, S. G.; Xia, J.; Zhang, X.; Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals. _Nature_ 2017, 546, 265.
  * 4 Huang, B.; Clark, G.; Navarro-Moratalla, E.; Klein, D. R.; Cheng, R.; Seyler, K. L.; Zhong, D.; Schmidgll, E.; McGuire, M. A.; Cobden, D. H.; Yao, W.; Xiao, D.; Jarillo-Herrero, P.; Xu, X.; Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit. _Nature_ 2017, 546, 270.
  * 5 Li, W.; Qian, X.; Li, J.; Phase transitions in 2D materials. _Nat. Rev. Mater._ 2012, 6, 829.
  * 6 Li, P.; Zhou, X. S.; Guo, Z. X.; Intriguing magnetoelectric effect in two-dimensional ferromagnetic/perovskite oxide ferroelectric heterostructure. _npj Comput. Mater._ 2022, 8, 20.
  * 7 Fei, Z.; Zhao, W.; Palomaki, T. A.; Sun, B.; Miller, M. K.; Zhao, Z.; Yan, J.; Xu, X.; Cobden, D. H.; Ferroelectric switching of a two-dimensional metal. _Nature_ 2018, 560, 336.
  * 8 Xiao, J.; Wang, Y.; Wang, H.; Pemmaraju, C. D.; Wang, S.; Muscher, P.; Sie, E. J.; Nyby, C. M.; Devereaux, T. P.; Qian, X.; Zhang, X.; Lindenberg, A. M.; Berry curvature memory through electrically driven stacking transitions. _Nat. Phys._ 2020, 16, 1028.
  * 9 Wang, X.; Yasuda, K.; Zhang, Y.; Liu, S.; Watanabe, K.; Taniguchi, T.; Hone, J., Fu, L.; Jarillo-Herrero, P.; Interfacial ferroelectricity in rhombohedral-stacked bilayer transition metal dichalcogenides. _Nat. Nanotechnol._ 2022, 17, 367.
  * 10 Li, L.; Wu, M.; Binary compound bilayer and multilayer with vertical polarizations: two-dimensional ferroelectrics, multiferroics, and nanogenerators. _ACS Nano_ 2017, 11, 6382.
  * 11 Yasuda, K.; Wang, X.; Watanabe, K.; Taniguchi, T.; Jarillo-Herrero, P.; Stacking-engineered ferroelectricity in bilayer boron nitride. _Science_ 2021, 372, 1458.
  * 12 Sivadas, N.; Okamoto, S.; Xu, X.; Fennie, C. J.; Xiao, X.; Stacking-Dependent Magnetism in Bilayer CrI3. _Nano Lett._ 2018, 18, 7658.
  * 13 Chen, W.; Sun, Z.; Wang, Z.; Gu, L.; Xu, X.; Wu, S.; Gao, C.; Direct observation of van der Waals stacking-dependent interlayer magnetism. _Sicence_ 2019, 366, 983.
  * 14 Li, X.; Xu, X.; Zhou, H.; Jia, H.; Wang, E.; Fu, H.; Sun, J. T.; Meng, S.; Tunable topological states in stacked chern insulator bilayers. _Nano Lett._ 2023, 23, 2839.
  * 15 Peng, R.; Zhang, T.; He, Z.; Wu, Q.; Dai, Y.; Huang, B.; Ma, Y.; Intrinsic layer-polarized anomalous Hall effect in bilayer MnBi2Te4. _Phys. Rev. B_ 2023, 107, 085411.
  * 16 Feng, Y.; Dai, Y.; Huang, B.; Kou, L.; Ma, Y.; Layer Hall Effect in Multiferroic Two-Dimensional Materials. _Nano Lett._ 2023, 23, 5367.
  * 17 Zhang, T.; Xu, X.; Guo, J.; Dai, Y.; Ma, Y.; Layer-Polarized Anomalous Hall Effects from Inversion-Symmetric Single-Layer Lattices. _Nano Lett._ 2024, 24, 1009.
  * 18 Wu, Y.; Tong, J.; Deng, L.; Luo, F.; Tian, F.; Qin, G.; Zhang, X.; Coexisting Ferroelectric and Ferrovalley Polarizations in Bilayer Stacked Magnetic Semiconductors. _Nano Lett._ 2023, 23, 6226.
  * 19 Li, Y. Q.; Zhang, X.; Shang, X.; He, Q. W.; Tang, D. S.; Wang, X. C.; Duan, C. G.; Magnetic and Ferroelectric Manipulation of Valley Physics in Janus Piezoelectric Materials. _Nano Lett._ 2023, 23, 10013.
  * 20 Gong, C.; Kim, E. M.; Wang, Y.; Lee, G.; Zhang, X.; Multiferroicity in atomic van der Waals heterostructures. _Nat. Commun._ 2019, 10, 2657.
  * 21 Sun, W.; Wang, W.; Li, H.; Zhang, G.; Chen, D.; Wang, J.; Cheng, Z.; Controlling bimerons as skyrmion analogues by ferroelectric polarization in 2D van der Waals multiferroic heterostructures. _Nat. Commun._ 2020, 11, 5930.
  * 22 Li, P.; Wu, C.; Peng, C.; Yang, M.; Xun, W.; Multifield tunable valley splitting in two-dimensional MXene Cr2COOH. _Phys. Rev. B_ 2023, 108, 195424.
  * 23 Tong, W. Y.; Gong, S. J.; Wan, X.; Duan, C. G.; Concepts of ferrovalley material and anomalous valley Hall effect. _Nat. Commun._ 2016, 7, 13612.
  * 24 Li, P.; Liu, B.; Chen, S.; Zhang, W. X.; Guo, Z. X.; Progress on two-dimensional ferrovalley materials. _Chin. Phys. B_ 2024, 33, 017505.
  * 25 Li, P.; Yang, X.; Jiang, Q. S.; Wu, Y. Z.; Xun, W.; Built-in electric field and strain tunable valley-related multiple topological phase transitions in VSiXN4 (X = C, Si, Ge, Sn, Pb) monolayers. _Phys. Rev. Mater._ 2023, 7, 064002.
  * 26 Zhang, A.; Gong, Z.; Zhu, Z.; Pan, A.; Chen, M.; Effects of the substrate-surface reconstruction and orientation on the spin valley polarization in MoTe2/EuO _Phys. Rev. B_ 2020, 102, 155413.
  * 27 Wang, K.; Li, Y.; Mei, H.; Li, P.; Guo, Z. X.; Quantum anomalous Hall and valley quantum anomalous Hall effects in two-dimensional d0 orbital XY monolayers. _Phys. Rev. Mater._ 2022, 6, 044202.
  * 28 Zhang, C.; Guo, P.; Zhou, J.; Tailoring bulk photovoltaic effects in magnetic sliding ferroelectric materials. _Nano Lett._ 2022, 22, 9297.
  * 29 Ji, J.; Yu, G.; Xu, C.; Xiang, H. J.; General Theory for Bilayer Stacking Ferroelectricity. _Phys. Rev. Lett._ 2023, 130, 146801.
  * 30 Goodenough, J. B.; Theory of the Role of Covalence in the Perovskite-Type Manganites [La, M(II)]MnO3. _Phys. Rev._ 1955, 100, 564.
  * 31 Kanamori, J.; Superexchange interaction and symmetry properties of electron orbitals. _J. Phys. Chem. Solids_ 1959, 10, 87.
  * 32 Anderson, P. W.; Antiferromagnetism. Theory of superexchange interaction. _Phys. Rev._ 1950, 79, 350.

## Graphical TOC Entry


# The Astrophysical Distance Scale III:
Distance to the Local Group Galaxy WLM using Multi-Wavelength
Observations of the Tip of the Red Giant Branch, Cepheids, and JAGB Stars

Abigail J. Lee  Department of Astronomy & Astrophysics, University of Chicago,
5640 South Ellis Avenue, Chicago, IL 60637<EMAIL_ADDRESS>Wendy L.
Freedman  Department of Astronomy & Astrophysics, University of Chicago, 5640
South Ellis Avenue, Chicago, IL 60637  Barry F. Madore  Department of
Astronomy & Astrophysics, University of Chicago, 5640 South Ellis Avenue,
Chicago, IL 60637  Observatories of the Carnegie Institution for Science 813
Santa Barbara St., Pasadena, CA 91101  Kayla A. Owens  Department of Astronomy
& Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL
60637  Andrew J. Monson  Department of Astronomy & Astrophysics, Pennsylvania
State University, 525 Davey Lab, University Park, PA 16802  Taylor J. Hoyt
Department of Astronomy & Astrophysics, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL 60637

###### Abstract

The local determination of the Hubble Constant sits at a crossroad. Current
estimates of the local expansion rate of the Universe differ by about
1.7-$\sigma$, derived from the Cepheid and TRGB based calibrations, applied to
type Ia supernovae. To help elucidate possible sources of systematic error
causing the tension, we show in this study the recently developed distance
indicator, the J-region Asymptotic Giant Branch (JAGB) method (Madore &
Freedman, 2020), can serve as an independent cross-check and comparison with
other local distance indicators. Furthermore, we make the case that the JAGB
method has substantial potential as an independent, precise and accurate
calibrator of type Ia supernovae for the determination of $H_{0}$. Using the
Local Group galaxy, WLM we present distance comparisons between the JAGB
method, a TRGB measurement at near-infrared ($JHK)$ wavelengths, a TRGB
measurement in the optical $I$ band, and a multi-wavelength Cepheid period-
luminosity relation determination. We find:

$\mu_{0}~{}(JAGB)=24.97\pm 0.02\text{ (stat)}\pm 0.04\text{ (sys) mag}$
$\mu_{0}~{}(TRGB_{NIR})=24.98\pm 0.04\text{ (stat)}\pm 0.07\text{ (sys) mag}$
$\mu_{0}~{}(TRGB_{F814W})=24.93\pm 0.02\text{ (stat)}\pm 0.06\text{ (sys)
mag}$ $\mu_{0}~{}(Cepheids)=24.98\pm 0.03\text{ (stat)}\pm 0.04\text{ (sys)
mag}$

All four methods are in good agreement, confirming the local self-consistency
of the four distance scales at the 3% level, and adding confidence that the
JAGB method is as accurate and as precise a distance indicator as either of
the other three astrophysically-based methods.

distance scale – stars: Population II – galaxies: individual (WLM) – galaxies:
stellar content – stars: AGB and post-AGB – stars: carbon – stars: variables:
Cepheids – cosmology: observations

††software: daophot (Stetson, 1987), allframe (Stetson, 1994), daogrow
(Stetson, 1990), TinyTim (Krist et al., 2011), topcat (Taylor, 2005), Astropy
(Astropy Collaboration et al., 2013, 2018), NumPy (van der Walt et al., 2011),
Matplotlib (Hunter, 2007), scipy (Virtanen et al., 2020)

##  1 Introduction

Measuring both precise and accurate distances to galaxies persists as a
challenging problem in astronomy. Attempts to discover, develop and refine new
distance indicator techniques in the past two decades have only increased as
the need for greater accuracy has grown. One goal of these endeavors, a robust
calibration of type Ia supernovae, would provide a more secure determination
of $H_{0}$ in the nearby universe, and a better understanding of the magnitude
and significance of the current tension between local (late-time) values of
$H_{0}$ and those derived from (early-time) cosmological modeling of cosmic
microwave background (CMB) observations. Resolving the Hubble tension will
require having systematic effects constrained at the 1-2% percent-level (or
better) in their individual accuracies. The two current front-runners for
high-precision measurements of local distances are the Tip of the Red Giant
Branch (TRGB) method used by the Carnegie-Chicago Hubble Program (CCHP)
(Freedman et al., 2019, 2020), and Cepheid Period-Luminosity based method (the
Leavitt Law) used by the Hubble Space Telescope Key Project for the past 30
years (Freedman et al., 2001, 2012) and more recently by the SHoES group
(Riess et al., 2016; Riess et al., 2019). Currently however, these two local
(TRGB vs Cepheid) measurements of $H_{0}$ differ by almost 1.7-$\sigma$, and
this tension continues. Is there an independent distance indicator, equal in
precision and accuracy to the Cepheids and TRGB methods, that can serve as a
cross-check to shed light on differences in systematics between the two more
classical methods? We are now exploring that possibility.

The aim of this paper is to further develop an independent distance
calibrator, the near-infrared J-region Asymptotic Giant Branch (JAGB) method,
recently (re)introduced into the distance scale arena (Madore & Freedman,
2020; Freedman & Madore, 2020; Ripoche et al., 2020). JAGB stars are
Thermally-Pulsating Asymptotic Giant Branch (TP-AGB) stars with considerable
amounts of carbon in their atmospheres (i.e., carbon-to-oxygen ratios of
C/O$>1$) (Habing & Olofsson, 2003). JAGB stars also undergo ‘dredge-up’
events, where with each thermal pulse the convective envelope penetrates into
deeper and deeper layers of the star’s interior. During the third dredge-up
phase (and beyond), the convective envelope of the star penetrates deep enough
that it encounters carbon, produced by the He-burning shell, and brings that
enriched material to the surface (see Habing & Olofsson (2003) for a general
review and Marigo & Girardi (2007); Marigo et al. (2017) for detailed
descriptions of TP-AGB evolution).

For the younger, more massive AGB stars ($1.6\times 10^{8}$ years), the
temperature at the bottom of the convective envelope becomes so hot that
during the third dredge-up phase, the recently-formed carbon is quickly
converted into nitrogen before it can be convected to the surface. This places
an upper limit on the luminosity of a JAGB star (Iben, 1973; Sackmann et al.,
1974). For older, less massive AGB stars, the third dredge-up phase (and
beyond) becomes ineffective in transporting carbon to the surface. This
imposes a lower limit on the luminosities of stars entering the JAGB region
(Iben & Renzini, 1983). These theoretical limits are reflected in the observed
properties JAGB stars, which are well-defined in their near-infrared (NIR)
colors and luminosities. JAGB stars are redder and quite distinct from the
bluer oxygen-rich AGB stars, especially at near and mid-infrared wavelengths:
a direct result of the large amounts of carbon in their atmospheres (Marigo &
Girardi, 2007). Thus, JAGB stars have been observed to have well-defined
limits on their absolute magnitudes and color, making them distinct and easily
identified in color-magnitude space. For a more comprehensive review of the
dredge-up theory for carbon stars, see Habing & Olofsson (2003).

A review of carbon stars by Battinelli & Demers (2005) demonstrated that the
mean luminosity of the distribution of JAGB stars is relatively constant from
galaxy to galaxy, and preliminary tests carried out by Freedman & Madore
(2020) have also demonstrated the JAGB absolute magnitude has little, if any,
dependence on metallicity or age. Empirically, the mean of the JAGB luminosity
distribution has all of the necessary attributes of an excellent distance
indicator.

Both the TRGB and JAGB methods provide a standard candle. To compare the two
methods, Freedman & Madore (2020) measured and compiled distances to 14
galaxies, and found the TRGB and JAGB distances to be in good agreement.
However, this study was based on previously-published data, not optimized for
the detection of JAGB stars. No studies to date have compared distance moduli
using both the optical and near-infrared TRGB and the JAGB method derived from
the same dataset, using self-consistent processing, calibration and analysis
procedures. Accordingly, in this study we undertake just such an in-depth
study of the galaxy Wolf-Lundmark-Melotte (WLM).

WLM is a dwarf irregular galaxy located in the outer edges of the Local Group
(Wolf, 1909; Melotte, 1926). We have measured its distance modulus using the
TRGB in both the near-infrared $JHK$ bands and optical $I$ band, a multi-
wavelength Cepheid period-luminosity relation determination, and the JAGB
method in the $J$ band. We have obtained high-precision ground-based data,
which we have calibrated using archival HST observations for the optical and
Two-Micron All-Sky Survey (2MASS) data (Skrutskie et al., 2006) for the near
infrared.

The first goal of this paper is to test the degree to which the JAGB method
can be used as a cross check and comparison with other local distance
indicators, such as Cepheids and the TRGB, in order to reveal potentially
obfuscated systematics. The second goal of this paper is to demonstrate that
with future development and testing, using both ground-based and space-based
data, the JAGB method is a distance indicator that is comparable to the TRGB
and Cepheids in accuracy and precision, and can eventually become the basis
for an independent determination of the Hubble Constant.

Our paper outline is as follows. In Section 2, we describe the imaging
datasets and the data-reduction methods used in this study. In Section 3, we
give a brief overview of the multi-wavelength application of the TRGB method,
and then present our measurements of the distance modulus to WLM in the near-
infrared and independently in the optical. In Section 4, we explain the newly
developed JAGB method, and then give our measured distance modulus. In Section
5, we calculate a Cepheid-based distance modulus, as well as compare our TRGB
distance moduli and JAGB distance modulus values with previous measurements
for WLM. Finally, we state our conclusions and look to the future with the
JAGB method in Section 6. We also provide supplementary information on our
photometric procedure in Appendix A.

##  2 Data and Photometry

Figure 1: Images of the Local Group Galaxy WLM data. In Figure (a), $I$-band
greyscale images for the eight IMACS chips are shown. The yellow rectangles
show the locations of the two HST ACS/WFC chips (F814W) used to calibrate our
$I$-band ground-based photometry. Figure (b) shows the positions of all the
queried 2MASS stars (red) overlaid onto the FourStar data (black) taken
September 7, 2011. To give our IMACS and FourStar data right ascension /
declination ($\alpha/\delta$) coordinates, we matched the data with HST/2MASS
data that did have $\alpha/\delta$ positions.

We now describe the five imaging datasets used in this study. A summary of the
data can be found in Table 1. We first describe our process for obtaining
final $VI$ photometry in Section 2.1, and then describe our process for
obtaining $JHK$ photometry in Section 2.2.

Table 1: A Summary of Observations used in this work Dates | Instrument | Filter(s) | $\alpha$ | $\delta$ | Exp. (sec) | Target
---|---|---|---|---|---|---
2019-11-16 | IMACS | $VI$ | 00:01:58.0 | -15:27:39.4 | 1,800 sec | $TRGB_{F814W}$
2015-07-17,18,19 | ACS/WFC | F814W | 00:01:57.3 | -15:31:22.7 | 34,320 sec |  $I$ calibration
2011-09-07 | FourStar | $JHK$ | 00:01:57:0 | -15:27:36.0 |  $625$ sec |  $TRGB_{JHK}$, $JAGB$
2011-10-05 | FourStar | $JHK$ | 00:01:57:0 | -15:27:36.0 |  $1,083$ sec | $JAGB$
2019-11-13 | FourStar | $JHK$ | 00:01:57:0 | -15:27:36.0 |  $39$ sec |  $JHK$ Calibration

###  2.1 VI Photometry

####  2.1.1 Magellan-Baade Telescope: IMACS

Ground-based $VI$ data of WLM were obtained at the _Inamori-Magellan Areal
Camera and Spectrograph_ (Dressler et al., 2011, IMACS) on the $6.5$m
Magellan-Baade telescope using the f/4 imaging mode at Las Campanas
Observatory on November 16, 2019. The field has a dimension of
$15.4\arcmin\times 15.4\arcmin$ at a scale of 0.111 arcsec/pixel. The raw CCD
data were then reduced using standard procedures: the overscan region was
first subtracted by column from the raw science image; then the data were
flat-fielded.

####  2.1.2 HST Archival Data

We utilized $HST$ archival data in order to calibrate the $I$-band ground-
based photometry to the Vegamag magnitude system. We found one HST/ACS dataset
of WLM in the HST archives that overlapped with the IMACS chips: _Completing
the Census of Isolated Dwarf Galaxy Star Formation Histories_ (Weisz, 2014,
PID: GO13768, PI: Weisz). The goal of this dataset was to measure the star
formation histories of isolated dwarf galaxies by obtaining deep optical data.
This dataset contained F814W data (the $HST$ filter corresponding to the Kron-
Cousins $I$ band) which we used to tie our $I$-band ground-based photometry
onto the HST flight magnitude system. We note that color was used only to
obtain our $VI$ color magnitude and visually check that the tip detection was
reliable; the $V$-band are not used in the actual measurement of the TRGB in
the $I$-band luminosity function.

The $HST$ dataset contained 13 pointings, each 2,640 seconds, for a total
exposure time of 34,320 seconds. The location of this pointing relative to the
8 IMACS chips is shown in Figure 1. The STScI processed and charge-transfer-
efficiency-corrected and flat-fielded individual frames (the .flc data type)
were pulled from the Mikulski Archive for Space Telescopes (MAST). We
corrected for pixel-to-pixel geometric distortions caused by the optical
design of the ACS camera by multiplying each image by its respective Pixel
Area Map, which was pulled from the stsci.skypac
module.111https://www.stsci.edu/hst/instrumentation/acs/data-analysis/pixel-
area-maps

####  2.1.3 $VI$ Instrumental Photometry

To obtain instrumental magnitudes from our imaging datasets, we followed the
standard photometry procedure for obtaining PSF-photometry outlined in the
daophot-II User Manual (Stetson, 2000). The procedures were identical for both
the HST archival data and the IMACS data with the exception that for the IMACS
data, empirical point-spread functions were manually created, and for the HST
data, theoretical model-based TinyTim point-spread functions (Krist et al.,
2011) were generated. TinyTim PSFs were used here to be consistent with the
$I$-band absolute magnitude TRGB calibration undertaken by the CCHP (Freedman
et al., 2019), which also used TinyTim PSFs.

The basic photometry procedure is as follows: for a given band, photometry was
first performed on each frame using the daophot suite, resulting in
instrumental magnitudes for each frame. Then, for a given chip, the images
were aligned using daomatch/daomaster, and then a medianed image was created
with montage2. Subsequently, the medianed image was photometered to create a
master source list, and then all of the individual frames were simultaneously
photometered using allframe and the master source list. Finally, all of the
final photometry files were matched using daomaster again to give final
instrumental magnitudes for a given chip. This process was repeated for all
chips and bands.

####  2.1.4 Calibrating HST Archival Data

In order to calibrate the instrumental HST magnitudes onto the Vegamag
magnitude system, we followed the procedure outlined in Sirianni et al.
(2005). First, the flight magnitude zero point was obtained from the STScI on-
line calculator222https://acszeropoints.stsci.edu/. For our HST observations
taken July 17, 2015, the ACS zeropoint reported was $ZP_{F814W}=25.517$ mag.
Second, the infinite aperture correction, the correction from a
$0.5\arcsec$-to-infinity aperture, was obtained from Bohlin (2016), determined
to be 0.098 mag. We adopted a $\pm 0.02$ mag error for both the zeropoint and
infinite aperture corrected consistent with the CCHP program (Beaton et al.,
2019). Third, we needed the correction from the PSF magnitudes to the
$0.5\arcsec$ aperture. This correction was obtained by first identifying
bright and isolated stars in the medianed image, and then re-locating them all
in 26 frames. We computed growth curves for each star using daogrow (Stetson,
1990), and eliminated the few stars where the growth curve showed evidence
that the star’s flux was being contaminated from neighboring stars. Now left
with 30 stars for a given frame, the aperture magnitude at $0.5\arcsec$ minus
the PSF magnitude was calculated for all the stars. After eliminating outliers
more than 2-$\sigma$ from the median value, we averaged all the stars’
individual aperture corrections to give the final aperture correction for a
given frame. Aperture corrections were then applied on a frame-by-frame basis.
The average aperture correction across all frames was $-0.107\pm 0.009$ mag.

####  2.1.5 Calibration of $I$-band Photometry onto HST Flight Magnitude
System

The last step of the calibration process was to bring the IMACS instrumental
magnitudes onto the Vegamag magnitude system. This was done by finding stars
in both the calibrated HST data and the IMACS instrumental data and computing
the mean magnitude offset between them, and then applying that offset to the
IMACS data. In order to prioritize the accuracy of the TRGB stars’ magnitudes,
stars were chosen based on their $I$-band magnitude near the TRGB. We chose a
range of magnitudes based off of the reported WLM TRGB value of
$m_{F814W}=20.91$ mag from Albers et al. (2019), a study that used the same
HST data as this study. Thus, we adopted a range of 20.80 and 21.30 mag for
our calibration stars. In total, 23 overlapping TRGB stars were found in both
the HST ACS/WFC footprint and the IMACS data. We show all of the stars’
calculated magnitude offsets in Figure 2. The mean magnitude offset between
the HST calibrated magnitudes and IMACS instrumental magnitudes was determined
to be $8.350\pm 0.003$ mag. This offset was applied to the instrumental
magnitudes from IMACS.

Figure 2: Figure (a) shows the magnitude offsets between our IMACS
instrumental photometry and HST calibrated photometry. Previous studied using
the same HST data showed the $TRGB_{F814W}$ to be around 20.91 mag (Albers et
al., 2019), so we restricted our search for stars in both the IMACS $I$-band
data and the HST archival data to be between $20.80<m_{F814W}<21.30$ mag in
order to most accurately measure the magnitudes of our targeted TRGB stars.
The average magnitude offset was determined to be $<F814W-I>=8.350\pm.003$
mag. Figure (b) shows the magnitude offsets between our FourStar photometry
and 2MASS standards.

The photometry was also cleaned of extended objects through the use of
photometric quality cuts on the following daophot parameters: the photometric
uncertainty $\sigma$, the sharp-parameter ($sharp$), and the chi-parameter
$\chi$. Cuts were made based on a constant+exponential function as a function
of magnitude. The specific parameters for the exponential functions can be
found in the Appendix A.

After making photometric quality cuts, we plotted a calibrated $I$ vs. $(V-I)$
color magnitude diagram. In Figure 3, we show this final color magnitude
diagram, $I$-band luminosity function, and TRGB edge detection. The plotted
stars are all in the halo of WLM, eliminating crowding effects and AGB
contamination from the disk. In Section 3, we describe our methods for the
TRGB detection process for the optical $I$ band, and our subsequent tip
detections in the near-infrared.

Figure 3: Color Magnitude Diagram for the optical $I$ band. The TRGB detection
is clear and distinct. Error bars show uncertainties at 1 mag intervals. The
middle panel shows the $I$-band GLOESS-smoothed luminosity function, and the
right panel shows the edge response function.

###  2.2 JHK Photometry

####  2.2.1 Magellan-Baade Telescope: FourStar

Ground-based $JHK$ data were taken at the $6.5$m Magellan-Baade telescope at
Las Campanas Observatory with the FourStar Infrared Camera (Persson et al.,
2013). The FourStar imager has a field of view of $10.8\arcmin\times
10.8\arcmin$ and a resolution of 0.159 arcsec/pixel. We took three separate
sets of observations, detailed in Table 1. The first set of observations,
taken September 7, 2011, targeted the near-infrared TRGB and JAGB stars. The
second set of observations, taken on October 5, 2011, also targeted JAGB
stars. Finally, the third set of observations, used for calibration purposes,
was taken November 13, 2019.

####  2.2.2 JHK Photometry

To obtain PSF photometry on the $JHK$ data, daophot was run on each frame
using the standard procedure outlined in Section 2.1.3. The next step was then
to calibrate our instrumental $JHK$ photometry onto the 2MASS filter system,
adopted for the $TRGB_{JHK}$ calibrations by Hoyt et al. (2018). We utilized
the 2MASS single-epoch photometry in the $JHK$ filters.

The infrared $JHK$ data taken September 7, 2011 and October 5, 2011 were first
moved onto the same instrumental magnitude system as the shallower dataset
taken November 13, 2019 as the original datasets were too deep and thus
contained only one or two unsaturated stars in common with the 2MASS data.
These data were then calibrated using data from 2MASS by finding stars in
common with the 2MASS data and the FourStar dataset in a similar fashion as
that described in Section 2.1.5. See Figure 2 for a plot of the FourStar to
2MASS calibration. We applied photometric quality cuts in an identical fashion
to our $VI$ data to clean it of extended objects. The specific exponential
function forms used can be found in Appendix A.

In Section 3.2 and Section 4, we describe how we use our $JHK$ data for our
near-infrared TRGB detections and JAGB measurement, respectively.

###  2.3 Galactic Foreground Extinction and Internal Reddening Corrections

WLM is at high galactic latitude ($b=-73.6^{\circ}$) and is therefore expected
to have a low line-of-sight extinction. We determined the foreground
extinction correction using the online IRSA Galactic Dust Reddening and
Extinction tool333https:/irsa.ipac.caltech.edu/applications/DUST/, which
queries the Schlegel et al. (1998) full-sky Milky Way extinction maps and
gives reddening estimates based on a location and rescaling from Schlafly &
Finkbeiner (2011). The foreground extinction corrections for this study in the
$I,J,H,K$ bands respectively amount to $A_{I}=0.057,A_{J}=0.027,A_{H}=0.017$
and $A_{K}=0.012$ mag, assuming the Cardelli et al. (1989) reddening law,
$R_{v}=3.1$. As per Beaton et al. (2019), we choose to adopt a systematic
uncertainty on the extinction correction equal to half its value.

To measure possible internal reddening differences resulting from in situ dust
of WLM, we split our field into three non-overlapping subsections. For each
field, the luminosity function was first smoothed using a smoothing parameter
of $\sigma_{s}=0.10$, which is slightly larger than that which was used for
the entire field, $\sigma_{s}=0.075$ (discussed in Section 3.1) in order to
account for the smaller numbers of stars. We then calculated the $TRGB_{I}$
for each section. We found a variation of 0.03 mag between the three fields.
Because the variations between fields is so small, we conclude that there is a
negligible differential internal reddening effect, as this variation is most
likely just due to small number statistics (e.g., see Hatt et al., 2017).

##  3 Tip of the Red Giant Branch

The use of the TRGB as a distance indicator can be understood within the
context of stellar evolution for low-mass red giant stars (Salaris & Cassisi,
1997). The beginning of helium burning by the triple-$\alpha$ process in the
star’s core is marked by a sharp discontinuity in the red giant branch
luminosity function (at the location of the helium flash) (Lee et al., 1993).

In the $I$ band, the luminosity of the TRGB is remarkably insensitive to
color, and can used as a standard candle (without additional color
corrections) (Freedman et al., 2019; Jang et al., 2020). The TRGB distance
indicator method is also well understood theoretically (Salaris & Cassisi,
1997) and is simple and empirically-based. Recent studies have also shown the
TRGB can also be be used as a distance indicator in the near-infrared $JHK$
bands even though the TRGB becomes upward-sloping with increasing color
(Madore et al., 2018; Hoyt et al., 2018).

For all of our TRGB detections, we made sure to only include stars in the halo
of WLM in order to minimize contamination of the TRGB by brighter
intermediate-age AGB stars, in situ reddening from gas and dust in the disk,
and to avoid crowding/blending in our photometry.

We first constructed four TRGB color magnitude diagrams: $J$ vs. $(J-K)$, $K$
vs. $(J-K)$, $H$ vs. $(J-K)$, and $I$ vs. $(V-I)$. As mentioned above, in the
optical $I$ band, the TRGB has a nearly constant magnitude, with little to no
dependence on age or metallicity in the color range found for the RGB in WLM,
whereas in the near-infrared, the TRGB is known to be upward sloping with
increasing color. In Section 3.1, we describe our measurement of the $I$-band
TRGB, and in Section 3.2, we describe our subsequent measurements of the near-
infrared TRGB zero point.

###  3.1 Optical $I$-band TRGB

Table 2: Error Budget for WLM  Source of Uncertainty | Value [mag] |  $\sigma_{stat}$ [mag] |  $\sigma_{sys}$ [mag]
---|---|---|---
Edge Detection | $~{}...$ | 0.01 | 0.01
STScI ACS F814W ZP | 25.517 | $~{}...$ | 0.02
STScI ACS F814W Infinite Aperture Correction | 0.098 | $~{}...$ | 0.02
Differential Aperture Correction | -0.107 | 0.009 | $~{}...$
$<$ACS $-$ IMACS$>$ | 8.350 | 0.003 | $~{}...$
$F814W(TRGB)$ | 20.93 | 0.01 | 0.03
$A_{F814W}$ | 0.057 | $~{}...$ | 0.029
$M^{TRGB}_{F814W}$ (Freedman et al., 2020) | -4.054 | 0.022 | 0.039
$\mu_{0}(TRGB)_{F814W}$ | 24.93 | 0.02 | 0.06
$J(TRGB)$ | 19.85 | 0.036 | $0.017$
$A_{J}$ | 0.027 | $~{}...$ | 0.014
$M^{TRGB}_{J}$ (Hoyt et al., 2018) | -5.14 | 0.01 | 0.06
$\mu_{0}(TRGB)_{J}$ | 24.96 | 0.04 | 0.06
$H(TRGB)$ | 19.05 | 0.038 | 0.022
$A_{H}$ | 0.017 | $~{}...$ | 0.029
$M^{TRGB}_{H}$ (Hoyt et al., 2018) | -5.94 | 0.01 | 0.06
$\mu_{0}(TRGB)_{H}$ | 24.97 | 0.04 | 0.07
$K(TRGB)$ | 18.85 | 0.036 | 0.030
$A_{K}$ | 0.012 | $~{}...$ | 0.029
$M^{TRGB}_{K}$ (Hoyt et al., 2018) | -6.14 | 0.01 | 0.06
$\mu_{0}(TRGB)_{K}$ | 24.98 | 0.04 | 0.08
$J(JAGB)$ | 18.80 | 0.02 | 0.007
$A_{J}$ | 0.027 | $~{}...$ | 0.014
$M_{J}^{JAGB}$ (Madore & Freedman, 2020) | -6.20 | 0.01 | 0.04
$\mu_{0}(JAGB)_{J}$ | 24.97 | 0.02 | 0.04

To detect the $I$-band TRGB, we followed the CCHP procedure standardized in
Hatt et al. (2017). First, the $I$-band luminosity function was binned at the
0.01 mag level, and then smoothed using a non-truncated Gaussian kernel
smoothing algorithm called GLOESS (Gaussian-windowed, Locally-weighted
Scatterplot Smoothing) (Persson et al., 2004; Monson et al., 2017) in order to
reduce Poisson noise peaks (see the middle panel of Figure 3). We chose a
smoothing parameter of $\sigma_{s}=0.075$ for the luminosity function. Then
the smoothed luminosity function was convolved with a Sobel kernel [-1,0,+1]
(first derivative) edge detection filter, producing a gradient approximation
which has also been weighted by signal-to-noise. The resulting edge response
shows where the slope of the luminosity function is the steepest, or where the
first derivative of the luminosity function is the largest (see the right-most
panel of Figure 3). The maximum of this edge response function marks the TRGB.

Figure 3 demonstrates our extremely clear tip detection in the $I$ versus
$(V-I)$ data (distinct enough to be visually apparent to a few hundredths of a
magnitude). We base our systematic and statistical error associated with our
TRGB measurement on the artificial stars experiments performed in a previous
CCHP study, (Hatt et al., 2017). From their analyses, we have adopted a
conservative estimate of the systematic error on our measurement of the TRGB,
$\sigma_{sys}=0.01$ mag, and a conservative estimate of $\sigma_{stat}=0.01$
mag based on the width of our narrow edge-response function. After including
errors from the calibration procedure, we present a TRGB measurement in the
optical of $m_{F814W}^{TRGB}=20.93$ $\pm 0.01$ (stat) $\pm 0.03$ (sys) mag.
Adopting the extinction correction of $A_{I}=0.057$ $\pm 0.02$ (sys) mag and
an absolute magnitude of $M_{F814W}=-4.054$ $\pm 0.022$ (stat) $\pm 0.039$
(sys) mag from Freedman et al. (2020), we determine the distance modulus based
on the $I$-band TRGB to be $\mu_{0}(TRGB_{F814W})=24.93$ $\pm 0.02$ (stat)
$\pm 0.06$ (sys) mag.

###  3.2 Near-Infrared TRGB

Figure 4: Color Magnitude Diagrams for WLM in the Optical and Near-Infrared.
We cross-identified 15 stars (blue points) near the TRGB in the $I$ band where
our tip detection was the clearest, and then mapped them into $JHK$. We then
fit zeropoints using pre-determined slopes by minimizing the scatter for $J$
vs. $(J-K)$, $H$ vs. $(J-K)$, $K$ vs. $(J-K)$

Since the conception of the TRGB method with the use of a Sobel filter as
described in Lee et al. (1993), the majority of subsequent studies have used
the optical $I$ band in their endeavors as the $I$-band TRGB is ‘flat’ in
magnitude and is largely independent in metallicity and color. On the other
hand, whereas the $I$-band TRGB only has a unique zero point because its color
dependence is 0, the near-infrared (NIR) TRGB is upward-sloping with color,
and therefore has both a zero-point and slope term.

The NIR TRGB is a relatively recently developed distance indicator. Only a few
studies to date have employed the use of the NIR TRGB to measure distances to
galaxies. Dalcanton et al. (2012) first published NIR CMDs, indicating a
promising future for the NIR TRGB as a distance indicator. Wu et al. (2014)
subsequently undertook a similar analysis and came to the same conclusion.
Górski et al. (2011) derived NIR TRGB distances to five galaxies in the local
group. Hoyt et al. (2018) derived a NIR absolute calibration for the NIR TRGB
for the 2MASS system, and Madore et al. (2018) subsequently used their
calibration to publish a distance to IC 1613 in a companion paper. Finally,
Durbin et al. (2020) recently developed a new method for determining the NIR
TRGB by fitting an n-dimensional Gaussian to photometry of TRGB stars.

Although use of the NIR TRGB is relatively new, there are two distinct
advantages of employing the NIR TRGB over the optical. The first is that the
TRGB is about 2 magnitudes brighter in the NIR than the optical, and the
second is that line-of-sight extinction effects decrease for redder
wavelengths.

To detect the NIR TRGB, we follow the procedure in Freedman et al. (2020): we
first selected a set of “tracer tip stars” in the $I$ band, as it is flat and
our detection is well-defined. We then mapped those stars into the $JHK$ bands
to obtain NIR CMDs. Finally, we fit a TRGB to the tracer stars in the NIR
using predetermined slopes from (Hoyt et al., 2018) to determine the TRGB
zeropoint by minimizing the scatter for the zeropoint fit. This method is
effective because the same stars defining the tip in the flat $I$ band are
still the same stars that define the TRGB at other wavelengths.

We selected 15 tracer stars with an $I$-band mean equal to that of the peak
magnitude of the $TRGB_{I}$. These stars were then identified in the $JHK$
data by using topcat to match the catalogs by their $\alpha/\delta$. Next,
these stars were plotted in color-magnitude diagrams for $J,H$ and $K$ vs
$(J-K)$. Finally, a line was fit to the points for each CMD using the
predetermined slopes from the absolute calibration of Hoyt et al. (2018) using
a least-squares regression to find the zeropoint. We present the NIR TRGB fits
overlaid on their respective CMDs in Figure 4, along with the stars used for
mapping in the $I$ band. We found:

$m_{J}^{TRGB}=19.85-0.85[(J-K)_{o}-1.00]$ (1)
$m_{H}^{TRGB}=19.05-1.62[(J-K)_{o}-1.00]$ (2)
$m_{K}^{TRGB}=18.85-1.85[(J-K)_{o}-1.00]$ (3)

Our fitted zeropoints had errors on the mean of 0.036, 0.038, and 0.036 mag
for the fits on $JH$ and $K$, respectively.

We determined distance moduli based on absolute calibrations from Hoyt et al.
(2018), which we repeat here:

$M_{J}^{TRGB}=-5.14-0.85[(J-K)_{o}-1.00]$ (4)
$M_{H}^{TRGB}=-5.94-1.62[(J-K)_{o}-1.00]$ (5)
$M_{K}^{TRGB}=-6.14-1.85[(J-K)_{o}-1.00]$ (6)

The quoted errors on the zeropoints are $\pm 0.01$ mag (stat) and $\pm 0.06$
mag (sys).

We then applied extinction corrections, and factored in the average
photometric error for each of the filters near the TRGB: 0.017, 0.022, and
0.030 mag for the $JH$ and $K$ filters respectively. Our final NIR TRGB
distance moduli are: $\mu_{0}(TRGB_{J})=24.96\pm 0.04$ (stat) $\pm 0.06$ (sys)
mag, $\mu_{0}(TRGB_{H})=24.97\pm 0.04$ (stat) $\pm 0.07$ (sys) mag,
$\mu_{0}(TRGB_{K})=24.98\pm 0.04$ (stat) $\pm 0.08$ (sys) mag. Averaging these
values gives a NIR TRGB value of $\mu_{0}(TRGB_{NIR})=24.98\pm 0.04$ (stat)
$\pm 0.07$ (sys) mag.

##  4 The JAGB Method

The JAGB method is a precise and accurate distance indicator recently
introduced in Madore & Freedman (2020), which followed on that of Nikolaev &
Weinberg (2000); Weinberg & Nikolaev (2001). Recently, an independent study by
Ripoche et al. (2020) also indicated that the peak of the JAGB luminosity
function had to the potential to be a powerful standard candle. JAGB stars are
thermally-pulsating AGB stars. During pulsations, the convective envelope
penetrates deep into the star, bringing material to the surface. Each dredge-
up phase penetrates increasingly deeper, and by the star’s third-dredge-up
phase (and later), the convective envelope brings carbon from the He-burning
shell to the stellar surface (Habing & Olofsson, 2003).

Because of the carbon in JAGB stars’ atmospheres, they are photometrically
distinct from and redder than oxygen-rich AGBs in color-magnitude space, and
in the NIR they have well-defined limits in color and magnitude. To delineate
JAGB stars in color, Madore & Freedman (2020) chose conservative color cuts of
$1.30<(J-K)<2.00$ mag, following Nikolaev & Weinberg (2000). Stars bluer than
$(J-K)=1.30$ mag are primarily oxygen-rich AGB stars, but also include carbon
core-burning TP-AGBs that have not yet ‘dredged-up’ enough carbon to become
JAGB stars. Stars redder than $(J-K)=2.00$ mag are extreme carbon stars with
winds in which dust is being formed, leading to the star being increasingly
obscured by its own circumstellar envelope (Nikolaev & Weinberg, 2000). The
third and later dredge-up phases of older and less-massive stars are not
sufficiently deep to dredge up carbon; while the younger and more massive
stars have such hot convective envelopes that the carbon is burned into
nitrogen before it has time to reach the surface (Marigo & Girardi, 2007;
Habing & Olofsson, 2003). From an empirical perspective, there exists an
intermediate range of masses (and corresponding luminosities) for AGB stars
that have carbon-enriched atmospheres that we see as JAGB stars. As
demonstrated in Madore & Freedman (2020), the above physical constraints
manifest themselves in the NIR as JAGB stars, having a low-dispersion mean
magnitude in the $J$ band, which is independent of color, and therefore
qualifies as a standard candle.

The JAGB method has the potential to be an important and widely applicable
distance indicator. Weinberg & Nikolaev (2001) found that most stars within
the JAGB color limits of $1.30<(J-K)<2.00$ are variable. The variability of
JAGB stars increases the scatter in the JAGB luminosity function, but does not
contribute a systematic error (Freedman & Madore, 2020). The statistical error
can be lowered significantly by using a large number of JAGB stars.
Empirically, the distribution of JAGB stars was found to have an intrinsic
dispersion of only $\pm$0.20 mag (Weinberg & Nikolaev, 2001; Freedman &
Madore, 2020). The statistical error on the mean magnitude then decreases as
$\pm 0.20/\sqrt{N_{c}}$ mag, where $N_{c}$ is the number of JAGB stars in the
observed JAGB luminosity function. Additionally, in the $J$ band, systematic
errors due to metallicity and reddening effects appear to be small (Freedman &
Madore, 2020).

Furthermore, the JAGB method has several advantages over other local distance
indicators, including the TRGB method and the Cepheid Leavitt law. First, in
the near-infrared, JAGB stars are about one magnitude brighter than the TRGB
in the J band, so they can be detected out to appreciably greater distances.
Second, the JAGB method needs only one epoch of observations (like the TRGB),
whereas the Cepheids require upward of a dozen observations, spread out over
many months, in order to detect variability and measure periods, amplitudes,
mean magnitudes and mean colors. Finally, JAGB stars are found in all types of
galaxies, unlike Cepheids, which are only found in the star-forming disks of
late-type spiral and irregular galaxies. Therefore, the JAGB method can be
applied to a heterogeneous sample of type Ia supernova host galaxies,
including lenticulars and ellipticals.

###  4.1 Measurement of the JAGB Magnitude for WLM

Figure 5: Color Magnitude Diagram (left) and the smoothed luminosity function
from a bin size of 0.05 mag (right) of the JAGB star distribution. This data
is based on two epochs of observing and are temporally averaged. The
horizontal line represents the mean $J=18.80$ mag. Vertical lines show the
color cuts of $1.30<(J-K)<2.00$ mag made on the LF in order to avoid oxygen-
rich AGB stars on the blue side, and K branch stars on the red side. The two
horizontal lines at $J=19.18$ and $J=18.48$ mag represent upper and lower
magnitude limits for the JAGB stars distribution. The blue points represent
the 30 stars time-averaged with the Valcheva et al. (2007) near-infrared
dataset of AGB/carbon stars in their study of WLM.

WLM was identified as a good candidate for the JAGB method because of its high
number of carbon stars, quantified by a high carbon to M star ratio (C/M).
Galaxies with high C/M ratios are also known to have low metallicities (Cook
et al., 1986). WLM’s low metallicity of [Fe/H]$=-1.45\pm 0.2$ dex (Minniti &
Zijlstra, 1997) made it a good candidate for measuring both a TRGB distance in
the halo and a JAGB distance in the outer disk.

In Figure 5, we show the color-magnitude diagram of our $J$ vs $(J-K)$ data.
The two epochs of FourStar data, which were taken about a month apart (on
September 7, 2011 and October 5, 2011), have been averaged together. Grey
lines show the adopted color cut-offs used to select JAGB stars; i.e.,
$1.30<(J-K)<2.00$ mag. In Paper II, Freedman & Madore (2020) measured distance
moduli using the mean J-band magnitudes of the JAGB populations. We adopt that
same measure here, while noting that future studies of galaxies with larger
statistics will be needed to explore the possible effects of star formation
histories and metallicities, for example, on the higher moments and multi-
modal substructure in the overall JAGB luminosity function. In WLM we only
have 59 JAGB stars selected from Figure 5, centered on a mean value of $J=$
18.83 $\pm$ 0.35 mag, where the upper and lower limits are adopted from
Weinberg & Nikolaev (2001) based on several thousand, single-epoch
observations of JAGB stars in the LMC.

Finally, we note that Battinelli & Demers (2005) identified 77 carbon stars in
WLM observations from Valcheva et al. (2007), which we cross-matched to our
dataset. After applying our JAGB color and magnitude cuts to their data, we
identified 30 stars in common between our two datasets. In Figure 5, we have
plotted these stars in blue. We have averaged our FourStar $J$ and $K$
magnitudes with the Valcheva et al. (2007) data.

The mean apparent magnitude of the JAGB distribution is $m_{J}=18.80$ mag,
with a dispersion of $\pm 0.17$ mag. Thus, the 59 JAGB stars contribute an
error on the mean of $\pm 0.17/\sqrt{59}=0.02$ mag. The mean photometric error
of stars in the $J$ band in the JAGB distribution yields an error of
$\sigma_{sys}=0.007$ mag. After applying an extinction correction of
$A_{J}=0.027$ mag, and using the absolute magnitude calibration from Madore &
Freedman (2020) of $M_{J}=-6.20\pm 0.01$ (stat) $\pm 0.04$ (sys) mag, we
determine a distance modulus of $\mu_{0}$ (JAGB) $=24.97\pm 0.02$ (stat) $\pm
0.04$ (sys) mag.

##  5 Independent Distance Comparisons

Figure 6: Previous distance modulus estimates to WLM. Filled circles represent
as published values and open circles represent re-calculated standardized
values using $M_{F814W}=-4.054$ (Freedman et al., 2020) and $A_{I}-0.057$
(NED). Papers from which these values were taken from are compiled in Table 4.
Gold stars represent values from this study.

In this section, we compare our JAGB and TRGB distances against previously
published, independent measurements. A compilation of these distance moduli
and their references can be found in Table 4. A visual compilation can be
found in Figure 6; filled circles represent original published values and
their error bars, and open circles represent re-calculated values standardized
to the Freedman et al. (2020) TRGB calibration in the $I$ band and a
foreground extinction value of $A_{I}=0.057$ mag. Finally, we present a
Cepheid distance comparison in Section 5.3.

###  5.1 TRGB Distance Comparisons

Since 1993, 11 TRGB distances have been measured in the $I$ band, one in the
$J$ band, and one in the $K$ band for WLM. In Figure 6, we show our TRGB
distance moduli compared with previous papers’ measurements. The open dots
represent re-calculated distance moduli normalized to the Freedman et al.
(2020) $I$-band calibration and a foreground extinction correction of
$A_{I}=0.057$ mag (NED). The filled dots represent the original values and
their errors quoted in the paper.

We find the average value of the 11 published $TRGB_{I}$ distance moduli to be
$\mu_{0}=24.87\pm 0.08$ mag. This is in agreement, at the 1-$\sigma$ level,
with our measured value of $\mu_{0}(TRGB_{F814W})=24.93\pm 0.02$ (stat) $\pm
0.06$ (sys) mag.

For the NIR TRGB, Górski et al. (2011) measured a near-infrared TRGB distance
to WLM; their average value of $25.13\pm 0.16$ mag differs by about 1-$\sigma$
from our measured value of $24.98\pm 0.04$ (stat) $\pm 0.08$ (sys) mag.

###  5.2 JAGB Distance Comparison

Only one previous study has used the JAGB method to determine the distance to
WLM (Freedman & Madore, 2020). They located known carbon stars found by
Battinelli & Demers (2005) in $JK$ data obtained by Valcheva et al. (2007)
using ESO’s New Technology Telescope. They found $J=18.80\pm 0.25$ mag, with a
true distance modulus of $\mu_{0}$ (JAGB) $=24.97\pm 0.05$ (stat) mag agreeing
(likely fortuitously, exactly) with our value of $\mu_{0}$ (JAGB) $=24.97\pm
0.02$ (stat) $\pm 0.04$ (sys) mag.

###  5.3 Cepheid Distance Comparison

Figure 7: GLOESS-fit light curves in $VIJHK$ when available for the 60
Cepheids from Pietrzyński et al. (2007). Bands are distinguished by color, as
labelled next to C001 at the top left. Cepheid numbers from the original
publication are labeled in the upper-left corner of each subplot, and our
quality rankings (where 5 = highest quality; 0 = lowest quality) are shown in
the lower left. Red quality numbers indicate Cepheids excluded from the final
sample. The light curve data shown in this figure are available as the Data
behind the Figure.  Figure 6. cont:

Sandage & Carlson (1985) (SC85) first reported the discovery of 15 short
period Cepheids in WLM. These discoveries were based on 30 epochs of deep
photographic $B$ band observations from the Palomar 5m and duPont 2.5m
telescopes. Despite substantial effort spanning over 30 years of observation,
their longest-discovered period was 9.6 days, and many Cepheids were close to
the $B=22.8$ mag detection limit of their photographic plates. It was another
two decades before WLM was subject to a variable star search using modern
detectors. Pietrzyński et al. (2007) (P07) commenced a follow up search in the
$V$ and $I$ bands using the Warsaw 1.3m telescope, and Gieren et al. (2008)
(G08) followed up in the $J$ and $K$ bands using infrared cameras on the New
Technology and Baade Magellan telescopes. P07 found 60 Cepheids, recovering
all but 2 of the original 15 discovered in SC85. Notably, they reported the
first long-period Cepheid in WLM. SC85 originally reported Cepheid V12 with a
period of 7.94 days, but P07 revised this to a significantly longer period of
54.7 days. G08 reported sparsely sampled (typically 3) $J$ and $K$
observations for 24 of the Cepheids in the P07 sample. We were able to add 6
new $JHK$ observations to each of 12 Cepheids taken on two nights in 2011 and
2019, as discussed above in Section 2.2.1. Additionally, Valcheva et al.
(2007) reported single phase points in the $J$ band for 5 Cepheids.

Merging these data sets and adopting the periods determined by P07, we GLOESS
fit the light curves. We automated the selection of the GLOESS smoothing
parameter $\sigma_{s}$, assigning it based on the number of points and the
average error of those points in each band. These initial fits were visually
inspected, and bands with $\sigma_{s}$ that significantly under-fit the data
were refined by hand.

Figure 9: Period Luminosity relations in $VIJHK$ for WLM Cepheids with quality
rankings of 4 or 5. The legend displays fits to the fixed slopes of Monson et
al. (2012) in each band. Open circles indicate excluded data. It is clear that
our quality cut closely resembles a period cut at $\log P\sim 0.8$.

Each of the Cepheids was ranked on a quality scale of 0 to 5 according to the
amount of phase coverage, overall data quality, and the GLOESS goodness-of-fit
across the bands. No period cuts were applied, but the quality rankings
tracked period and apparent magnitude fairly closely, as might be expected.
All 60 Cepheids are shown with our GLOESS fits and quality rankings in Figure
7.

From the final GLOESS fits, we calculated intensity-averaged mean magnitudes
in all five bands. However, we should note that the archival $V$ band data
from P07 was discovered to be incomplete when compared to the published light
curve plots. Due to a server failure, the rest of the full data set was
unavailable from the authors when contacted. We elected to still use the
incomplete $V$ data since it is sufficient for calculating mean magnitudes,
and we obtain similar values to the original publication.

We constructed five sets of multi-wavelength PL relations by performing
progressive quality cuts. The largest set included all Cepheids, and the
smallest data set contained only the 14 best Cepheids (ranked 4 and 5). We
used this highest-quality subsample for determining our final distance modulus
relations.

We fit the fixed slopes of Monson et al. (2012) to our Cepheid subsample and
used their Milky Way intercepts to determine wavelength-dependent apparent
distance moduli, plotted in Figure 9. Moduli and corresponding errors for all
quality cutoffs are given in Table 3. The Cardelli et al. (1989) extinction
curve was fit to the distance moduli using a chi-squared grid minimizer,
allowing the true distance modulus and reddening to vary.

Table 3: Cepheid Distance Moduli Dependence on Quality Cutoffs Quality cutoff | $\mu_{V}$ | $\mu_{I}$ | $\mu_{J}$ | $\mu_{H}$ | $\mu_{K}$ | $\mu_{0}$ | $E(B-V)$
---|---|---|---|---|---|---|---
all data | $25.03\pm 0.08$ | $24.90\pm 0.08$ | $24.97\pm 0.06$ | $24.82\pm 0.05$ | $24.93\pm 0.05$ | $24.96\pm 0.05$ | $-0.02\pm 0.07$
$\geq 1$ | $25.03\pm 0.09$ | $24.91\pm 0.09$ | $24.97\pm 0.07$ | $24.85\pm 0.05$ | $24.97\pm 0.05$ | $24.96\pm 0.05$ | $-0.01\pm 0.07$
$\geq 2$ | $25.03\pm 0.07$ | $24.89\pm 0.10$ | $24.96\pm 0.06$ | $24.94\pm 0.05$ | $25.04\pm 0.06$ | $24.94\pm 0.04$ | $0.03\pm 0.06$
$\geq 3$ | $25.00\pm 0.08$ | $24.89\pm 0.10$ | $24.97\pm 0.06$ | $25.00\pm 0.04$ | $25.10\pm 0.04$ | $24.92\pm 0.03$ | $0.06\pm 0.05$
$\boldsymbol{\geq 4}$ | $\boldsymbol{24.98\pm 0.07}$ | $\boldsymbol{25.02\pm 0.08}$ | $\boldsymbol{25.05\pm 0.06}$ | $\boldsymbol{25.06\pm 0.03}$ | $\boldsymbol{25.13\pm 0.04}$ | $\boldsymbol{24.98\pm 0.03}$ | $\boldsymbol{0.05\pm 0.04}$
Figure 10: Extinction curve fit to the wavelength-dependent distance moduli
derived from WLM Cepheids. Error bars on the individual points give the
standard error on the mean as determined from the PL-relation. The solid line
shows the best fit, and the dashed lines give the 1$\sigma$ error. The inset
in the lower right corner shows the 1, 2, and 3$\sigma$ contours of the error
space of the fit, with the true distance modulus on the x-axis and the
reddening on the y-axis.

Due to the relatively large errors on the individual distance moduli compared
to the (minimized) scatter in the finally-adopted fit (see Figure 10), we
chose to estimate the fit error from the data point uncertainties alone. We
did this by assuming a Gaussian distribution for each error, with the width
equal to the uncertainty on the point and the mean equal to the distance from
the fit. We summed the distributions and calculated the plus and minus offsets
from the mean that would each contain 34% of the area. The offsets were equal
in our chosen subsample, indicating a symmetrical distribution, and we
determined a statistical error on the true distance modulus to be $\pm 0.03$
mag. Using this error value, the 1, 2, and 3$\sigma$ contours of the error
space were calculated (see inset in Figure 10). We determined the $E(B-V)$
error by dividing half the full 1$\sigma$ contour range by the square root of
the number of data points. This results in an error of $\pm 0.04$ mag for our
subsample.

The Cepheid Leavitt law distance modulus and reddening determined most
recently by the Araucaria Group (Gieren et al., 2008) were $\mu_{0}=24.92\pm
0.04$ mag and $E(B-V)=0.08\pm 0.02$ mag. Our finally-adopted distance modulus
agrees well with this earlier determination, with $\mu_{0}=24.98\pm 0.03$
(stat) $\pm 0.04$ (sys) mag. With this fit, we obtain a total line-of-sight
reddening value of $E(B-V)=0.05\pm 0.04$ mag. Adopting an approximate value of
0.03 mag for the foreground galactic reddening to WLM (Schlafly & Finkbeiner,
2011), we obtain an intrinsic reddening of $0.02\pm 0.04$ mag.

Table 4: Previously Published TRGB and JAGB Distance Moduli to WLM

Study | Filter & Method |  $\mu_{0}$ (published) [mag] |  $\mu_{0}$ (standardized) [mag] | Notes
---|---|---|---|---
Lee et al. (1993) | $I_{~{}TRGB}$ | $24.81\pm 0.10$ | $24.85$ |
Minniti & Zijlstra (1997) | $I_{~{}TRGB}$ | $24.75\pm 0.10$ | 24.80 |
Salaris & Cassisi (1997) | $I_{~{}TRGB}$ | $24.97$ | $24.85$ | No errors given
Ferrarese et al. (2000) | $I_{~{}TRGB}$ | $24.84$ | $-$ | No errors given, no original $I_{TRGB}$ given
McConnachie et al. (2005) | $I_{~{}TRGB}$ | $24.85\pm 0.08$ | $24.85$ |
Rizzi et al. (2007) | $I_{~{}TRGB}$ | $24.93\pm 0.04$ | $24.92$ |
Jacobs et al. (2009) | $I_{~{}TRGB}$ | $24.95^{+0.03}_{-0.09}$ | $24.95$ |
McCall et al. (2012) | $I_{~{}TRGB}$ | $24.82\pm 0.10$ | $24.85$ |
McQuinn et al. (2017) | $I_{~{}TRGB}$ | $24.94\pm 0.03$ | $24.91$ |
Albers et al. (2019) | $I_{~{}TRGB}$ | $24.93\pm 0.07$ | $24.91$ |
Górski et al. (2011) | $J_{~{}TRGB}$ | $25.14\pm 0.15$ | $-$ | No original $J_{~{}TRGB}$ given
| $K_{~{}TRGB}$ | $25.12\pm 0.17$ | $-$ | No original $K_{~{}TRGB}$ given
Freedman & Madore (2020) | $JAGB$ | $24.97\pm 0.05$ | $24.97$ |

Note. — The $\mu_{0}$ ($TRGB_{I}$) has been standardized to
$M_{I}^{TRGB}=-4.05$ mag from Freedman et al. (2020) and an extinction of
$A_{I}=0.057$ mag.

###  5.4 Additional Distance Comparisons

Finally, we compare our four measured distance moduli to the most precise
previously published distances to WLM (obtained via NED) calculated with
methods not used in this study. These measurements are based on the Horizontal
Branch (HB) feature and CMD fitting techniques.

Rejkuba et al. (2000) measured a distance modulus of $\mu_{0}=24.95\pm 0.13$
mag based on the $V$-band magnitude of the HB, which is in 1-$\sigma$
agreement with all four of our distance measurements.

There have also been three measured distance moduli to WLM based on CMD-
fitting analysis techniques (Hodge et al., 1999; Dolphin, 2000; Weisz et al.,
2014). Averaging them yields a distance modulus of $\mu_{0}=24.83\pm 0.11$
mag, which is $1-2$ $\sigma$ from all four of our measured distance moduli.

##  6 Summary and Conclusions

In this paper we have presented consistent measured distances to WLM using the
optical TRGB, NIR TRGB, JAGB method, and $VIJHK$ observations of Cepheids,
demonstrating promise for the JAGB method as a competitive supernovae type Ia
calibrator. We have also showed the JAGB’s future effectiveness as an
independent cross-check for systematics with the Cepheid and TRGB methods at
NIR wavelengths.

We determined the optical TRGB distance modulus, measured in the $HST$ $F814W$
filter to be $\mu_{0}(TRGB_{F814W})=24.93\pm 0.02$ (stat) $\pm 0.06$ (sys)
mag. Our second distance modulus, determined using the NIR TRGB, was
$\mu_{0}(TRGB_{NIR})=24.98\pm 0.04$ (stat) $\pm 0.07$ (sys) mag. Third, we
calculated our value for the JAGB distance modulus to be $\mu_{0}$ (JAGB)
$=24.97\pm 0.02$ (stat) $\pm 0.04$ (sys) mag. Finally, we calculated a Cepheid
distance modulus of $\mu_{0}$ (Cepheid)$=24.98\pm 0.03\text{ (stat)}\pm
0.04\text{ (sys) mag}$. These correspond to distances of $968\pm 27$ kpc,
$991\pm 37$ kpc, $986\pm 21$ kpc, and 991 $\pm$ 23 kpc for the optical TRGB,
NIR TRGB, JAGB, and Cepheid-based methods, respectively.

The consistency of the distance moduli calculated across the optical TRGB,
near-infrared TRGB, JAGB methods, and the Cepheid Leavitt law adds to the
growing evidence that the JAGB method is an excellent distance indicator and
has promising potential for future work. The JAGB method is accurate, as shown
by its consistency with the TRGB and Cepheids in its distance modulus
measurement. In addition, the JAGB method is also precise, with error bars
comparable to those based on optical measurements of the TRGB and multi-
wavelength measurements of Cepheids. Finally, the JAGB method is simple; only
the mean magnitude of the JAGB distribution needs to be measured.

The agreement in the measured TRGB and Cepheid distances for WLM is consistent
with that found for other nearby galaxies; however, the two methods
increasingly disagree out to further distances (Freedman et al., 2019),
indicating possible systematic errors in one or both methods. In the future,
it it will be helpful to have independent JAGB distances for galaxies at
farther distances for which the agreements starts to break down.

The full potential of the JAGB method as a distance indicator is still to be
fully realized. Future studies using HST and eventually JWST data are capable
of measuring distances well beyond the Local Group. We note that JAGB stars
are about one mag brighter than the TRGB in the infrared (see Figure 5).
Although long-period Cepheids and Mira variables are brighter, they need
increasingly longer baselines in order to secure their periods, ranging from
100 to up to 1000 days or more, which can become observationally expensive.
Our expectation is that JAGB stars will be able to calibrate type Ia
supernovae hosts out to 100 Mpc.

With future work and development and a large sample of JAGB distances to
supernovae type Ia hosts, the JAGB method will eventually result in its own
determination of the Hubble Constant using near-infrared data from HST and
JWST. Furthermore, the JAGB method can be used in tandem with the NIR TRGB and
Cepheid-based methods on galaxies with low-reddening and low surface-
brightness such as WLM; as shown in this study, the two methods can be used as
a cross-checks and comparison for each other for revealing systematic errors.

We thank Peter Stetson for providing us with his copies of daophot, daogrow,
and allframe. We thank Violet Mager for providing useful comments on flat-
fielding IMACS data. We thank Dylan Hatt for his help with the preparation of
a figure. Finally, we thank the anonymous referee for their constructive
suggestions. Some of the data presented in this paper were obtained from the
Mikulski Archive for Space Telescopes (MAST). STScI is operated by the
Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-26555. This research has made use of the NASA/IPAC
Extragalactic Database (NED) and the NASA/IPAC infrared Science Archive
(IRSA), both of which are operated by the Jet Propulsion Laboratory,
California Institute of Technology, under contract with the National
Aeronautics and Space Administration. This paper includes data gathered with
the 6.5 meter Baade/Magellan Telescope located at Las Campanas Observatory,
Chile. This publication makes use of data products from the Two Micron All Sky
Survey, which is a joint project of the University of Massachusetts and the
Infrared Processing and Analysis Center/California Institute of Technology,
funded by the National Aeronautics and Space Administration and the National
Science Foundation. HST (ACS/WFC), Magellan-Baade (IMACS, FourStar)

## References

  * Albers et al. (2019) Albers, S. M., Weisz, D. R., Cole, A. A., et al. 2019, MNRAS, 490, 5538
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Battinelli & Demers (2005) Battinelli, P. & Demers, S. 2005, A&A, 434, 657
  * Beaton et al. (2019) Beaton, R. L., Seibert, M., Hatt, D., et al. 2019, ApJ, 885, 141
  * Bohlin (2016) Bohlin, R. C. 2016, AJ, 152, 60
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Cook et al. (1986) Cook, K. H., Aaronson, M., & Norris, J. 1986, ApJ, 305, 634
  * Dalcanton et al. (2012) Dalcanton, J. J., Williams, B. F., Melbourne, J. L., et al. 2012, ApJS, 198, 6
  * Dolphin (2000) Dolphin, A. E. 2000, ApJ, 531, 804.
  * Dressler et al. (2011) Dressler, A., Bigelow, B., Hare, T., et al. 2011, PASP, 123, 288
  * Durbin et al. (2020) Durbin, M. J., Beaton, R. L., Dalcanton, J. J., et al. 2020, ApJ, 898, 57
  * Ferrarese et al. (2000) Ferrarese, L., Mould, J. R., Kennicutt, R. C., et al. 2000, ApJ, 529, 745
  * Freedman et al. (2001) Freedman, W. L., Madore, B. F., Gibson, B. K., et al. 2001, ApJ, 553, 47.
  * Freedman et al. (2019) Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34
  * Freedman et al. (2012) Freedman, W. L., Madore, B. F., Scowcroft, V., et al. 2012, ApJ, 758, 24.
  * Freedman & Madore (2020) Freedman, W. L. & Madore, B. F. 2020, ApJ, 899, 67
  * Freedman et al. (2020) Freedman, W. L., Madore, B. F., Hoyt, T., et al. 2020, ApJ, 891, 57
  * Gieren et al. (2008) Gieren, W., Pietrzyński, G., Szewczyk, O., et al. 2008, ApJ, 683, 611
  * Górski et al. (2011) Górski, M., Pietrzyński, G., & Gieren, W. 2011, AJ, 141, 194
  * Habing & Olofsson (2003) Habing, H. J. & Olofsson, H. 2003, Asymptotic giant branch stars
  * Hatt et al. (2017) Hatt, D., Beaton, R. L., Freedman, W. L., et al. 2017, ApJ, 845, 146
  * Hodge et al. (1999) Hodge, P. W., Dolphin, A. E., Smith, T. R., et al. 1999, ApJ, 521, 577.
  * Hoyt et al. (2018) Hoyt, T. J., Freedman, W. L., Madore, B. F., et al. 2018, ApJ, 858, 12
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Iben (1973) Iben, I. 1973, ApJ, 185, 209.
  * Iben & Renzini (1983) Iben, I. & Renzini, A. 1983, ARA&A, 21, 271.
  * Jacobs et al. (2009) Jacobs, B. A., Rizzi, L., Tully, R. B., et al. 2009, AJ, 138, 332
  * Jang et al. (2020) Jang, I. S., Hoyt, T., Beaton, R., et al. 2020, arXiv:2008.04181
  * Krist et al. (2011) Krist, J. E., Hook, R. N., & Stoehr, F. 2011, in Proc. SPIE, Vol. 8127, Optical Modeling and Performance Predictions V, 81270J
  * Lee et al. (1993) Lee, M. G., Freedman, W. L., & Madore, B. F. 1993, ApJ, 417, 553
  * Madore et al. (2018) Madore, B. F., Freedman, W. L., Hatt, D., et al. 2018, ApJ, 858, 11
  * Madore & Freedman (2020) Madore, B. F. & Freedman, W. L. 2020, ApJ, 899, 66
  * Madore et al. (2009) Madore, B. F., Mager, V., & Freedman, W. L. 2009, ApJ, 690, 389
  * Marigo & Girardi (2007) Marigo, P. & Girardi, L. 2007, A&A, 469, 239.
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77.
  * McCall et al. (2012) McCall, M. L., Vaduvescu, O., Pozo Nunez, F., et al. 2012, A&A, 540, A49
  * McConnachie et al. (2005) McConnachie, A. W., Irwin, M. J., Ferguson, A. M. N., et al. 2005, MNRAS, 356, 979
  * McQuinn et al. (2017) McQuinn, K. B. W., Boyer, M. L., Mitchell, M. B., et al. 2017, ApJ, 834, 78
  * Melotte (1926) Melotte, P. J. 1926, MNRAS, 86, 636.
  * Minniti & Zijlstra (1997) Minniti, D., & Zijlstra, A. A. 1997, AJ, 114, 147
  * Monson et al. (2017) Monson, A. J., Beaton, R. L., Scowcroft, V., et al. 2017, AJ, 153, 96
  * Monson et al. (2012) Monson, A. J., Freedman, W. L., Madore, B. F., et al. 2012, ApJ, 759, 146
  * Nikolaev & Weinberg (2000) Nikolaev, S., & Weinberg, M. D. 2000, ApJ, 542, 804
  * Persson et al. (2004) Persson, S. E., Madore, B. F., Krzemiński, W., et al. 2004, AJ, 128, 2239
  * Persson et al. (2013) Persson, S. E., Murphy, D. C., Smee, S., et al. 2013, PASP, 125, 654
  * Pietrzyński et al. (2007) Pietrzyński, G., Gieren, W., Udalski, A., et al. 2007, AJ, 134, 594
  * Rejkuba et al. (2000) Rejkuba, M., Minniti, D., Gregg, M. D., et al. 2000, AJ, 120, 801
  * Riess et al. (2019) Riess, A. G., Casertano, S., Yuan, W., et al. 2019, ApJ, 876, 85.
  * Riess et al. (2016) Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56
  * Ripoche et al. (2020) Ripoche, P., Heyl, J., Parada, J., et al. 2020, MNRAS, 495, 2858
  * Rizzi et al. (2007) Rizzi, L., Tully, R. B., Makarov, D., et al. 2007, ApJ, 661, 815
  * Sackmann et al. (1974) Sackmann, I.-J., Smith, R. L., & Despain, K. H. 1974, ApJ, 187, 555.
  * Salaris & Cassisi (1997) Salaris, M., & Cassisi, S. 1997, MNRAS, 289, 406
  * Salaris & Cassisi (2005) Salaris, M., & Cassisi, S. 2005, Evolution of Stars and Stellar Populations, 1st edn. (Wiley)
  * Sandage & Carlson (1985) Sandage, A. R. & Carlson, G. 1985, AJ, 90, 1464
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
  * Sirianni et al. (2005) Sirianni, M., Jee, M. J., Benítez, N., et al. 2005, PASP, 117, 1049
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191
  * Stetson (1990) Stetson, P. B. 1990, PASP, 102, 932
  * Stetson (1994) —. 1994, PASP, 106, 250
  * Stetson (2000) —. 2000
  * Taylor (2005) Taylor, M. B. 2005, Astronomical Data Analysis Software and Systems XIV, 29
  * Valcheva et al. (2007) Valcheva, A. T., Ivanov, V. D., Ovcharov, E. P., et al. 2007, A&A, 466, 501
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Weinberg & Nikolaev (2001) Weinberg, M. D. & Nikolaev, S. 2001, ApJ, 548, 712
  * Weisz et al. (2014) Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014, ApJ, 789, 147.
  * Weisz (2014) Weisz, C. 2014, Completing the Census of Isolated Dwarf Galaxy Star Formation Histories, HST Proposal
  * Wolf (1909) Wolf, M. 1909, Astronomische Nachrichten, 183, 187.
  * Wu et al. (2014) Wu, P.-F., Tully, R. B., Rizzi, L., et al. 2014, AJ, 148, 7

##  Appendix A Photometric Quality Cuts

In order to limit our photometry to stellar sources and eliminate extended
sources, we applied cuts on the following photometry parameters provided by
daophot: the photometric uncertainty $\sigma$, the sharpness parameter, and
the chi parameter $\chi$. To be considered a stellar source, every source had
to pass all three cuts. In Figure A1, the black points show sources that
passed all three cuts, and the grey points are sources that failed one or more
of the cuts. We made cuts based on the following exponential functions as
function of magnitude $m$, where the chosen values for each given function and
filter can be found in Table 5:

$\sigma<a+b\times e^{m-c}$ (A1) $|sharp|<a+b\times e^{m-c}$ (A2)
$\chi<a+e^{-(m-b)}$ (A3)

Table 5: Photometric Quality Cuts  Filter |  $\sigma$ (a,b,c) |  $sharp$ (a,b,c) |  $\chi$ (a,b)
---|---|---|---
$I$ | 0.01,0.01,21.25 | 0.20, 70.0, 25.5 | 1.15, 20.5
$J$ | 0.01,0.003,17.5 | 0.15,0.01,15.5 | .35, 17.5
$K$ | 0.01,0.003,16.1 | 0.1,0.01,14.5 | .1,16
$H$ | 0.01,0.003,16.5 | 0.1,0.01,16 | .13,15
Figure A1: Photometric quality cuts made on $IJHK$ data. For the $JHK$ data, a
source had to pass all three $\sigma,sharp$, and $\chi$ cuts for all three
bands to be considered a stellar source.