# Structural parameters and possible association of the Ultra-Faint Dwarfs
Pegasus III and Pisces II from deep Hubble Space Telescope photometry

Hannah Richstein1, Ekta Patel2,3, Nitya Kallivayalil1, Joshua D. Simon4, Paul
Zivick5, Erik Tollerud6, Tobias Fritz1, Jack T. Warfield1, Gurtina Besla7,
Roeland P. van der Marel6,8, Andrew Wetzel9, Yumi Choi2,6, Alis Deason10,11,
Marla Geha12, Puragra Guhathakurta13, Myoungwon Jeon14, Evan N. Kirby15,16,
Mattia Libralato17, Elena Sacchi18,19, and Sangmo Tony Sohn6 1University of
Virginia, Department of Astronomy, 530 McCormick Road, Charlottesville, VA
22904, USA
2Department of Astronomy, University of California, Berkeley, 501 Campbell
Hall, Berkeley, CA 94720, USA
3Miller Institute for Basic Research in Science, 468 Donner Lab, Berkeley, CA
94720, USA
4Observatories of the Carnegie Institution for Science, 813 Santa Barbara
Street, Pasadena, CA 91101, USA
5Mitchell Institute for Fundamental Physics and Astronomy, Department of
Physics and Astronomy, Texas A & M University, 578 University Drive, College
Station, TX 77843, USA
6Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD
21218, USA
7Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson,
AZ 85721, USA
8Center for Astrophysical Sciences, The William H. Miller III Department of
Physics & Astronomy, Johns Hopkins University, 3400 N. Charles Street,
Baltimore, MD 21218, USA
9Department of Physics & Astronomy, University of California, Davis, One
Shields Avenue, Davis, CA 95616, USA
10Institute for Computational Cosmology, Department of Physics, University of
Durham, South Road, Durham DH1 3LE, UK
11Centre for Extragalactic Astronomy, Department of Physics, University of
Durham, South Road, Durham DH1 3LE, UK
12Department of Astronomy, Yale University, 52 Hillhouse Avenue, New Haven, CT
06520, USA
13UCO/Lick Observatory, Department of Astronomy & Astrophysics, University of
California Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA
14Department of Astronomy & Space Science, Kyung Hee University, 1732
Deogyeong-daero, Yongin-si, Gyeonggi-do 17104, Korea
15Department of Astronomy, California Institute of Technology, 1200 E.
California Blvd., MC 249-17, Pasadena, CA 91125, USA
16Department of Physics, University of Notre Dame, 225 Nieuwland Science Hall,
Notre Dame, IN 46556, USA
17AURA for the European Space Agency (ESA), Space Telescope Science Institute,
3700 San Martin Drive, Baltimore, MD 21218, USA
18Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16, 14482
Potsdam, Germany
19INAF–Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via
Gobetti 93/3, I-40129 Bologna, Italy

###### Abstract

We present deep Hubble Space Telescope (HST) photometry of the ultra-faint
dwarf (UFD) galaxies Pegasus III (Peg III) and Pisces II (Psc II), two of the
most distant satellites in the halo of the Milky Way (MW). We measure the
structure of both galaxies, derive mass-to-light ratios with newly determined
absolute magnitudes, and compare our findings to expectations from UFD-mass
simulations. For Peg III, we find an elliptical half-light radius of
$a_{h}{{=}}1\farcm 88^{+0.42}_{-0.33}$ ($118^{+31}_{-30}$ pc) and
$M_{V}{=}{-4.17}^{+0.19}_{-0.22}$; for Psc II, we measure $a_{h}{=}1\farcm
31^{+0.10}_{-0.09}$ ($69\pm 8$ pc) and $M_{V}{=}{-4.28}^{+0.19}_{-0.16}$. We
do not find any morphological features that indicate a significant interaction
between the two has occurred, despite their close separation of only $\sim$40
kpc. Using proper motions (PMs) from Gaia early Data Release 3, we investigate
the possibility of any past association by integrating orbits for the two UFDs
in a MW-only and a combined MW and Large Magellanic Cloud (LMC) potential. We
find that including the gravitational influence of the LMC is crucial, even
for these outer-halo satellites, and that a possible orbital history exists
where Peg III and Psc II experienced a close ($\sim$10–20 kpc) passage about
each other just over $\sim$1 Gyr ago, followed by a collective passage around
the LMC ($\sim$30–60 kpc) just under $\sim$1 Gyr ago. Considering the large
uncertainties on the PMs and the restrictive priors imposed to derive them,
improved PM measurements for Peg III and Psc II will be necessary to clarify
their relationship. This would add to the rare findings of confirmed pairs of
satellites within the Local Group.

###### Subject headings:

galaxies: dwarf; galaxies: stellar content; Local Group

##  1\. Introduction

Ultra-faint dwarf (UFD) galaxies are the most dark-matter-dominated systems
discovered, and thus a preferred laboratory for studying how well cosmological
models predict behavior on small scales. The widely accepted Lambda Cold Dark
Matter ($\Lambda$CDM) model uses the hierarchical accretion of low-mass
systems to explain the growth of dark matter halos (e.g., Navarro et al.,
1997). As UFDs occupy the least-massive dark matter halos discovered, they
could be considered remnants of this hierarchical process while themselves
having the simplest assembly histories.

The dark matter halos in which UFDs reside have extrapolated virial masses of
approximately 109 M⊙ (e.g., Strigari et al., 2008), about two orders of
magnitude smaller than where the $\Lambda$CDM theory predicts central dark
matter densities in apparent contrast with observations. For example, in the
mid-1990s, the cusp-core issue arose when observations of dwarf spiral
galaxies showed cored mass distributions rather than the cusps predicted by
dark-matter-only simulations (e.g., Flores & Primack, 1994; Moore, 1994). This
was also found to be the case for Local Group dwarf spheroidals such as Draco
and Fornax (e.g., Read & Gilmore, 2005; Goerdt et al., 2006), although other
analyses have suggested that dSphs are consistent with the expected dark
matter density profiles (e.g., Jardel & Gebhardt, 2013; Breddels & Helmi,
2013; Strigari et al., 2017; Read et al., 2018, 2019).

A limiting factor for observational constraints is that we are often
restricted to line-of-sight (LOS) velocity data of the stars residing in dwarf
galaxies, and with no information on the tangential velocity components, we
suffer from the mass-anisotropy degeneracy. It will thus take more
information, such as the shape of the velocity distribution or galactic
internal proper motions (PMs), in addition to radial velocities to be able to
distinguish a dark matter central core or cusp (e.g., Strigari et al., 2007;
Read et al., 2021; Guerra et al., 2021). Once we have full kinematic
information, current virial mass estimators can be extended to further
constrain UFD dark matter halo properties (Errani et al., 2018). In the
meantime, we can use more easily measured UFD properties such as half-light
radius and luminosity to explore the population as a whole, as well as their
simulated analogs.

Defined by having $M_{V}$ values fainter than $-7.7$ (e.g., Bullock & Boylan-
Kolchin, 2017; Simon, 2019), UFDs went undetected until the advent of large-
scale digital sky surveys, beginning with the Sloan Digital Sky Survey in 2005
(SDSS; York et al. 2000; Willman et al. 2005). Currently, over 21 UFDs have
been spectroscopically verified, and more than 20 other candidates have been
identified (e.g., Simon, 2019). These numbers are expected to increase further
with the beginning of data collection at the Vera C. Rubin Observatory.

While there is much to learn from studying any of these faint satellites,
examining the relatively isolated UFDs at larger Galactocentric distances is
particularly useful, as their current environments are more comparable to
those produced in most simulations that resolve down to the UFD-mass level
($M_{\mathrm{vir}}\simeq 10^{9}$ M⊙ at $z{=}0$) (e.g., Jeon et al., 2017,
2021a, 2021b; Wheeler et al., 2019). More recently, Applebaum et al. (2021)
and Grand et al. (2021) were also able to resolve down to UFD-mass scale using
cosmological MW zoom-in simulations. Comparing the properties of observed UFDs
to both of these types of simulations could help us to disentangle the effects
that host galaxies may have on their satellite UFDs.

Pegasus III (Peg III) is located at a heliocentric distance of approximately
215 kpc ($R_{\mathrm{GC}}{\sim}213$ kpc; Kim et al. 2016) and is thus one of
the most distant MW UFDs known. Peg III was reported in Kim et al. (2015),
having been discovered in Data Release 10 of SDSS (Ahn et al., 2014) and
confirmed with the Dark Energy Camera (DECam). The discovery team noted the
proximity ($\sim$30 kpc projected and $\sim$32 kpc line-of-sight separation)
of Peg III to Pisces II (Psc II; $R_{\odot}{\sim}$183 kpc,
$R_{\mathrm{GC}}{\sim}$182 kpc; Belokurov et al. 2010; Sand et al. 2012) at
the time and suggested the possibility of an association. In a follow-up paper
using Magellan/IMACS for photometry and Keck/DEIMOS for spectroscopy, Kim et
al. (2016) derived a radial velocity for Peg III that, in the Galactic
standard-of-rest (GSR), only differed from that of Psc II by $\sim$10 km s-1
($-67.6\pm 2.6$ and $-79.9\pm 2.7$ km s-1 (Kirby et al., 2015), respectively),
and calculated a 3D-separation of $\sim$43 kpc. Their team also found Peg III
to have an irregular shape elongated in the direction of Psc II.

More recently, Garofalo et al. (2021) used the Large Binocular Telescope (LBT)
to study variable stars in both UFDs. Using isodensity contour maps, they
found no support for a physical connection between Peg III and Psc II, as
neither UFD appeared to have an irregular shape. They suggested that the
regular structures of both UFDs eliminate the notion of a stellar stream or
another clear link between them. However, even if the two UFDs themselves have
not interacted, the possibility of them infalling as a pair or as part of a
group is not precluded (e.g., Wetzel et al., 2015).

Here, we present new, deep Hubble Space Telescope (HST) imaging of Peg III and
Psc II, allowing further exploration of how these two UFDs may or may not be
associated. We produce photometric catalogs and derive structural parameters,
integrated $V$-band magnitudes, and mass-to-light ratios. We also conduct an
orbital analysis of the two UFDs using proper motions (PMs) from Gaia (Gaia
Collaboration et al., 2016) early Data Release 3 (Gaia Collaboration et al.
2021; eDR3). In Section 2, we present our data and describe how they were
processed. We measure the structural parameters and calculate mass-to-light
ratios in Section 3. In Section 4, we use $Gaia$ eDR3 to examine whether Peg
III and Psc II could have had a past interaction. We discuss our results and
conclude in Sections 5 and 6, respectively.

##  2\. Hubble Space Telescope Data

###  2.1. Observations

The observations of Peg III and Psc II were performed using the F606W and
F814W filters of the HST Advanced Camera for Surveys (ACS) Wide Field Channel
(WFC) as part of Treasury program GO-14734 (PI: Kallivayalil). Parallel, off-
target fields were simultaneously taken with the Wide Field Camera 3 (WFC3)
using the same filters on the UV/visible (UVIS) channel. The Peg III
observations were taken on 2017 April 26 and 2017 May 2 using two orbits for
F606W and two orbits for F814W. Psc II was observed between 2017 June 19 and
2017 July 12 with two orbits dedicated to each filter. Each pair of ACS
exposures totaled 4744 s, while each WFC3 pair totaled 5117 s. The long
integration times allowed us to reach a signal-to-noise ratio (S/N) of 250 at
$m_{\rm{F606W}}{=}23$. The observations for both galaxies covered a single ACS
$202^{\prime\prime}\times 202\arcsec$ field and a single WFC3
$162^{\prime\prime}\times 162\arcsec$ off-field. The visits for each UFD were
performed within a restricted orientation range and over a four-point dither
pattern, optimized for the astrometric goals of this Treasury program, but
suitable for photometry as well.

###  2.2. Reduction and Photometry

Figure 1.— Color-magnitude diagrams of Pegasus III (left) and Pisces II
(right) in VEGAMAG. The typical color and magnitude errors are shown to the
left. Confirmed spectroscopic members in our catalog are shown in circles for
Peg III, along with a star that has ambiguous membership in a square (Kim et
al., 2016). There is one Peg III spectroscopic member that is not included in
our catalog due to it falling within the ACS chip gap. The Psc II
spectroscopic members with matches in our catalog are shown in circles (Kirby
et al., 2015). Only 3 of the 7 member stars fall within the ACS field-of-view
(FOV). The spectroscopic member with a Gaia proper motion (PM) is in the
diamond. The other stars used in the Psc II PM measurement were either too
bright or not in the ACS FOV. Only sources flagged as stars are shown. Figure
2.— Density contour maps of Peg III (left) and Psc II (right). Each map shows
10 levels defined by the color bar and in units of stars per $4\farcs 5\times
4\farcs 5$ bin. Note that the two color bars have different ranges and are not
integer values as they are from kernel density estimates. The white lines show
the borders of the ACS FOV, and the white arrows represent
$\mu_{\mathrm{tot}}$ over a timespan of 0.5 Myr from McConnachie & Venn
(2020). In this and all subsequent plots, the galaxies are oriented such that
North is in the direction of the increasing $y$-axis and East is in the
direction of the decreasing $x$-axis.

The images were processed and corrected for charge-transfer inefficiency (CTI)
using the current ACS and WFC3 pipelines. In each filter, the four dithered
images were combined using the DRIZZLE package (Fruchter & Hook, 2002) to
create drc fits files. Jackknife resampling was performed on the separate
dither images to create four three-dither-combined images for deriving
empirical errors. We used the photutils (Bradley et al., 2020) routines
DAOStarFinder and aperture_photometry to detect sources and calculate the flux
inside four- and six-pixel radii circular apertures.

We used the four- and six-pixel instrumental magnitudes ($m_{\mathrm{inst}}$)
to create a flag differentiating stars from galaxies. First, we calculated the
median magnitude difference between the four- and six-pixel radius values and
then determined the uncertainty on the difference between the four- and six-
pixel binned $m_{\mathrm{inst}}$ values as a function of the four-pixel
$m_{\mathrm{inst}}$. After accounting for an error floor, we fit a linear
relation to these uncertainties as a function of flux. We considered a source
to be a star if its magnitude difference from the median was within two sigma
of the fitted functional value. This criterion excluded galaxies, which showed
larger magnitude differences as well as deviation from the fitted relation
between the four- and six-pixel $m_{\mathrm{inst}}$. After applying the flags,
we accounted for the encircled energy corrections on the four-pixel radius
flux values, converted the flux to STMAG, and matched sources between the
F606W and F814W images. If a source was flagged as a star in either filter, it
was used in the analysis for this paper.

Sources in the three-dither combined images went through the same aperture
photometry pipeline and were matched across the four combinations in each
filter using a 6-parameter linear transformation. To derive empirical errors,
we took the standard deviation of the magnitudes of the sources found across
all of the combined images. The sources were matched in the F606W and F814W
filters, then matched to the drc source list using the same 6-parameter
transformation.

The observed magnitudes were corrected for dust extinction and reddening using
the dustmaps module (Green, 2018) with the Schlegel et al. (1998) maps and
Schlafly & Finkbeiner (2011) recalibration. We converted the ACS photometry
from STMAG filters to VEGAMAG $V$ and $I$ using the conversions and zero-
points given in Sirianni et al. (2005). To convert the WFC3 photometry, we
employed synphot (STScI Development Team, 2018) to generate representative
blackbodies of field stars (5200 K and 5100 K for the Peg III and Psc II off-
fields). Finally, we used stsynphot (STScI Development Team, 2020) to
calculate the filter conversion terms for these blackbodies and the WFC3
VEGAMAG zero-points for the observation dates. The resulting color-magnitude
diagrams (CMDs) and typical errors are shown in Figure 1, where cuts have been
made based on the “star-galaxy” flag.

To present the data and to explore the existence of the elongation in Peg III
measured by Kim et al. (2016), we created a contour map of stellar number
density by performing a kernel density estimate on the data using 67
4$\farcs$5 by 4$\farcs$5 bins. The left panel of Figure 2 shows the 10 levels
as filled contours for Peg III. The contour map for Psc II is illustrated in
the right panel of Figure 2, created using the same process as for the Peg III
map. This is consistent with density contour maps from Belokurov et al.
(2010), Sand et al. (2012), Muñoz et al. (2018), and Garofalo et al. (2021),
as all show Psc II with more regular levels. Comparing Psc II to Peg III, we
see that Psc II is more compact and has a higher peak smoothed surface
density.

##  3\. Structural Analysis

Figure 3.— Corner plot for the Pegasus III structural parameters from the
exponential fit. Positive correlations exist between the richness (number of
stars) and elliptical half-light radius as well as the elliptical half-light
radius and ellipticity. There is a negative correlation between the richness
and background density.

###  3.1. Fitting 2D-Profiles

To model the spatial structure of Peg III and Psc II, we followed the
technique described in Drlica-Wagner et al. (2020) and Simon et al. (2021),
largely based on the method shown in Martin et al. (2008). We modeled each UFD
with exponential and Plummer (1911) profiles and performed binned Poisson
maximum likelihood fits to the probability density functions with the
following free parameters: center position ($x_{0},y_{0}$), richness (number
of stars), 2D, projected semimajor axis of the ellipse that contains half of
the total integrated surface density of the galaxy (elliptical half-light
radius; $a_{h}$), ellipticity ($\epsilon$), position angle of the semimajor
axis measured East through North ($\theta$), and background surface density
($\Sigma_{b}$). In past literature on these UFDs, the elliptical half-light
radius $a_{h}$ has often been referred to as $r_{h}$, but here we have chosen
to make the explicit distinction between $a_{h}$ and the azimuthally-averaged
half-light radius $r_{h}$ (equal to $a_{h}\sqrt{1-\epsilon}$) for clarity and
the purpose of our comparison to simulations.

The normalized functional forms of the exponential and Plummer profiles are as
follows:

$\Sigma_{\mathrm{exp}}(r_{i})=\frac{1}{2\pi
r_{e}^{2}(1-\epsilon)}\textnormal{exp}\bigg{(}{-}\frac{r_{i}}{r_{e}}\bigg{)}$
(1)

$\Sigma_{p}(r_{i})=\frac{r_{p}^{2}}{2\pi(1-\epsilon)}(r_{i}^{2}+r_{p}^{2})^{-2}.$
(2)

Here, $r_{e}$ and $r_{p}$ are the scale lengths for each respective model,
with $r_{e}$=$1.68a_{h}$ and $r_{p}$=$a_{h}$, and the first term in each
equation is the normalization term, set to integrate to unity over all space.
The $r_{i}$ term is the elliptical radius of source $i$, defined as

$r_{i}=\bigg{\\{}\bigg{[}\frac{1}{1-\epsilon}(X_{i}\textnormal{cos}\theta-
Y_{i}\textnormal{sin}\theta)\bigg{]}^{2}+(X_{i}\textnormal{sin}\theta+Y_{i}\textnormal{cos}\theta)^{2}\bigg{\\}}^{1/2}.$
(3)

$X_{i}$ and $Y_{i}$ are the spatial offsets from the centroid, where
$X_{i}=x_{i}-x_{0}$ and $Y_{i}=y_{i}-y_{0}$.

To calculate the best-fitting parameters, we used the Markov chain Monte Carlo
(MCMC) ensemble sampler emcee (Foreman-Mackey et al., 2013). After adding the
background surface density term, we fit the following two functions:

$\Sigma_{\mathrm{exp,tot}}(r_{i})=\Sigma_{\mathrm{exp}}(r_{i})+\Sigma_{b}$ (4)

$\Sigma_{p\mathrm{,tot}}(r_{i})=\Sigma_{p}(r_{i})+\Sigma_{b}.$ (5)

We created $4\farcs 5$ by $4\farcs 5$ bins across the ACS field-of-view (FOV),
masking the area of the chip gap and outside the FOV, and counted the number
of stars in each bin. The corner plots (Foreman-Mackey, 2016) and posterior
distributions for the exponential fits are shown in Figures 3 and 4 for Peg
III and Psc II, respectively. For these and many subsequent figures, we choose
to show plots using the exponential fit values as there are not considerable
differences between the two profiles, and it allows easy comparison to what
has already been presented in the literature. The values from this work and
past works that measured structural parameters are shown in Tables 1 and 2,
with the center position transformed from pixel space to coordinates in right
ascension (RA) and declination (DEC). The spatially-binned, smoothed data, the
smoothed best-fit exponential model, and the residuals are shown in Figures 5
and 6.

Figure 4.— Same as Fig. 3, for Pisces II. Pegasus III
---
Parameter | Kim et al. (2015) | Kim et al. (2016) | This Work
$M_{V}$ | $-4.1\pm 0.5$ | $-3.4\pm 0.4$ |  $-4.17$ ${}^{+0.19}_{-0.22}$
Exponential
RA (h:m:s) | 22:24:22.6 $\pm 15\arcsec$ | 22:24:24.48 | 22:24:25.82 $\pm 5\arcsec$
DEC (d:m:s) | +05:25:12 $\pm 14\arcsec$ | +05:24:18.0 | +05:24:54.01 $\pm 3\arcsec$
$\theta_{\mathrm{exp}}$ (deg) |  $133\pm$ 17 | $114^{+19}_{-17}$ | $85\pm 8$
$\epsilon_{\mathrm{exp}}$ | $0.46^{+0.18}_{-0.27}$ | $0.38^{+0.22}_{-0.38}$ | $0.36^{+0.09}_{-0.10}$
$a_{h,\mathrm{exp}}$ (arcmin) | $1.3^{+0.5}_{-0.4}$ |  $0.85\pm$ 0.22 | $1.88^{+0.42}_{-0.33}$
$a_{h,\mathrm{exp}}$ (pc) | $78^{+30}_{-24}$ |  $53\pm$ 14 | $118^{+31}_{-30}$
Plummer
RA (h:m:s) | - | - | 22:24:25.78 $\pm 5\arcsec$
DEC (d:m:s) | - | - | +05:24:54.17 $\pm 3\arcsec$
$\theta_{p}$ (deg) | - | - | $83^{+8}_{-7}$
$\epsilon_{p}$ | - | - | $0.37^{+0.08}_{-0.09}$
$a_{h,p}$ (arcmin) | - | - | $1.67^{+0.26}_{-0.21}$
$a_{h,\mathrm{exp}}$ (pc) | - | - | $104^{+20}_{-23}$
Table 1Absolute magnitude and structural properties for Pegasus III, with the
top six lines reporting the absolute $V$-band magnitude and best-fit
exponential values and the bottom five listing the best-fit Plummer values.
The uncertainties reported for RA, DEC, $\theta$, $\epsilon$, and $a_{h}$ in
arcminutes are the 16th and 84th percentiles from the MCMCs. The uncertainties
on $M_{V}$ and $a_{h}$ in parsecs are the 16th and 84th percentiles from Monte
Carlo simulations that took into account the errors on the distance modulus.
Kim et al. (2015) and Kim et al. (2016) did not fit Plummer profiles to their
data, and Kim et al. (2016) did not provide uncertainties on their central
positions. Pisces II
---
Parameter | Belokurov et al. (2010) | Sand et al. (2012) | Muñoz et al. (2018) | This Work
$M_{V}$ | $-5.0\pm 0.5$ | $-4.1\pm 0.4$ | $-4.22\pm 0.38$ | $-4.28^{+0.19}_{-0.16}$
Exponential
RA (h:m:s) | - | 22:58:32.33 $\pm 5\arcsec$ | 22:58:32.28 $\pm 9\farcs 15$ | 22:58:32.76 $\pm 2\arcsec$
DEC (d:m:s) | - | +05:57:17.7 $\pm 4\arcsec$ | +05:57:09.36 $\pm 5\farcs 7$ | +05:57:20.36 $\pm 1\arcsec$
$\theta_{\mathrm{exp}}$ (deg) | - | 107 | $98\pm 13$ | $97\pm 3$
$\epsilon_{\mathrm{exp}}$ | - | <0.28 | $0.39\pm 0.10$ | $0.37\pm 0.04$
$a_{h,\mathrm{exp}}$ (arcmin) | - | $1.09\pm 0.19$ | $1.18\pm 0.20$ | $1.31^{+0.10}_{-0.09}$
$a_{h,\mathrm{exp}}$ (pc) | - | $58\pm 10$ | $62.5\pm 10.6$ | $69\pm 8$
Plummer
RA (h:m:s) | 22:58:31 $\pm 6\arcsec$ | 22:58:32.20 $\pm 5\arcsec$ | 22:58:32.28 $\pm 9\farcs 15$ | 22:58:32.75 $\pm 2\arcsec$
DEC (d:m:s) | +05:57:09 $\pm 4\arcsec$ | +05:57:16.3 $\pm 4\arcsec$ | +05:57:09.36 $\pm 5\farcs 7$ | +05:57:19.96 $\pm 1\arcsec$
$\theta_{p}$ (deg) | $77\pm 12$ | $110\pm 11$ | $78\pm 20$ | $98\pm 3$
$\epsilon_{p}$ | $0.4\pm 0.1$ | $0.33\pm 0.13$ | $0.34\pm 0.10$ | $0.37^{+0.03}_{-0.04}$
$a_{h,p}$ (arcmin) | $1.1\pm 0.1$ | $1.12\pm 0.18$ | $1.12\pm 0.16$ | $1.34^{+0.08}_{-0.07}$
$a_{h,p}$ (pc) |  $58\pm 5$* | $60\pm 10$ | $59.3\pm 8.5$ | $71\pm 8$
Table 2Same as Table 1, for Pisces II. Belokurov et al. (2010) did not fit an
exponential model to their data. Sand et al. (2012) had unconstrained
uncertainties for the position angle in their exponential model and their
reported ellipticity is the 68% upper confidence limit. *The uncertainty on
the $a_{h}$ in parsecs was derived for this work from a Monte Carlo simulation
using the uncertainty on the $a_{h}$ in arcminutes and the distance modulus of
21.3 reported in Belokurov et al. (2010).

The parameters of elliptical half-light radius and background surface density
appear both correlated with the richness (number of stars) and with each
other. The elliptical half-light radius has a positive correlation with the
number of stars, while the background density has a negative correlation with
the richness parameter. The elliptical half-light radius and background
density parameters also have a negative correlation. One would expect the
richness to increase with the elliptical half-light radius for a galaxy with
the same central surface density. If the elliptical half-light radius is
smaller, however, then a greater percentage of stars within the FOV would be
expected to belong to the background. Similarly, for a given number of stars
within a FOV, if a higher percentage of stars belongs to the galaxy, then a
lower number will be attributed to the background.

We show the best-fitting exponential and Plummer profiles for each galaxy
(left: Peg III, right: Psc II) in Figure 7. Surface density measurements taken
at 0.1 increments of the elliptical half-light radius are plotted against the
elliptical radius $R_{e}$ of each annular bin.

Figure 5.— Left: Peg III data, smoothed using a Gaussian kernel with a FWHM of
0$\farcs$5\. Middle: Best-fitting exponential model, smoothed with a
0$\farcs$5 FWHM Gaussian kernel. Right: The residuals between the smoothed
data and the model. Figure 6.— Left: Psc II data, smoothed using a Gaussian
kernel with a FWHM of 0$\farcs$5\. Middle: Best-fitting exponential model,
smoothed with a 0$\farcs$5 FWHM Gaussian kernel. Right: The residuals between
the smoothed data and the model. Figure 7.— Best-fitting models for the
surface density profiles of Pegasus III (left) and Pisces II (right) plotted
against the elliptical radius. The curves show the best-fit one dimensional
exponential (solid) and Plummer (dashed) profiles, which are fit to the full
stellar distribution, not the annularly binned data shown here. The gray
dashed line marks the approximate point where the elliptical annuli used for
the surface density measurements begin to cover area outside the FOV in pixel
space without source information in our catalogs. This was corrected for by
dividing the number of stars in each elliptical annulus by only the area
overlapping the FOV. The diamonds (squares) represent the surface density
measurements in bins using the exponential (Plummer) model, with elliptical
radii in increments of 0.1$a_{h}$. The errors come from a Poisson
distribution. The $a_{h}$ marks the data points corresponding to the best-fit
elliptical half-light radius.

While the structural parameters we measured for Psc II are consistent with
previous literature values, the $a_{h}$ value of Peg III
(1$\farcm$88${}^{+0.42}_{-0.33}$) is much larger ($\sim$2.2$\times$) than the
most recent literature value from Kim et al. (2016) ($0\farcm 85\pm 0.22$,
referred to as $r_{h}$ in their paper). This can be seen in the leftmost panel
of Figure 8. Here, we illustrate the best-fit 2D-exponential models projected
onto the sources in our FOV. The best-fit Plummer model from this work is
shown as the dashed ellipse. The same comparisons are shown for Psc II, with
the exponential fits shown in the middle panel and the Plummer fits in the
right panel.

Figure 8.— Comparison of past literature fits. Individual sources are colored
based on the binned surface densities of the best-fitting exponential models
in the left and center panels and the best-fitting Plummer model in the right
panel. Left: Pegasus III. Exponential profile fits from Kim et al. (2015) and
Kim et al. (2016) are shown in green and orange, respectively. The exponential
profile fit from this work is shown in solid blue. The Plummer fit from this
work is shown as the dashed blue ellipse. The central position is shown with
error bars for Kim et al. (2015) in green and for this work’s exponential fit
in blue. The central position from this work’s Plummer fit is not shown as it
almost completely overlaps with the exponential value. The central position
from Kim et al. (2016) is marked with an orange star, as no errors were
reported. Middle: Exponential profile fits and central positions with error
bars for Pisces II. Belokurov et al. (2010) is shown in fuchsia, Sand et al.
(2012) in green, Muñoz et al. (2018) in orange, and this work in blue. Right:
Plummer profile fits, with the same assigned colors as in the middle panel.

As discussed in Muñoz et al. (2012), the FOV should be at least $3\times$ the
half-light radius (assuming circular symmetry) in order to measure the
morphology to 10% accuracy. With the ACS FOV width of $202\arcsec$, or
$3\farcm 36$, we cover approximately 3.96$\times$ the $0\farcm 85$ value for
the elliptical half-light radius that we expected based on Kim et al. (2016).
However, when compared to our measured values for Peg III, the ACS FOV is
$\sim$1.8$\times$ our exponential $a_{h}$ and $\sim$2$\times$ our Plummer
$a_{h}$. For Psc II, the ratios are $\sim$2.6$\times$ and $\sim$2.5$\times$
for the exponential and Plummer $a_{h}$ values, respectively. We explore
trends in the structural parameters as a function of $a_{h}$ to FOV ratio in
the next subsection.

###  3.2. Mock Galaxy Tests

To test how well our fitting procedure could recover the elliptical half-light
radius for a galaxy given a limited FOV, we created simulated galaxies in
Cartesian space based on the best-fit values for Peg III. After generating
distributions based on the probability density functions of the Plummer and
exponential models, we randomly drew stars to match the length of the final
Peg III catalog and added a set number of background stars matching the best-
fit background density value. These simulated data were then masked with the
ACS FOV and chip gap and run through the MCMC. The results of these tests are
shown in Figure 9.

From this, we see that the code can recover $a_{h}$ up to $\sim$1$\farcm 4$
before it begins to return larger than the input $a_{h}$ values. Even up to
$1\farcm 8$ however, the expected $a_{h}$ values do fall within the
uncertainties for both models, and remain close up through $2\farcm 5$. Our
measured $a_{h}$ values for Peg III of $1\farcm 88$ and $1\farcm 67$ for the
exponential and Plummer profiles, respectively, could therefore be taken as
upper limits on the true value. The the maximum well-recovered value of
$\sim$1$\farcm 4$ extending beyond the previously reported $a_{h}=0\farcm 85$
(Kim et al., 2016) suggests that our modeling would have been able to capture
if Peg III was indeed so compact.

Figure 9.— Average difference of the MCMC-fit output $a_{h}$ and set input $a_{h}$ for 100 mock galaxies. The exponential fit differences are shown in cyan and the Plummer in orange. The black horizontal line at 0 represents the ideal case, where the fit $a_{h}$ matches the input $a_{h}$. The error bars represent the 16th and 84th percentiles from the 100 fits. The analysis for each mock galaxy is performed over a FOV equivalent to that of the HST ACS FOV, with the gray dashed vertical line marking where the input $a_{h}$ begins to be greater than one-half of the ACS FOV. Parameter | Peg III | Psc II
---|---|---
$(m-M)_{0}$ | $21.66\pm 0.12$ | $21.31\pm 0.17$
$v_{\mathrm{GSR}}$ (km s-1) | $-67.6\pm 2.6$ | $-79.9\pm 2.7$
$\sigma_{v}$ (km s-1) | $5.4^{+3.0}_{-2.5}$ | $5.4^{+3.6}_{-2.4}$
$\mu_{\alpha}$cos$\delta$ (mas yr-1) | $0.06\pm 0.1$ | $0.11\pm 0.11$
$\mu_{\delta}$ (mas yr-1) | $-0.2\pm 0.1$ | $-0.24^{+0.12}_{-0.11}$
$r_{h}$ (arcmin) | $1.51^{+0.35}_{-0.29}$ | $1.04\pm 0.08$
$r_{h}$ (pc) | $94^{+25}_{-24}$ | $55\pm 6$
$m_{V}$ | $17.50^{+0.15}_{-0.21}$ | $17.03\pm 0.04$
$M_{1/2}$ ($10^{6}$ M⊙) | $3.2^{+4.3}_{-2.1}$ | $1.9^{+3.3}_{-1.3}$
($M/L_{V}$)1/2 (M${}_{\odot}/$L⊙) | $1600^{+480}_{-580}$ | $850^{+570}_{-260}$
Table 3Adopted and derived values for Peg III and Psc II. The Peg III distance
modulus is from Kim et al. (2016), and the distance modulus for Psc II is from
Sand et al. (2012). The $\sigma_{v}$ values are from Kim et al. (2016) and
Kirby et al. (2015) for Peg III and Psc II, respectively, while the PMs are
both taken from McConnachie & Venn (2020). The azimuthally-averaged half-light
radius $r_{h}$ is calculated using the relation
$r_{h}{=}a_{h}\sqrt{1-\epsilon}$. ($M/L_{V}$)1/2 is the mass-to-light ratio
within the elliptical half-light radius. The processes for the derivation of
$m_{V}$, $M_{1/2}$, and ($M/L_{V}$)1/2 are described in Sections 3.3 and 3.4.
The uncertainties on all derived quantities are the 16th and 84th percentiles
of Monte Carlo simulations including the full-error space of all the relevant
terms.

###  3.3. Apparent and Absolute V-Band Magnitude Calculations

To derive the integrated magnitudes of each UFD, we opted for a probabilistic
model approach as opposed to using discrete stars. To begin, we selected a box
in color-magnitude space from $0\leq(V-I)\leq 1.2$ and $20\leq m_{V}\leq 29$.
For both the on-field and off-field, we created a Gaussian kernel using
scikit-learn (Pedregosa et al., 2011) Kernel Density and fit the kernel to the
stars inside the color-magnitude box (CMD box). We gridded this space into 100
bins along each dimension, resulting in color bins $\sim$0.01 dex and
magnitude bins $\sim$0.1 dex in width. We then generated random samples from
the on- and off-field kernels and computed the log-likelihood of each sample
under the model. The off-field was used to create a probabilistic background
model from which we estimated the excess flux. Integrating in both color and
magnitude space, we calculated the stellar density in the CMD box by
multiplying the log-likelihoods by the area of the CMD box and subtracting the
off-field model from the on-field. The integral returned the flux of the stars
inside the CMD box. We multiplied this flux by a correction factor derived
from the exponential models to account for the flux outside the FOV.
Converting this flux back to magnitude space yielded the integrated apparent
magnitude.

This calculation was performed within a Monte Carlo simulation (MC) that
included the individual source magnitude errors and FOV corrections calculated
from different sets of model parameters. We report the median integrated
$m_{V}$ values of $17.50^{+0.15}_{-0.21}$ and $17.03\pm 0.04$ for Peg III and
Psc II, respectively, with the uncertainties representing the 16th and 84th
percentiles. For the integrated $M_{V}$, we performed the same MC, this time
also including errors on the distance modulus. We used the Kim et al. (2016)
value of $21.66\pm 0.12$ for Peg III and the Sand et al. (2012) $21.31\pm
0.17$ value for Psc II. The most recent distance estimates (Garofalo et al.,
2021) indicate that the two UFDs have more similar heliocentric distances (Peg
III$\sim$174 kpc; Psc II$\sim$175 kpc), however, our CMDs in Figure 1 suggest
that Peg III is farther away than Psc II, as the horizontal branch and main-
sequence turnoff for Peg III are both $>$0.5 magnitudes less bright than those
of Psc II. As such, we chose to move forward with the previously reported
literature values in our MC, which gave the median
$M_{V}{=}-4.17^{+0.19}_{-0.22}$ for Peg III and $-4.28^{+0.19}_{-0.16}$ for
Psc II, with the 16th and 84th percentiles quoted as the uncertainties.

###  3.4. Mass-to-Light Ratios

Using the newly derived elliptical half-light radii values, we calculate
updated mass values using velocity dispersion values from Kirby et al. (2015)
and Kim et al. (2016). We use Equation 6 derived by Wolf et al. (2010), which
was also employed by Kim et al. (2016) to estimate the mass within the
elliptical half-light radius of Peg III. $R_{e}$ is defined as the
2D-projected half-light radius from elliptical fits of surface brightness
profiles and aligns with our $a_{h}$ values.

$M_{1/2}\simeq\frac{4}{G}\sigma^{2}_{v}R_{e}$ (6)

The Kim et al. (2016) $\sigma_{v}{=}5.4_{-2.5}^{+3.0}$ km s-1, determined
using seven member stars, and our elliptical half-light radius measurement of
118 parsecs at 215 kpc gives a $M_{1/2}{=}3.2^{+4.3}_{-2.1}\times 10^{6}$ M⊙
for Peg III. Converting our $M_{V}$ value to luminosity, we obtain the mass-
to-light ratio within the elliptical half-light radius
$(M/L_{V})_{1/2}{=}1600^{+480}_{-580}$ M${}_{\odot}/$L⊙. This is within one
sigma of the previously derived $1470^{+5660}_{-1240}$ M${}_{\odot}/$L⊙ from
Kim et al. (2016). Using the Kirby et al. (2015)
$\sigma_{v}{=}5.4^{+3.6}_{-2.4}$ km s-1, which was also calculated using seven
member stars, and our elliptical half-light radius of 69 parsecs at 183 kpc,
we calculate a $M_{1/2}{=}1.9^{+3.3}_{-1.3}\times 10^{6}$ M⊙ for Psc II. The
$(M/L_{V})_{1/2}$ is then $850^{+570}_{-260}$ M${}_{\odot}/$L⊙. This agrees
within two sigma with the value of $370^{+310}_{-240}$ M${}_{\odot}/$L⊙ that
Kirby et al. (2015) derived.

##  4\. Orbital Analysis

###  4.1. Previous Literature

Kim et al. (2015) and Kim et al. (2016) hypothesized a connection between Peg
III and Psc II immediately upon the discovery of the former, due at first to
their small on-sky separation (8.5∘), line-of-sight separation of $\sim$32
kpc, and similar heliocentric distances, and later from their similar radial
velocities and calculated physical separation of $\sim$36 kpc. Kim et al.
(2016) also noted the presence of an irregularity in the isodensity lines for
Peg III in the direction of Psc II, but included the caveat of small number
statistics.

Focusing on the ellipticity of Peg III, which we find to be
$0.36^{+0.09}_{-0.10}$, Kim et al. (2016) considered that, if not caused by
the formation process, this could be the result of tidal interactions with the
MW, a more likely scenario if Peg III is not in dynamic equilibrium or if
unresolved binaries are inflating the velocity dispersion. In such a case,
they suggested that Peg III might be a tidally disrupted dwarf galaxy remnant
but would need to have had a highly eccentric orbit to achieve close enough
Galactocentric distances to experience tidal effects. Alternatively, the
ellipticity of Peg III could have been caused by a tidal interaction with Psc
II, though at the time, there was not enough orbital information for Kim et
al. (2016) to further investigate this possibility.

###  4.2. Proper Motion Measurements

Recently, several authors have measured the proper motions (PMs) of Peg III
and Psc II using data from Gaia eDR3 (McConnachie & Venn, 2020; Li et al.,
2021; Battaglia et al., 2021), now making it possible to explicitly calculate
the orbital histories of both galaxies simultaneously to assess whether or not
they could be or have been a bound pair.

Li et al. (2021) reported a PM for Psc II using three stars and found a highly
eccentric orbital history such that Psc II only reaches its closest distance
relative to the MW at its current distance. By integrating both forward and
back for 10 Gyr, they also concluded that Psc II shows a high probability of
being unbound from the MW’s halo (95.85%).

Battaglia et al. (2021) independently measured the PM of Psc II using two
member stars from Gaia eDR3, and though they measured a different PM for Psc
II and included the infall of a massive LMC in their models, their
measurements still implied that Psc II is not and has never been bound to the
MW. Both Battaglia et al. (2021) and Li et al. (2021) tested a similar range
of MW mass models to arrive at these conclusions.

McConnachie & Venn (2020) reported the first PM value for Peg III measured
using only one member star, in addition to their own estimate for the PM of
Psc II using four spectroscopic member stars, three of which had radial
velocities more than 3$\sigma$ away from the mean velocity of the system. The
reported uncertainties were smaller than those of Li et al. (2021) and
Battaglia et al. (2021) by a factor of two or more. However, unlike Li et al.
(2021) and Battaglia et al. (2021), McConnachie & Venn (2020) imposed a prior
on their proper motion determinations such that the galaxies are assumed to be
bound to the MW (in direct contrast to the orbital histories suggested by the
Battaglia et al. (2021) and Li et al. (2021) PMs). Given that Peg III and Psc
II are outer MW dwarfs, such a restrictive prior biases the measured PMs
towards zero. In this work, we calculated orbital histories using the
McConnachie & Venn (2020) PMs (see Table 3) for both galaxies given the
uniform methodology that yields PMs with the lowest reported uncertainties for
both dwarfs, keeping in mind that the prior could affect the resulting orbital
histories and subsequent interpretations.

###  4.3. Orbital Model and Results

MW only potential
---
Orbit | $\rm f_{peri}$ |  $\rm t_{peri}$ [Gyr] |  $\rm r_{peri}$ [kpc] | $\rm f_{apo}$ |  $\rm t_{apo}$ [Gyr] |  $\rm r_{apo}$ [kpc]
MW1
PegIII-MW | 0.53 | 2.56 [2.61,4.53] | 20 [24,115] | 0.67 | 0.58 [0.64,2.61] | 217 [205,279]
PscII-MW | 0.58 | 2.56 [2.54,4.65] | 9 [22,103] | 0.70 | 0.71 [0.71,2.42] | 203 [196,270]
PegIII-PscII | 0.75 | 1.8 [0.06,3.65] | 12 [20,113] | 0.61 | 2.4 [2.21,4.63] | 17 [159,572]
MW2
PegIII-MW | 0.7 | 2.07 [2.14,4.12] | 15 [26,130] | 0.82 | 0.42 [0.46,2.24] | 213 [201,278]
PscII-MW | 0.73 | 2.06 [2.05,3.85] | 11 [22,112] | 0.84 | 0.5 [0.51,2.23] | 197 [189,269]
PegIII-PscII | 0.84 | 2.1 [0.07,3.47] | 7 [22,128] | 0.82 | 5.31 [1.8,4.05] | 256 [146,523]
Table 4Orbital parameters for Peg III and Psc II using the McConnachie & Venn
(2020) proper motions for both galaxies. The orbital parameters listed here
are for the representative cases shown in Figure 10, rather than the median
orbital history. The listed uncertainties correspond to the 16th and 84th
percentiles around the median of the distributions calculated using 1,000
Monte Carlo samples in the MW-only potential to illustrate the range of
orbital uncertainty. $\rm f_{peri}$ ($\rm f_{apo}$) is the fraction of 1,000
orbits where a pericenter (apocenter) is recovered, $\rm t_{peri}$ ($\rm
t_{apo}$) is the time at which the pericenter (apocenter) occurs on average,
and $\rm r_{peri}$ ($\rm r_{apo}$) is the distance of the most recent
pericenter (apocenter). Pericenter and apocenter are defined as the critical
minima and maxima along orbital trajectories in Galactocentric distance. MW-
LMC potential
---
Orbit | $\rm f_{peri}$ |  $\rm t_{peri}$ [Gyr] |  $\rm r_{peri}$ [kpc] | $\rm f_{apo}$ |  $\rm t_{apo}$ [Gyr] |  $\rm r_{apo}$ [kpc]
MW1
PegIII-MW | 0.71 | 1.88 [1.81,3.12] | 95 [29,145] | 0.76 | 0.3 [0.3,0.74] | 211 [194,239]
PegIII-LMC | 0.71 | 0.85 [1.81,3.12] | 32 [29,145] | 0.76 | 5.92 [0.3,0.74] | 817 [194,239]
PscII-MW | 0.72 | 2.08 [1.75,3.3] | 99 [31,125] | 0.81 | 0.35 [0.33,1.05] | 193 [182,225]
PscII-LMC | 0.72 | 0.81 [1.75,3.3] | 31 [31,125] | 0.81 | 4.8 [0.33,1.05] | 685 [182,225]
PegIII-PscII | 0.8 | 1.16 [0.06,2.7] | 10 [21,97] | 0.75 | 2.18 [1.4,4.05] | 42 [106,441]
MW2
PegIII-MW | 0.8 | 1.64 [1.59,2.64] | 55 [31,150] | 0.89 | 0.25 [0.25,0.82] | 210 [193,248]
PegIII-LMC | 0.8 | 0.85 [1.59,2.64] | 56 [31,150] | 0.89 | 3.35 [0.25,0.82] | 476 [193,248]
PscII-MW | 0.82 | 1.64 [1.5,2.79] | 66 [32,136] | 0.91 | 0.29 [0.28,1.22] | 191 [180,227]
PscII-LMC | 0.82 | 0.81 [1.5,2.79] | 49 [32,136] | 0.91 | 3.38 [0.28,1.22] | 464 [180,227]
PegIII-PscII | 0.88 | 1.16 [0.07,2.84] | 18 [21,113] | 0.9 | 1.7 [1.28,3.49] | 21 [109,440]
Table 5Same as Table 4 but for the MW-LMC potential. Here, the orbital
parameters of Peg III and Psc II calculated relative to the LMC’s orbit are
also included.

Here, we explore whether the two UFDs could have interacted with each other,
with the MW, or with the LMC. Radial velocities are adopted from Kim et al.
(2016; Peg III) and Kirby et al. (2015; Psc II) and distance moduli111During
the course of the analysis for this paper, new distance measurements were
published by Garofalo et al. (2021), derived using RR Lyrae. The range of
errors we considered in the distance moduli for the orbital analysis spans a
large portion of their values of $174\pm 18$ kpc for Peg III and $175\pm 11$
kpc for Psc II, thus we did not perform separate calculations here because the
PM error space is by far the dominant source of measurement uncertainties. are
adopted from Kim et al. (2016; Peg III) and Sand et al. (2012; Psc II),
respectively. Direct orbital histories are calculated using initial phase
space coordinates that equal the direct transformation of PMs, LOS velocity,
and distance converted to Galactocentric phase-space coordinates.

Following the methodology of Patel et al. (2020), we integrate orbits for Peg
III and Psc II in a three-component (Hernquist bulge + Miyamoto-Nagai disk +
NFW halo) MW-only potential and a combined three-component MW and two-
component (Miyamoto-Nagai disk + NFW halo) LMC potential using a fixed LMC
mass of $1.8\times 10^{11}$ M⊙. Here, the MW is allowed to move in response to
the LMC’s gravitational influence. Two MW masses are considered: MW1 has a
virial mass of $10^{12}$ M⊙, while MW2 has a virial mass of $1.5\times
10^{12}$ M⊙. Peg III and Psc II are modeled as Plummer spheres with a total
mass of $10^{9}$ M⊙ and a Plummer scale radius of 1 kpc (Jeon et al., 2017).
Dynamical friction owing to the MW and the LMC is also included, as in Patel
et al. (2020). See Patel et al. (2020) Tables 3 and 4 for all MW and LMC model
parameters.

One important difference with respect to the methodology of Patel et al.
(2020) is that we have explicitly accounted for the gravitational forces that
Peg III and Psc II exert on each other since it has been speculated that they
may be a bound pair. Thus, we calculate the joint, 4-body orbital history of
the MW, LMC (when applicable), Peg III, and Psc II. We have also checked
whether including the additional gravitational influence of the SMC is
significant and find that the SMC negligibly affects the orbits of Peg III and
Psc II, thus we do not incorporate it in the rest of this analysis.

The direct orbital histories of Peg III (cyan) and Psc II (magenta) using the
McConnachie & Venn (2020) PM values are shown in Figure 10. These orbits do
not represent the uncertainties on the measured quantities (i.e., distance,
proper motion, LOS velocity). Uncertainties on the orbital parameters that do
account for such are provided in Tables 4 and 5.

Figure 10.— Top: Direct orbital histories of Peg III (cyan) and Psc II
(magenta) relative to the MW, the LMC, and to each other are shown in the top,
middle, and bottom panels, respectively. Gaia eDR3 proper motions from
McConnachie & Venn (2020) are used for both Peg III and Psc II. MW1 (left) has
a virial mass of $10^{12}$ M⊙ while MW2 (right) has a virial mass of
$1.5\times 10^{12}$ M⊙. Dashed lines represent orbits computed in a MW-only
potential and solid lines indicate a MW+LMC potential with an LMC mass of
$1.8\times 10^{11}$ M⊙. The presence of the LMC decreases the distance at
pericenter relative to the MW-only model by $\sim$40–80 kpc. Middle: The
orbital history of each UFD relative to the LMC. For both MW1 and MW2, the
satellites complete a percentric passage about the LMC at $\sim$1 Gyr ago.
Bottom: The orbit of Peg III and Psc II relative to one another. In both MW
potentials, when the LMC is included, the two UFDs remain within $\sim$60 kpc
of each other in the last 4 Gyr, indicating that the satellites may be an
interacting pair that could have originated in a similar environment.

Orbital histories with respect to the MW (top panel, Figure 10) in the MW-only
potential (dashed lines) indicate that both Peg III and Psc II complete a very
close ($\leq$20 kpc), pericentric approach around the MW at $\sim$2.5 Gyr
(MW1) or $\sim$2 Gyr ago (MW2). However, when the influence of the LMC is
included, the satellites complete a pericentric approach around the MW at
$\sim$1.5 Gyr ago in both MW mass models at distances of $\sim$50–90 kpc.
Thus, the influence of the LMC significantly increases the distance at
pericenter about the MW and changes the timing at which the pericenter occurs.
This effect is two-fold as it is caused by the gravitational influence of the
LMC itself in addition to the LMC causing the MW’s center of mass to move in
response to the passage of the LMC. We also test LMC masses in the range
8–$25\times 10^{10}$ M⊙ and find no significant differences in the resulting
orbital histories.

These results are in contrast to the role the LMC plays on the orbit of Crater
2, for example, where including the LMC plunges Crater 2 further into the halo
of the MW, providing a possible explanation for its puffed up morphology
(Torrealba et al., 2016; Caldwell et al., 2017; Fritz et al., 2018; Fu et al.,
2019; Ji et al., 2021). Peg III and Psc II are much more compact than Crater
2, however, and are thus much less vulnerable to tidal effects. In upcoming
work, we will further disentangle which physical effects associated with the
passage of the LMC (i.e., the formation of the DM wake, the large-scale
density perturbations introduced to the MW DM halo, or LMC tidal debris) are
most significant for altering the orbital histories of low mass satellite
companions in the MW’s halo (Garavito-Carmargo et al, in prep.). Indeed, in
our density contour maps (Figure 2) and in the recent Garofalo et al. (2021)
paper, there appear to be no signs of tidal effects, and while our Peg III
elliptical half-light radius is larger than previous literature values, it is
still an order of magnitude smaller than that of Crater 2 ($\sim$118 pc
compared to $\sim$1100 pc; Torrealba et al. 2016).

The middle panels of Figure 10 show the orbits of Peg III and Psc II with
respect to the LMC as the LMC also falls into the halo of the MW. The orbits
of the two UFDs are very similar to each other and exhibit a pericentric
passage about the LMC at about 0.8 Gyr ago at $\sim$30 kpc for MW1 and
$\sim$60 kpc for MW2. This yields further evidence that the dynamics of the
UFDs have been strongly impacted by the gravitational influence of the LMC,
even if they are not in close proximity to the LMC at present (Peg III is 219
kpc from the LMC and Psc II is 195 kpc from the LMC). It is unclear whether
the satellites were ever bound to the LMC based on the calculations presented
here. Note that for MW1, the LMC’s median infall time is approximately 1.4 Gyr
ago and about 7 Gyr ago for MW2 (Patel et al., 2017). Despite this large
variation in infall time for the LMC, the UFDs behave fairly similarly over
the last $\sim$3 Gyr in both MW mass potentials.

In addition to computing the direct orbital histories, we also assess how much
the measurement uncertainties in PM, LOS velocity, and distance affect the
corresponding orbital histories for Peg III and Psc II. To this end, we
calculate 1,000 orbits for each orbital model: light MW-only, heavy MW-only,
light MW+LMC, heavy MW+LMC. These 1,000 orbits use initial conditions drawn in
a Monte Carlo fashion from the 1$\sigma$ uncertainty on PM, LOS velocity, and
distance for each dwarf galaxy (plus the LMC when included), thus the phase
space uncertainties for both galaxies are jointly sampled in the 3-body and
4-body calculations.

In Tables 4 and 5 we list the distance and timing of pericenter and apocenter
for the representative direct orbital histories in Figure 10. We use the
direct values because the mean and median obtained from the Monte Carlo
distribution can be biased towards larger values, especially for the apocenter
(see Fritz et al., 2018). The uncertainties provided correspond to the 16th
and 84th percentiles around the median of the distributions resulting from
1,000 orbit calculations spanning the uncertainty in phase space coordinates.
We report orbital parameters for each satellite with respect to the MW, LMC,
and each other for both MW mass models. Note that orbital parameters are only
listed for the fraction of orbits where a pericenter222Note that some studies
will adopt the minimum relative distance when no critical minimum is recovered
(e.g., Battaglia et al., 2021) however, we do not follow this methodology as
the minimum distance does not always correspond to a true pericentric passage.
(or apocenter) is recovered along the orbital trajectory. Therefore, the
uncertainties do not encompass the orbital histories where no such critical
minima and maxima are recovered in the last 6 Gyr.

As we do not account for the evolving MW potential or large scale structure,
we only integrate orbits for the last 6 Gyr. Furthermore, it has been shown
that the LMC significantly perturbs the shape and density of the MW’s dark
matter halo, loses a significant portion of its own mass, and creates a
trailing wake of dark matter over the last $\sim$2 Gyr (Garavito-Camargo et
al., 2019, 2021). Thus, the most physically accurate model of Peg III and Psc
II’s orbital histories would be calculated in a time-evolving MW+LMC potential
accounting for such features (Garavito-Carmargo et al., in prep.). Given the
significant PM uncertainties and the complex dynamics of the MW-LMC system, we
note that the orbits presented here are only an approximation of the most
plausible orbital histories within the measured phase space of these dwarf
galaxies, assuming they are indeed bound to the MW (as imposed by the prior in
McConnachie & Venn, 2020).

Figure 11.— Comparison of simulated UFDs to observed MW satellites. The
simulated UFDs (Jeon et al. 2017, blue diamonds; Jeon et al. 2021a, pink
diamonds; Jeon et al. 2021b, purple diamonds; Wheeler et al. 2019, orange
inverted triangles) are the same in both panels, except for the Applebaum et
al. (2021) (green triangles), which are split between the two panels based on
their simulated Galactocentric distances. Simulated UFDs can be distinguished
from observed satellites as the latter all have error bars. The observed data
are also split based on their Galactocentric distance, with satellites within
100 kpc shown in the left panel and satellites beyond 100 kpc in the right.
Confirmed MW UFDs are shown as filled squares, candidate MW UFDs as open
squares, MW globular clusters (GCs) as open circles, and ambiguous MW
satellites as x’s. The purple symbols use data from the updated McConnachie
(2012) table, while the gray symbols (and Leo T, in orange) are from Muñoz et
al. (2018). The dashed line represents a constant surface brightness of 32 mag
arsec-2, approximately the current observational limit. Our measured Peg III
(cyan star) and Psc II (magenta star) half-light radius and $M_{V}$ values lie
in the right panel, within the $r_{h}$ and $M_{V}$ range of other MW UFDs as
well as the Jeon et al. (2021a, b) simulated field UFDs. Some inner-halo
satellites from the left panel also fall within the Jeon et al. (2021a, b)
range.

##  5\. Discussion

###  5.1. Comparison to Simulations and Other UFDs

Peg III and Psc II are among the few known UFDs in the outer halo of the MW,
at Galactocentric distances of $\sim$213 and $\sim$182 kpc, respectively. To
see how our measured Peg III and Psc II azimuthally-averaged half-light radius
($r_{h}$) and $M_{V}$ values compare with those of other faint MW satellites,
we place them in the size-luminosity plane (Figure 11). In addition to
observed MW satellites (with $a_{h}$ converted to $r_{h}$ where necessary;
McConnachie 2012; Muñoz et al. 2018), we have also included values from five
sets of simulated UFDs: Jeon et al. (2017, 2021a, 2021b), Wheeler et al.
(2019), and Applebaum et al. (2021). The observed satellites are split between
two panels to explore the difference between those closer (left; $<$100 kpc)
and farther (right; $>$100 kpc) from the Galactic Center. The simulated
galaxies are shown in both panels, except for the Applebaum et al. (2021)
values, which are split between the two according to the Galactocentric
distances determined within their framework.

To derive $M_{V}$ values for the Jeon et al. and Wheeler et al. (2019)
simulations, we used Starburst99 (Leitherer et al., 1999) to convert from
stellar mass, based on the Padova evolutionary tracks (Bressan et al., 1993;
Fagotto et al., 1994a, b; Girardi et al., 2000), and the given half-stellar-
mass radii as (circular) half-light radii. The half-light radii that Applebaum
et al. (2021) report are circular and derived from the summation of particle
luminosities.

The Wheeler et al. (2019) simulations were unable to produce UFDs with half-
light radii lower than 200 pc, which the authors suggested could be in tension
with current observations because telescopes might only be sensitive to the
“bright” cores of diffuse and relatively massive objects. Similarly, the
Applebaum et al. (2021) simulations did not produce any galaxies in the UFD
magnitude range with smaller $r_{h}$ than $\sim$300 pc, which may have been
due to their force softening. Among the simulations included in our
comparison, the Wheeler et al. (2019) and Jeon et al. (2021a, b) have the
highest resolution, with Wheeler et al. using $m_{\mathrm{gas}}{=}30$ and
$m_{gas}{=}250$ M⊙ in their high- and median-resolution simulations,
respectively, and Jeon et al. using $m_{\mathrm{gas}}{\sim}60$ M⊙. As
simulation techniques improve and are able to resolve a broader diversity of
UFDs in a MW-environment, it will be interesting to see whether they more
closely reproduce the scatter shown by McConnachie (2012) and Muñoz et al.
(2018).

For Figure 11, we have imposed an upper magnitude limit of $M_{V}{=}{-}7.7$ to
only include simulated galaxies in the UFD range, as Jeon et al. (2017),
Wheeler et al. (2019), and Applebaum et al. (2021) produced galaxies in the
dSph regime as well. As seen in Figure 11, Peg III and Psc II are well within
the observed $M_{V}$ and $r_{h}$ ranges for other MW UFDs. They also fall in
the range covered by the Jeon et al. (2021a, b) simulated field UFDs. Some
closer MW UFDs and MW UFD candidates also fall within the simulated field UFD
area, suggesting against this agreement being unique to outer-halo UFDs.
Additionally, there is a higher number of observed satellites with smaller
$r_{h}$ and fainter $M_{v}$ within 100 kpc, likely due to observational
constraints.

As new observatories come online with deeper detection limits, smaller and
fainter satellites beyond 100 kpc could be discovered that have no analogs in
current simulations. One might expect outer-halo satellites that have never
been within 100 kpc of the MW to be more compact than UFDs with similar masses
that have been closer and possibly subjected to strong tidal forces. If there
are no or very few distant UFDs found in this smaller and fainter regime, the
impact of the LMC and how it might have drawn UFDs closer ($\sim$50 kpc) at
some point in their orbital history could be considered. While our orbital
histories for Peg III and Psc II show the LMC increasing the pericenter
distance, it is possible that other MW satellites (e.g., Crater 2, Ji et al.
2021) could have experienced very different effects.

###  5.2. A Bound Pair?

Several proposed pairs of galaxies have been reported within the Local Group
(Pawlowski, 2021). Speculation that Peg III and Psc II are associated began in
the Peg III discovery paper (Kim et al., 2015), due to their on-sky proximity
and spatial separation of $\sim$30 kpc. The support for this association
increased with the Kim et al. (2016) measurement of a radial velocity for Peg
III within $\sim$10 km s-1 of that of Psc II. Additionally, they found the
same irregularity in structure in their deeper imaging as in the discovery
paper.

The next work to investigate the possible connection between Peg III and Psc
II was Garofalo et al. (2021). Using LBT data, they found no irregular shape
and concluded that there were no signs of past tidal interactions. In our
density contour maps (Figure 2), we find no signs of an irregular shape for
Peg III either. While previously limited to only looking at the morphologies
for hints of past interactions, Gaia eDR3 has provided us with the opportunity
to use kinematics to study the UFDs’ possible shared history. With proper
motions in hand, Peg III and Psc II are amongst the first pairs of UFDs to be
investigated further in this context.

To calculate a rough estimate of the likelihood that two satellite galaxies at
similarly large distances and with radial velocities comparable to Peg III and
Psc II would appear together by chance (i.e., not associated with each other
prior to entering the MW halo) with a similar relative distance and radial
velocity difference between them, we drew 1 million UFDs from a Gaussian
velocity distribution centered around 0 km s-1 with $\sigma{=}100$ km s-1 and
a distance distribution of $r^{3-1-\gamma}$ from 30 to 300 kpc, with
$\gamma{=}2.11$ taken from Fritz et al. (2020). We compared the drawn distance
and velocity values to those of Psc II, calculating the percentage of these
UFDs with $\Delta v\leq 15$ km s-1 and $\Delta r\leq 40$ kpc with respect to
the properties of Psc II (${\sim}{-}80$ km s-1; $\sim$183 kpc). Of the 1
million draws, 2.5% met the criteria, suggesting a low probability that two
galaxies with such similar properties to each other as Peg III and Psc II
share would exist in the MW halo by chance.

Furthermore, we can use the orbital histories computed using the McConnachie &
Venn (2020) PMs to speculate about whether Peg III and Psc II have a shared
orbital history or one that indicates that they have been bound in the recent
past. The bottom panel of Figure 10 shows the orbits of Peg III and Psc II
with respect to one another in the MW-only (dashed lines) and MW-LMC potential
(solid lines). In the MW-only potential, the satellites reach as close 20 kpc
to each other in the last 2–3 Gyr. For the heavy MW-only potential, Peg III
and Psc II remain significantly far from each other until $\sim$2 Gyr ago,
when they make a close ($<$20 kpc), pericentric passage about both the MW and
then subsequently each other. However, for the light MW-only potential, the
UFDs have remained within 80 kpc of each other for at least the last few
billion years.

When the influence of the LMC is included, an opposite trend is observed. For
the light MW model, the UFDs first interact with one another about 3 Gyr ago
(prior to the infall of the LMC) and then pass about each other again 1 Gyr
ago (after the infall of the LMC). They are otherwise significantly separated.
In the heavy MW+LMC model, the UFDs are within about 80 kpc of each other for
the last 6 Gyr. They pass around each other twice at a distance of $\sim$20
kpc at $\sim$1 and $\sim$2 Gyr ago. These results, summarized in Tables 4 and
5, illustrate that Peg III and Psc II may have had an intricate orbital
history that includes close interactions with the MW, LMC, and each other when
the McConnachie & Venn (2020) values are considered. These conclusions are
based on the direct orbital histories, yet the results in Table 4 and 5
illustrate that there is significant statistical uncertainty on the properties
of this system due to large measurement uncertainties on distance and proper
motion.

We also calculate a simple metric to test whether Peg III and Psc II are
currently bound to each other as was done in Geha et al. (2015) and Sohn et
al. (2020) for NGC 147 and NGC 185, two dwarf elliptical galaxies located
about 1∘ from each other. It is known that for two point masses to be
gravitationally bound, the potential energy of the system must be greater than
the kinetic energy. This yields the criterion $b\equiv 2GM_{sys}/\Delta
r\Delta v^{2}$. Thus, when $b>1$, the system is considered bound.

Assuming a relative Galactocentric distance of $\Delta r{=}36.4$ kpc between
the dwarfs, a relative Galactocentric velocity of $\Delta v{=}20.6$ km s-1 at
present day, and a mass of $10^{9}$ M⊙ for each dwarf, we find $b{=}1.11\pm
0.31$. This $b$ corresponds to the initial conditions used to calculate the
direct orbital histories and the uncertainty represents the standard deviation
across $b$ calculated for 1,000 initial conditions encompassing the measured
phase space of these two dwarf galaxies. Thus, Peg III and Psc II are
consistent with being bound to each other, as illustrated by the orbital
histories in the bottom panels of Figure 10. However, the uncertainty on $b$
indicates that some unbound orbits are possible within the measured phase
space. This, in combination with the non-zero transverse motions of these
dwarfs, indicates that improved distances and PMs, in addition to a more
precise understanding of the MW’s mass, are therefore necessary to
definitively conclude whether Peg III and Psc II have experienced a shared
orbital history over the last few billion years.

If Peg III and Psc II are indeed bound to each other, as suggested by the most
common orbits recovered from the McConnachie & Venn (2020) PMs, this would add
to the rare findings of confirmed pairs of satellites within the Local Group.
It has also been shown that pairs and groups of satellites that fall into the
halo of their hosts together tend to merge with one another, sometimes only
1–3 Gyr after infall (e.g., Deason et al., 2014; Wetzel et al., 2015). Peg III
and Psc II could thus be in the process of or on the way to the early phases
of merging. Conversely, another possibility would be for the two UFDs to
disperse in phase space soon after falling into the MW, depending on how
strongly the two were bound during infall (e.g., Deason et al., 2015).
Improved PMs will continue to shed light on the potential future trajectories
of these UFDs.

###  5.3. Current and Future Precision

To formally investigate the precision needed for determining an association,
we perform a Monte Carlo simulation, drawing the transverse velocities in both
dimensions (X and Y) for two galaxies in both an unassociated and associated
case. To do this, we must consider the dispersion of the halo ($\sim$100 km
s-1; Correa Magnus & Vasiliev 2021) as well as the proper motion
uncertainties. Gaia’s current measurement errors for the two UFDs are clearly
larger than the halo dispersion: without using a prior (which artificially
constrains the error), the best measurement determined for Psc II using Gaia
eDR3 has a precision of about 230 km s-1 (Li et al., 2021), while a proper
motion uncertainty cannot even be obtained for Peg III. Nonetheless, to
calculate an estimate of the necessary precision for determining the presence
of an association, we optimistically scale the Li et al. (2021) proper motion
error by the Kim et al. (2016) distance to obtain 270 km s-1 for Peg III.

In the unassociated case, transverse velocities for two galaxies are drawn
from a Gaussian distribution centered on 0 km s-1, with $1\sigma$ equal to the
halo dispersion of 100 km s-1 plus the measurement errors of each galaxy, as
described in the previous paragraph.

In the associated case, we again center the Gaussian distribution on 0 km s-1,
but instead set $1\sigma$ equal to only the current Gaia eDR3 measurement
errors of Peg III and Psc II (i.e., without including the halo dispersion). We
then calculate the velocity difference for each realization:

$V_{\mathrm{diff}}{=}\sqrt{(V_{\mathrm{Psc\,II,X}}-V_{\mathrm{Peg\,III,X}})^{2}+(V_{\mathrm{Psc\,II,Y}}-V_{\mathrm{Peg\,III,Y}})^{2}}.$
(7)

We take the median $V_{\mathrm{diff}}$ from the associated case and compare it
to the $V_{\mathrm{diff}}$ values calculated in the unassociated case and
determine how often the unassociated case produces a velocity difference less
than or equal to the median associated $V_{\mathrm{diff}}$. We obtain a median
$V_{\mathrm{diff}}$=418 km s-1 from the associated realizations. A velocity
difference of this size or smaller occurs with a 45% probability in the
unassociated case. Thus, it is not expected that Gaia eDR3 measurements
without a prior could give any constraints on whether or not Peg III and Psc
II are associated.

Now we investigate how future measurements could improve upon current
information to constrain an association between Peg III and Psc II. Gaia will
operate much longer, optimistically achieving a baseline of 11 years. Since
proper motion errors of continuously observing telescopes scale with time to
the power of $-1.5$ (e.g., Lindegren et al. 2021), the precision achievable
with the final Gaia data release would be approximately 36 and 30 km s-1 for
Peg III and Psc II, respectively. This precision could already be met and
surpassed using HST now. We estimate the achievable precision using the data
described here as a first epoch and with a second epoch taken as outlined in
Kallivayalil et al. (2015). A new measurement with HST imaging taken in 2023
(a six-year baseline) would measure the transverse velocity to a precision in
both dimensions of 30.3 km s-1 for Peg III and 24.5 km s-1 for Psc II.

Repeating the Monte Carlo simulation with the estimated HST errors, we obtain
a median $V_{\mathrm{diff}}$=46 km s-1 from the associated cases. A velocity
difference of this size or smaller occurs with a 4.7% probability in the
unassociated case. If a velocity difference less than $\sim$50 km s-1 was
observed, the small probability from the unassociated Monte Carlo cases,
combined with the already small probability of similar line-of-sight
velocities and distances (2.5%, see Section 5.2), would almost completely rule
out a chance association.

Given a longer baseline, HST(-like) observations taken with HST, JWST, or the
Roman Space Telescope could measure the velocity well enough to prove that the
two systems are associated with each other. For example, observations taken in
2034 would give a full transverse velocity error of 19.4 km s-1, smaller than
the current Galactocentric velocity difference ($\sim$20.6 km s-1).

##  6\. Conclusions

Based on deep HST imaging and Gaia eDR3 PMs, we have measured the structural
parameters and performed an orbital analyses for two distant MW UFD
satellites, Peg III and Psc II. For the structural analysis, we used a
2D-maximum likelihood MCMC to fit exponential and Plummer profiles, finding
central positions, elliptical half-light radii, ellipticities, position
angles, number of stars, and surface background densities.

For Peg III, we measured an elliptical half-light radius of
$1.88^{+0.42}_{-0.33}$ arcminutes, a position angle of $85\pm 8$ degrees, and
an ellipticity of $0.36^{+0.09}_{-0.10}$. The best-fit exponential model for
Psc II gave an elliptical half-light radius of $1.31^{+0.10}_{-0.09}$
arcminutes, a position angle of $97\pm 3$ degrees, and an ellipticity of
$0.37\pm 0.04$. Our Psc II measurements are within $1\sigma$ agreement when
compared to previous literature values, while we find a larger elliptical
half-light radius value for Peg III. See Tables 1 and 2 for full comparisons.
Future imaging of similar depth to ours paired with a larger FOV would help
clarify the true size of Peg III.

We compared the measured sizes and magnitudes of Peg III and Psc II to both
those of other observed faint MW satellites and simulations of isolated field
UFDs and satellite UFDs in a MW-like environment. While Peg III and Psc II are
more distant MW UFDs with $r_{h}$ and $M_{V}$ values that correspond well to
simulations of isolated field UFDs, they are not significantly distinct in
structure from other observed UFDs in the inner MW halo. This could point to
Peg III and Psc II (and other outer-halo UFDs) having been subject to tidal
forces from the MW and/or the LMC throughout their lifetimes.

The first orbital analysis, using solely the McConnachie & Venn (2020) PM
values, demonstrated the importance of including the LMC in such calculations,
as it affects the timing and pericenter of the orbits of both UFDs. In this
analysis, a statistically significant fraction of the computed orbital
histories in the measured phase space show Peg III and Psc II are consistent
with being a gravitationally bound pair today.

Comparing the vastly different results from the Li et al. (2021) study versus
our analysis using the McConnachie & Venn (2020) Psc II PM values and
considering that the Peg III PM measurement came from a single star, it is
also clear that improved measurements based on more member stars and a clearer
knowledge of the MW’s mass are necessary to narrow down the possible orbital
histories and allow a more definitive statement on whether Peg III and Psc II
are indeed a bound pair. These improved measurements could be taken from
imaging by HST now or measured from the final Gaia data release in years to
come. With the precision achievable from these baselines, measuring a velocity
difference between Peg III and Psc II of less than $\sim$50 km s-1 would
suggest with about 95% confidence that these UFDs are bound. This would
finally resolve the question and confirm if these two fascinating galaxies are
indeed bound, thereby providing a critical data point for understanding dwarf
galaxy evolution.

## Acknowledgements

These data are associated with the HST Treasury Program 14734 (PI:
Kallivayalil). Support for this program was provided by NASA through grants
from the Space Telescope Science Institute. This material is based upon work
supported by the National Science Foundation under grant No. AST-1847909. HR
acknowledges support from the Virginia Space Grant Consortium Graduate
Research STEM Fellowship. EP acknowledges support from the Miller Institute
for Basic Research in Science at University of California, Berkeley. This
research has made use of NASA’s Astrophysics Data System. Facility: HST
(ACS,WFC3) Software: Astrodrizzle (Fruchter & Hook, 2002); Astropy (Astropy
Collaboration et al., 2013, 2018); corner.py (Foreman-Mackey, 2016); dustmaps
(Green, 2018); emcee (Foreman-Mackey et al., 2013); Jupyter Notebook (Kluyver
et al., 2016); Matplotlib (Hunter, 2007); Numpy (Harris et al., 2020);
photutils (Bradley et al., 2020); scikit-learn (Pedregosa et al., 2011) Scipy
(Virtanen et al., 2020); stsynphot (STScI Development Team, 2020); synphot
(STScI Development Team, 2018)

## References

  * Ahn et al. (2014) Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2014, ApJS, 211, 17, doi: 10.1088/0067-0049/211/2/17
  * Applebaum et al. (2021) Applebaum, E., Brooks, A. M., Christensen, C. R., et al. 2021, ApJ, 906, 96, doi: 10.3847/1538-4357/abcafa
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Battaglia et al. (2021) Battaglia, G., Taibi, S., Thomas, G. F., & Fritz, T. K. 2021, arXiv e-prints, arXiv:2106.08819.  https://arxiv.org/abs/2106.08819
  * Belokurov et al. (2010) Belokurov, V., Walker, M. G., Evans, N. W., et al. 2010, ApJ, 712, L103, doi: 10.1088/2041-8205/712/1/L103
  * Bradley et al. (2020) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: 1.0.0, 1.0.0, Zenodo, doi: 10.5281/zenodo.4044744
  * Breddels & Helmi (2013) Breddels, M. A., & Helmi, A. 2013, A&A, 558, A35, doi: 10.1051/0004-6361/201321606
  * Bressan et al. (1993) Bressan, A., Fagotto, F., Bertelli, G., & Chiosi, C. 1993, A&AS, 100, 647
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343, doi: 10.1146/annurev-astro-091916-055313
  * Caldwell et al. (2017) Caldwell, N., Walker, M. G., Mateo, M., et al. 2017, ApJ, 839, 20, doi: 10.3847/1538-4357/aa688e
  * Correa Magnus & Vasiliev (2021) Correa Magnus, L., & Vasiliev, E. 2021, MNRAS, doi: 10.1093/mnras/stab3726
  * Deason et al. (2014) Deason, A., Wetzel, A., & Garrison-Kimmel, S. 2014, ApJ, 794, 115, doi: 10.1088/0004-637X/794/2/115
  * Deason et al. (2015) Deason, A. J., Wetzel, A. R., Garrison-Kimmel, S., & Belokurov, V. 2015, MNRAS, 453, 3568, doi: 10.1093/mnras/stv1939
  * Drlica-Wagner et al. (2020) Drlica-Wagner, A., Bechtol, K., Mau, S., et al. 2020, ApJ, 893, 47, doi: 10.3847/1538-4357/ab7eb9
  * Errani et al. (2018) Errani, R., Peñarrubia, J., & Walker, M. G. 2018, MNRAS, 481, 5073, doi: 10.1093/mnras/sty2505
  * Fagotto et al. (1994a) Fagotto, F., Bressan, A., Bertelli, G., & Chiosi, C. 1994a, A&AS, 104, 365
  * Fagotto et al. (1994b) —. 1994b, A&AS, 105, 29
  * Flores & Primack (1994) Flores, R. A., & Primack, J. R. 1994, ApJ, 427, L1, doi: 10.1086/187350
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fritz et al. (2018) Fritz, T. K., Battaglia, G., Pawlowski, M. S., et al. 2018, A&A, 619, A103, doi: 10.1051/0004-6361/201833343
  * Fritz et al. (2020) Fritz, T. K., Di Cintio, A., Battaglia, G., Brook, C., & Taibi, S. 2020, MNRAS, 494, 5178, doi: 10.1093/mnras/staa1040
  * Fruchter & Hook (2002) Fruchter, A. S., & Hook, R. N. 2002, PASP, 114, 144, doi: 10.1086/338393
  * Fu et al. (2019) Fu, S. W., Simon, J. D., & Alarcón Jara, A. G. 2019, ApJ, 883, 11, doi: 10.3847/1538-4357/ab3658
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Garavito-Camargo et al. (2019) Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2019, ApJ, 884, 51, doi: 10.3847/1538-4357/ab32eb
  * Garavito-Camargo et al. (2021) —. 2021, ApJ, 919, 109, doi: 10.3847/1538-4357/ac0b44
  * Garofalo et al. (2021) Garofalo, A., Tantalo, M., Cusano, F., et al. 2021, arXiv e-prints, arXiv:2105.11488.  https://arxiv.org/abs/2105.11488
  * Geha et al. (2015) Geha, M., Weisz, D., Grocholski, A., et al. 2015, ApJ, 811, 114, doi: 10.1088/0004-637X/811/2/114
  * Girardi et al. (2000) Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, A&AS, 141, 371, doi: 10.1051/aas:2000126
  * Goerdt et al. (2006) Goerdt, T., Moore, B., Read, J. I., Stadel, J., & Zemp, M. 2006, MNRAS, 368, 1073, doi: 10.1111/j.1365-2966.2006.10182.x
  * Grand et al. (2021) Grand, R. J. J., Marinacci, F., Pakmor, R., et al. 2021, MNRAS, 507, 4953, doi: 10.1093/mnras/stab2492
  * Green (2018) Green, G. 2018, The Journal of Open Source Software, 3, 695, doi: 10.21105/joss.00695
  * Guerra et al. (2021) Guerra, J., Geha, M., & Strigari, L. E. 2021, arXiv e-prints, arXiv:2112.05166.  https://arxiv.org/abs/2112.05166
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jardel & Gebhardt (2013) Jardel, J. R., & Gebhardt, K. 2013, ApJ, 775, L30, doi: 10.1088/2041-8205/775/1/L30
  * Jeon et al. (2017) Jeon, M., Besla, G., & Bromm, V. 2017, ApJ, 848, 85, doi: 10.3847/1538-4357/aa8c80
  * Jeon et al. (2021a) —. 2021a, MNRAS, 506, 1850, doi: 10.1093/mnras/stab1771
  * Jeon et al. (2021b) Jeon, M., Bromm, V., Besla, G., Yoon, J., & Choi, Y. 2021b, MNRAS, 502, 1, doi: 10.1093/mnras/staa4017
  * Ji et al. (2021) Ji, A. P., Koposov, S. E., Li, T. S., et al. 2021, ApJ, 921, 32, doi: 10.3847/1538-4357/ac1869
  * Kallivayalil et al. (2015) Kallivayalil, N., Wetzel, A. R., Simon, J. D., et al. 2015, arXiv e-prints, arXiv:1503.01785.  https://arxiv.org/abs/1503.01785
  * Kim et al. (2015) Kim, D., Jerjen, H., Mackey, D., Da Costa, G. S., & Milone, A. P. 2015, ApJ, 804, L44, doi: 10.1088/2041-8205/804/2/L44
  * Kim et al. (2016) Kim, D., Jerjen, H., Geha, M., et al. 2016, ApJ, 833, 16, doi: 10.3847/0004-637X/833/1/16
  * Kirby et al. (2015) Kirby, E. N., Simon, J. D., & Cohen, J. G. 2015, ApJ, 810, 56, doi: 10.1088/0004-637X/810/1/56
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Scmidt (IOS Press), 87–90.  https://eprints.soton.ac.uk/403913/
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3, doi: 10.1086/313233
  * Li et al. (2021) Li, H., Hammer, F., Babusiaux, C., et al. 2021, ApJ, 916, 8, doi: 10.3847/1538-4357/ac0436
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2, doi: 10.1051/0004-6361/202039709
  * Martin et al. (2008) Martin, N. F., de Jong, J. T. A., & Rix, H.-W. 2008, ApJ, 684, 1075, doi: 10.1086/590336
  * McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4, doi: 10.1088/0004-6256/144/1/4
  * McConnachie & Venn (2020) McConnachie, A. W., & Venn, K. A. 2020, Research Notes of the American Astronomical Society, 4, 229, doi: 10.3847/2515-5172/abd18b
  * Moore (1994) Moore, B. 1994, Nature, 370, 629, doi: 10.1038/370629a0
  * Muñoz et al. (2018) Muñoz, R. R., Côté, P., Santana, F. A., et al. 2018, ApJ, 860, 66, doi: 10.3847/1538-4357/aac16b
  * Muñoz et al. (2012) Muñoz, R. R., Padmanabhan, N., & Geha, M. 2012, ApJ, 745, 127, doi: 10.1088/0004-637X/745/2/127
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Patel et al. (2017) Patel, E., Besla, G., & Sohn, S. T. 2017, MNRAS, 464, 3825, doi: 10.1093/mnras/stw2616
  * Patel et al. (2020) Patel, E., Kallivayalil, N., Garavito-Camargo, N., et al. 2020, ApJ, 893, 121, doi: 10.3847/1538-4357/ab7b75
  * Pawlowski (2021) Pawlowski, M. S. 2021, arXiv e-prints, arXiv:2109.02654.  https://arxiv.org/abs/2109.02654
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Plummer (1911) Plummer, H. C. 1911, MNRAS, 71, 460, doi: 10.1093/mnras/71.5.460
  * Read & Gilmore (2005) Read, J. I., & Gilmore, G. 2005, MNRAS, 356, 107, doi: 10.1111/j.1365-2966.2004.08424.x
  * Read et al. (2018) Read, J. I., Walker, M. G., & Steger, P. 2018, MNRAS, 481, 860, doi: 10.1093/mnras/sty2286
  * Read et al. (2019) —. 2019, MNRAS, 484, 1401, doi: 10.1093/mnras/sty3404
  * Read et al. (2021) Read, J. I., Mamon, G. A., Vasiliev, E., et al. 2021, MNRAS, 501, 978, doi: 10.1093/mnras/staa3663
  * Sand et al. (2012) Sand, D. J., Strader, J., Willman, B., et al. 2012, ApJ, 756, 79, doi: 10.1088/0004-637X/756/1/79
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Simon (2019) Simon, J. D. 2019, ARA&A, 57, 375, doi: 10.1146/annurev-astro-091918-104453
  * Simon et al. (2021) Simon, J. D., Brown, T. M., Drlica-Wagner, A., et al. 2021, ApJ, 908, 18, doi: 10.3847/1538-4357/abd31b
  * Sirianni et al. (2005) Sirianni, M., Jee, M. J., Benítez, N., et al. 2005, PASP, 117, 1049, doi: 10.1086/444553
  * Sohn et al. (2020) Sohn, S. T., Patel, E., Fardal, M. A., et al. 2020, ApJ, 901, 43, doi: 10.3847/1538-4357/abaf49
  * Strigari et al. (2007) Strigari, L. E., Bullock, J. S., Kaplinghat, M., et al. 2007, ApJ, 669, 676, doi: 10.1086/521914
  * Strigari et al. (2008) —. 2008, Nature, 454, 1096, doi: 10.1038/nature07222
  * Strigari et al. (2017) Strigari, L. E., Frenk, C. S., & White, S. D. M. 2017, ApJ, 838, 123, doi: 10.3847/1538-4357/aa5c8e
  * STScI Development Team (2018) STScI Development Team. 2018, synphot: Synthetic photometry using Astropy.  http://ascl.net/1811.001
  * STScI Development Team (2020) —. 2020, stsynphot: synphot for HST and JWST.  http://ascl.net/2010.003
  * Torrealba et al. (2016) Torrealba, G., Koposov, S. E., Belokurov, V., & Irwin, M. 2016, MNRAS, 459, 2370, doi: 10.1093/mnras/stw733
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wetzel et al. (2015) Wetzel, A. R., Deason, A. J., & Garrison-Kimmel, S. 2015, ApJ, 807, 49, doi: 10.1088/0004-637X/807/1/49
  * Wheeler et al. (2019) Wheeler, C., Hopkins, P. F., Pace, A. B., et al. 2019, MNRAS, 490, 4447, doi: 10.1093/mnras/stz2887
  * Willman et al. (2005) Willman, B., Dalcanton, J. J., Martinez-Delgado, D., et al. 2005, ApJ, 626, L85, doi: 10.1086/431760
  * Wolf et al. (2010) Wolf, J., Martinez, G. D., Bullock, J. S., et al. 2010, MNRAS, 406, 1220, doi: 10.1111/j.1365-2966.2010.16753.x
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513