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

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 for the purpose of learning more about the UFD stars farther from their galactic centers (when applicable) and being able to better characterize the stellar background distributions. 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.

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. False-color images (Robitaille & Bressert, 2012) of the two drizzled fields are shown in Figure 2. 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 two sizes of circular apertures. After analyzing flux counts from different apertures, we found that a four-pixel radius was both large enough to capture the concentrated stellar light and small enough avoid the inclusion of light from neighboring sources. By extending the aperture radius by two pixels to create a six-pixel radius, we could determine whether a source was more extended and thus not a star.

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 3, 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 maps are orientated such that north is in the direction of the positive $y$-axis and east is in the direction of the negative $x$-axis. The left panel of Figure 4 shows the 10 levels as filled contours for Peg III. There is an overdensity in the northwest (NW), however it is not in the direction of Psc II, which is indicated by the white arrow pointing from the center to the east. The contour map for Psc II is illustrated in the right panel of Figure 4, 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.

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 from North through East ($\theta$), and background surface density (average density of stars in the field not belonging to the galaxy) ($\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:

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:

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. We chose this bin size as it was large enough that the highest surface density areas of these sparse galaxies contained more than one star per bin and small enough that it would not obscure any structural features. The corner plots (Foreman-Mackey, 2016) and posterior distributions for the exponential fits are shown in Figures 5 and 6 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 7 and 8. For Peg III, the two most prominent residuals are in the NW and southeast (SE), with another to the southwest (SW). These show that a pure elliptical model cannot describe the full morphological complexity of Peg III. Psc II has prominent residuals to the north and NW, as well as one near the center of the image. We expect such small residual variation to be present even for models that reproduce the galaxy morphology.

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 ($S_{0}$), as $S_{0}$ is proportional to the richness divided by the squared scale radius, which in this case is the elliptical half-light radius. As the radius increases, for the $S_{0}$ to remain the same, the richness must also increase. 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 9. 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.

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 10. 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. Both the Kim et al. (2015, 2016) structural fits show the semimajor axis extending in the SE to NW direction, while our fits extend more along the east to west axis. The NW overdensity in our Figure 4 and smoothed data and residuals in Figure 7 could be have affected the Kim et al. fits. The same best-fit 2D-model comparisons are shown for Psc II, with the exponential fits shown in the middle panel and the Plummer fits in the right 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.

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

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

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$, as this encompassed the majority of stars from the CMD that would be expected to belong to the galaxies. 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 3 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.

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.

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.

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. While our mass estimate for Peg III is greater than that of Psc II, Kim et al. (2016) measured Peg III to be $0\farcm 85$ and thus smaller than Psc II, which could explain why they suggested Peg III could have been tidally affected by Psc II.

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., 2022), 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. (2022) 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. (2022) 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. (2022) by a factor of two or more. However, unlike Li et al. (2021) and Battaglia et al. (2022), 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. (2022) 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.¹¹1Following the initial submission of this paper, Pace et al. (2022) measured new systemic proper motions for Peg III and Psc II using the structural parameters reported herein. Their measurements agree with the McConnachie & Venn (2020) values used in our orbital analysis within one to two sigma.

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 moduli²²2During 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_⊙, which was used as the fiducial model in Patel et al. (2020) and is consistent with abundance-matched halo estimates and recent Magellanic system models (e.g., Guo et al., 2011; Besla et al., 2012, 2013, 2016) based on empirical measurements and dynamical mass arguments (e.g., Kim et al., 1998; Majewski et al., 2009; Saha et al., 2010; van der Marel & Kallivayalil, 2014; Mackey et al., 2016). 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 12. 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.

Orbital histories with respect to the MW (top panel, Figure 12) 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 (Patel et al, in prep.). Indeed, in our density contour maps (Figure 4) 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 12 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 12. 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 pericenter³³3Note that some studies will adopt the minimum relative distance when no critical minimum is recovered (e.g., Battaglia et al., 2022) 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 (D’Souza & Bell, 2022, Patel, 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).

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 13). 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 13, 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 13, 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.

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 4), we find no signs of an irregular shape for Peg III (unlike Kim et al. 2015, 2016) or for Psc II, which we calculated to be the less massive of the two UFDs and therefore possibly more susceptible to tidal effects. 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 12 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 12. 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.

Based on Gaia eDR3 measurements alone, it is unlikely for the possibility of an association between Peg III and Psc II to be further constrained. Gaia will operate much longer, however, 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. If a velocity difference less than $\sim$50 km s^-1 was observed, we would be able to exclude the possibility of a chance association at a confidence level of at least 95% (see Appendix A for detailed calculations).