The Hubble Space Telescope Survey of M31 Satellite Galaxies. I. RR Lyrae–based Distances and Refined 3D Geometric Structure

In this section we summarize the main steps of our RR Lyrae detection and modeling algorithm and describe the broad properties of the RR Lyrae sample. We will present an in-depth description of the variable star sample and analysis (e.g., variable star completeness quantification, non–RR Lyrae variables) in the main survey paper. Here, we provide a description of our procedure as applied to putative RR Lyrae stars.

We developed an approach to address two main challenges provided by our data set: the high number of RR Lyrae candidates (several thousands) and the modest number of photometric epochs that many of our target galaxies have (between 10 and 15 per filter, for roughly half of our sample; see Table 1). This required the development of a procedure that could be carried out mostly in an automated fashion, while being robust enough to reliably model light curves in a sparsely sampled regime. We devised the following multistep approach:

• 1.  
  We first identified all possible variable stars candidates following the procedure described in Dolphin et al. (2001, 2004), which selects variable candidates on the basis of four criteria, namely: (1) the rms of the magnitude measurements, (2) the amount of crowding contamination reported by DOLPHOT, (3) the χ² of the PSF fits, and (4) a variability metric based on Lafler & Kinman (1965). This step in the procedure is commonly used in the literature (e.g., McQuinn et al. 2015, for further details).
• 2.  
  We then inspected the variable sources on the CMD, where we selected likely RR Lyrae candidates, i.e., all variable stars lying within roughly half a magnitude from the mean magnitude of the HB and within the approximate color range of the instability strip, i.e., 0.45 ≲ (F475W − F814W) ≲ 1.1 and 0.15 ≲ (F606W − F814W) ≲ 0.7, adjusted for extinction (Schlafly & Finkbeiner 2011).
• 3.  
  For each star, we made a first estimate of the pulsation period using the peak of the variability merit function defined by Saha & Vivas (2017). This method is optimized for the analysis of sparse multiband light curves and combines metrics from the two main families of periodicity analysis, Fourier-based periodograms and phase-dispersion minimization, to boost the pulsation signal relative to the artifacts introduced by the discrete data sampling.
• 4.  
  We used the period measured in the previous step as a first guess for a multiband fitting of the HST light curve using 57 empirical templates derived from the extensive ground-based observations presented in Monson et al. (2017). Our light-curve models are parameterized by the pulsation period P, as well as a set of filter-dependent pulsation phases (Φ_F606W/F475W and Φ_F814W), amplitudes (A_F606W/F475W and A_F814W), and intensity-averaged magnitudes (M_F606W/F475W and M_F814W). We performed the light-curve fit in the native HST bands. We constructed our models using ground-based templates in the most closely matching filter to our HST observations, i.e., Johnson B for F475W, Johnson V for F606W, and Johnson I for F814W. The effective wavelengths of the ground-based filters are sufficiently similar to the HST counterparts that they result in similar light-curve shapes, while most of the difference between ground-based and HST light curves is effectively bypassed by our free parameters (amplitudes and mean magnitudes). Using a Gaussian likelihood function and noninformative priors (see Table 2), we sampled the posterior probability distribution (PPD) for each star using the affine invariant ensemble Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013). We defined the convergence length of the MCMC chain as 50 times the autocorrelation length. This sampling provides constraints on each star's variability parameters (defined as the 50th percentile of the PPD), the pulsation mode (RRab or RRc, based on the period and the best-fitting template), and a full uncertainty characterization of our measurements (15.9th and 84.1th percentiles of the PPD).
• 5.  
  We inspected the output sample using several diagnostic diagrams such as the period–amplitude diagram and the amplitude–amplitude diagram. We flagged ∼10% of the sample as anomalous sources in at least one of these diagrams. We followed up with manual inspection and refinement of the fit, when possible. We discarded sources that show obvious modeling inconsistencies or ambiguous period solutions (e.g., significant multimodality in the variability function of step (iii)). For the purpose of distance determinations, we prioritized purity over completeness, but for forthcoming population studies, we will revisit completeness (e.g., using artificial variable stars). As our light-curve sampling is not always sufficient to recover secondary pulsation modes, we did not attempt to identify RRd variables, which were either assigned to the RRab/RRc sample (based on the dominant pulsation mode) or discarded due to the poor light-curve fitting. We discuss this in more depth in the main survey paper.

Step (i) of this procedure produced a catalog of 6190 candidate variable sources. Application of the subsequent steps ((ii)–(v)) refined this to 4775 bona fide RR Lyrae, which we use for the distance determinations. We detect RR Lyrae in each system analyzed, ranging from a minimum of four (And xix) to a maximum of 712 (Peg DIG).

Table 2. Parameters and Priors Used in the Light-curve Modeling

Parameter	Prior	Description
P (days)	${ \mathcal U }(0,1.2)$	Period (days)
Φλ	${ \mathcal U }(0,1)$	Pulsation phase
Aλ	${ \mathcal U }(0,1.5)$	Pulsation amplitude (mag)
Mλ	${ \mathcal U }(23,26)$	Intensity-averaged magnitude

Note. Note that Φ, A, and M are parameters for each filter.

Download table as:  ASCIITypeset image

Figure 1 shows CMDs for three galaxies (And vii, And xv, and And xi), and the RR Lyrae stars we identified in them. These examples illustrate the data quality and RR Lyrae population size across the full galaxy luminosity range of our sample. The CMDs highlight the excellent photometric quality of our data. At the magnitude of the HB (m_F606W ∼ 25), the typical photometric error for nonvariable sources is 0.01 mag, for single-epoch variables 0.03 mag, and the photometric completeness within our observed fields is virtually 100%.

Zoom In Zoom Out Reset image size

Figure 2 shows a range of example light curves and best-fit models, selected to illustrate the range of data and fit quality. Even in the limit of our sparsest sampling (e.g., And xii), our light-curve fits are generally good. We constrain periods to a typical precision of 0.01 day, pulsation amplitudes to within 0.11 mag, and the mean intensity-averaged magnitudes to 0.04 mag.

Zoom In Zoom Out Reset image size

Figure 3 shows that we generally have robust fits even in the sparsely sampled regime. Here, we show the PPDs for our fit to And xii–V004614 (bottom light curve in Figure 2), a typical example of our lowest sampling regime. The posteriors are generally well constrained and unimodal, except for the phase Φ. The purpose of this phase term is to ensure that, for any given period, the maximum-light epochs of model and data are matching. Therefore, Φ is effectively a nuisance parameter over which we marginalize. In our lowest light-curve sampling regime, multimodality in the other variability parameters occasionally occurs due to trade-offs between light-curve cadence and variable star period. Most of these stars were flagged as anomalous in our step (v) and re-fit or removed. As a guard against poorly modeled stars that made it past this step, we used a Gaussian mixture model (GMM; see Section 4.2) to mitigate the impact of contaminants on our distance determinations.

Zoom In Zoom Out Reset image size

Figure 4 shows the period–amplitude diagrams (Bailey diagrams; Bailey et al. 1919) for our sample. This diagram is a widely used diagnostic of RR Lyrae properties (e.g., Sandage 1981; Soszyński et al. 2009; Clementini et al. 2019), and our ability to recover well-defined fundamental mode (RRab, shown in purple) and first-overtone (RRc, in cyan) sequences illustrates the reliability of our light-curve modeling scheme.

Zoom In Zoom Out Reset image size

This gallery of Bailey diagrams is the largest homogeneously derived gallery of its type and the first for a nearly complete satellite system. It shows the sheer diversity of the RR Lyrae populations, both in size, RRab/RRc ratio, and morphology of the period–amplitude sequences. In terms of RR Lyrae population size, there is a general trend between decreasing galaxy luminosity and decreasing number of RR Lyrae stars. However, at fixed galaxy luminosity, there is substantial scatter in the number of RR Lyrae stars, and several factors may be contributing.

As RR Lyrae stars are primarily the manifestation of old stars that live in a narrow temperature range on the HB (e.g., Walker 1989; Savino et al. 2020), their number abundance and pulsation properties are strongly influenced by the SFH of their host galaxy, which in turn affects the morphology of the HB (e.g., Salaris et al. 2013; Savino et al. 2018, 2019). A clear demonstration of this effect can be appreciated, for instance, by comparing NGC 147 and NGC 185. While NGC 147 is about a magnitude brighter than NGC 185, we detect roughly five times fewer RR Lyrae (124 versus 579) in the former than in the latter. This is explained by NGC 147 having a substantially younger stellar population than NGC 185 (Geha et al. 2015), and therefore a much lower RR Lyrae specific frequency (number of RR Lyrae produced per unit luminosity). Similar SFH effects also have consequences on the RRab/RRc ratio, since efficient production of RRc pulsators requires a sufficient extension of the HB to high temperatures. On the other hand, the presence of a predominantly red HB (as in the case of some of our targets; Martin et al. 2017) will result in a deficiency of first-overtone pulsators.

Another feature we note is the degree of variation in the slope of the RRab period–amplitude sequence, which can range from substantially slanted (e.g., NGC 147) to almost vertical (e.g., And xiii). This may again be reflective of different SFHs, as galaxies with predominantly red HBs are expected to be deficient in short-period, high-amplitude pulsators and could manifest a narrower RRab period distribution than galaxies with more extended HBs. Peculiarities in the chemical enrichment history can also play a role. RRab stars of different metallicities are known to follow different period–amplitude loci, with metal-poor RR Lyrae stars being shifted to longer periods at a given amplitude (e.g., Clementini et al. 2022). Therefore, potential correlations between the metal abundance and pulsation amplitude, produced by a peculiar SFH and chemical enrichment history, could naturally alter the morphology of the period–amplitude sequence.

Of course, the populations shown in Figure 4 are also affected by observational effects. In the case of the less-populated Bailey diagrams, the inevitable stochasticity in the RR Lyrae population may introduce apparent differences, as the position of only a handful of stars can easily affect the Bailey diagram morphology. Another important effect is the limited coverage of the ACS field of view, which reduces the number of detected RR Lyrae. This is particularly relevant in galaxies with large apparent sizes. The prime example of this effect (other than the obvious cases of M31 and M33) is the strongly elongated And xix for which, due to the half-light radius of 142 (Martin et al. 2016), we detected only four RR Lyrae in spite of the relatively bright luminosity of this galaxy (M_V = −10.1 ± 0.3). The effect of a limited field of view, in the presence of spatial gradients in the stellar population properties, can also alter the measured RRab/RRc ratio. Such field-of-view effects, however, become less pronounced for fainter galaxies. For the smaller dwarfs, our ACS field generally covers two to three half-light radii and differences in the RR Lyrae demographic of galaxies with comparable luminosity (such as And xvi versus And x, or And xii versus And XXII) become representative of intrinsic differences in their stellar populations.

Because of the interplay among these numerous astrophysical and observational effects, the detailed interpretation of our Bailey diagrams requires a quantitative modeling framework that we will explore once we have measured SFHs for our full sample of galaxies (A. Savino et al. 2022, in preparation).

The analysis procedure described in Section 3.1 was designed to maximize fit robustness on sparsely sampled light curves. Even so, retrieving accurate pulsation properties of variable stars with few photometric observations remains a challenging task. The most common problem that can arise in a low light-curve sampling regime is the contamination from period aliases. These are spurious variability signals that arise from the discrete nature of the photometric observations and, if strong enough, can lead to an incorrect period being assigned to the star. A high enough incidence of period aliases would therefore be detrimental for our distance determination accuracy.

Quantification of period alias incidence is usually done by applying the variable star analysis to light curves of known period. In our case, we leveraged the large range of observation cadence in our sample galaxies and used stars with highly sampled light curves as templates, degrading them to simulate the range of observing conditions in our most sparsely sampled galaxies. Specifically, we selected three variable stars belonging to M31, chosen to represent different light-curve properties: a high-amplitude RRab star (V005030, P = 0.559 day, A_F606W = 1.008), a low-amplitude RRab star (V005674, P = 0.734 day, A_F606W = 0.436), and an RRc star (V005402, P = 0.283 day, A_F606W = 0.499). We took the well-measured light curves of these stars (58 epochs in F606W and 60 epochs in F814W) and selected random subsets of data with variable length, ranging from 9–20 photometric epochs per filter. The photometric time-series were selected to span a small range in observation epoch and without any significant time gap, to simulate typical HST data. For each number of photometric epochs, we generated 200 random light-curve realizations.

We then analyzed the simulated light curves as we did with our observed sample, following the procedure of Section 3.1, and quantified the amount of solutions affected by a period alias. Within this experiment, we defined a solution to be affected by an alias if the recovered logarithm of the period differs from the real one (i) in excess of measurement uncertainties and (ii) by more than 0.05 dex. We chose this value because it represents a period systematic that is comparable with the other sources of uncertainties in our distance determination. Aliases at smaller period separations may still be possible; however, they would not have any significant impact on the results of this paper.

Figure 5 shows the fraction of simulated light curves whose solution is affected by a period alias, as a function of the number of photometric epochs per filter, at different stages of our analysis pipeline. For light-curve samplings greater than approximately 15 epochs per filter, we efficiently recover the correct pulsation period, already during our initial guess (step (iii), circles in Figure 5). This highlights the robustness of the Saha & Vivas (2017) methodology, which was specifically developed to minimize period ambiguity in sparsely populated light curves. At lower sampling values, aliases in the first-guess period start increasing in frequency, reaching an incidence between 5% and 15% in the sparsest sampling regime of our data (9–11 epochs per filter). This fraction is only marginally reduced by the template fitting analysis (step (iv), squares), as the algorithm preferentially samples the parameter space in the vicinity of the initial period. Many of the spurious solutions, however, result in period–amplitude combinations that are in contrast to what is expected for RR Lyrae stars, and they are therefore flagged as anomalous (step (v), triangles). This post-processing step drastically reduces the incidence of period aliases to values below 5%, at all samplings. The efficiency of step (v) in identifying aliases is in agreement with results from our real sample. In fact, out of the ∼5000 RR Lyrae candidates processed in Section 3.1, roughly 10% were flagged in step (v). Of these, approximately 75% were flagged because of an evident period alias, in line with the decrease observed, in Figure 5, between step (iv) and step (v).

Zoom In Zoom Out Reset image size

The prevalence of aliases seems to be more important for long-period, low-amplitude RRab stars, for which it occurs at roughly double the rate of the other variable templates. This is reasonable as, at fixed light-curve sampling, observations will cover fewer pulsation cycles (or even just a part of the pulsation cycle) in long-period variables, therefore increasing the chance of misclassification. Furthermore, high-amplitude pulsators are known to exist predominantly as short-period RRab stars, and will manifest clearly unphysical aliased solutions (either high-amplitude RRc solutions or long-period high-amplitude RRab). In contrast, spurious solutions in long-period RRab stars could potentially manifest as low-amplitude RRc light curves, and therefore the alias would be less easily identified.

From the results of Figure 5, we estimated that the fraction of unidentified period aliases in our final sample is likely no larger than a few percent and mostly prevalent among long-period RRab stars. This low fraction of inaccurate solutions effectively acts as a contaminating population in our distance fit and, given the low prevalence, is efficiently taken into account by the GMM formalism described in Section 4.2.

To derive robust distances to our sample, we made use of our multiband photometry to construct Wesenheit magnitudes (Madore 1982), defined as:

$$\begin{eqnarray}&&W(X,X-Y)=X-R(X-Y)\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&R={A}_{X}/E(X-Y),\end{eqnarray} \tag{ 1b }$$

where R is the total-to-selective dust absorption ratio for the X, Y filter pair. Wesenheit magnitudes, by construction, have the advantage of being reddening-free and are widely used as distance indicators, with both RR Lyrae stars and Cepheids (e.g., Neeley et al. 2019; Riess et al. 2021). They are known to be tightly correlated with period and metallicity of the variable star, according to the following functional form:

$$\begin{eqnarray}&&W=\mu +a+b\,\mathrm{log}P+c\,[\mathrm{Fe}/{\rm{H}}],\end{eqnarray} \tag{ 2 }$$

where μ is the distance modulus to the star.

With the general form of the RR Lyrae period–Wesenheit–metallicity (hereafter PWZ) in hand, we must decide on what distance anchor to use: empirical, theoretical, or a mixture of the two (e.g., Kovacs & Jurcsik 1997; Braga et al. 2015; Marconi et al. 2015, 2021; Neeley et al. 2019). Empirical calibrations have the advantage that the RR Lyrae absolute W magnitude can be tied to independent distances (e.g., via geometric parallax, through cluster CMD fitting). A particularly relevant example is that of Nagarajan et al. (2021, hereafter N21), which used a sample of 36 local RR Lyrae stars with excellent ground-based light curves and metallicities to empirically calibrate a PWZ relation anchored to Gaia eDR3 parallaxes. This is among the first to tie an RR Lyrae–based Population II distance indicator to the Gaia eDR3 reference frame.

However, empirical anchors also have limitations, particularly for the work at hand. To date, PWZ optical calibrations essentially only exist in ground-based filters (e.g., Johnson B, V, I) or for Gaia magnitudes, whereas our data is taken in the native HST photometric system. Transforming magnitudes between different photometric systems is notoriously difficult, particularly when the two systems are sufficiently different from each other. Though some HST and Johnson filters are considered similar, a close inspection of available HST–Johnson filter transformations present in the literature shows that systematic differences of up to 0.1 mag can arise between different empirical and/or synthetic transformations (Sirianni et al. 2005; Bernard et al. 2009; Saha et al. 2011; McQuinn et al. 2015). Such differences are amplified by the use of W magnitudes, which are defined as a linear combination of different filters. Adopting R = 0.960 (F475W/F814W), a systematic of 0.1 mag in one of the two filters could translate to a difference in W of up to ∼0.2 mag, while a value of R = 1.785 (F606W/F814W) could result in systematics up to the ∼0.3 mag level. As a reference, our median uncertainty on W is 0.11 mag. As Riess et al. (2021) note, ground-to-HST filter transformations are one of the dominant sources of uncertainty on galaxy distances in the limit of excellent data, such as ours.

In comparison, theoretical PWZ calibrations provide a way to mitigate this problem. Bolometric luminosities of model RR Lyrae can be readily translated into any photometric system of choice, eliminating the systematics inherent to empirical transformations. One of the most widely used theoretical calibrations is provided by Marconi et al. (2015, hereafter M15), who use a grid of nonlinear pulsation models to explore the dependency of W magnitudes on a large range of stellar mass, period, and metallicity. Such calibrations, opportunely mapped onto the ACS/WFC photometric system (M. Marconi, private communication), constitute the basis of our distance determinations. The PWZ coefficients we used in this work are listed in Table 3. Our distances were determined using only data and PWZ coefficients for the RRab stars. In Section 5.1 we will examine the effect of this choice and discuss modeling of RRc stars.

A drawback of theoretical calibrations is that they are less easily informed by independent distance anchors, resulting in potentially substantial differences from what is measured through, e.g., parallaxes. The work from N21, for instance, reports that the M15 PWZ systematically overpredicts distances compared to Gaia eDR3. This can be easily appreciated by comparing the PWZ coefficients that both authors provide for the W(I,V-I)-based PWZ. As the transformations to convert bolometric luminosities into Johnson magnitudes are among the most well established, it is reasonable that the dominant source in the predicted W(I,V-I) difference, for a given period and metallicity, lies in the absolute luminosity of the theoretical pulsation models.

Under this assumption, we derived an empirical correction term, defined as the difference between the W(I,V-I)-based PWZ of M15 and that of N21, and used it to recalibrate the HST PWZs to be consistent with the Gaia eDR3 distance scale. The functional form for this correction term is:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}\mu & = & 0.23\mathrm{log}P-0.02[\mathrm{Fe}/{\rm{H}}]\\ & & -0.087(V-I).\end{array}\end{eqnarray} \tag{ 3 }$$

Figure 6 illustrates the overall size of the empirical correction term as a function of period and [Fe/H]. For periods typical of RRab pulsators (−$0.35\lt \mathrm{log}P\lt 0$), μ decreases by 0.02 to 0.1 mag, in line with the findings of N21. For the illustrative purposes of Figure 6, we have assumed a linear correlation between the period and the (V − I) color, so that a period range of −$0.6\lt \mathrm{log}P\lt 0$ would map onto the color range 0.3 < (V − I) < 0.6. The color term in Equation (3) arises from the fact that M15 and N21 use different values of R in their definition of W(I, V − I). To calculate this color term in Equation (3), we transformed the intensity-averaged F606W/F814W magnitudes to V/I using the prescription of Saha et al. (2011) and the F475W/F814W magnitudes to V/I using that of McQuinn et al. (2015). The presence of this color term implies that a transformation from HST magnitudes to Johnson was still involved in our procedure, which was otherwise performed entirely in the native HST photometric system. However, the small coefficient of this color term means that any systematic arising from the filter transformations affects the final distance only by a factor of the order of 0.01 mag. Similarly, the quoted random uncertainties in the filter transformations only contribute a term of the order of 0.003 mag to the final distance error budget.

Zoom In Zoom Out Reset image size

We measured the distance to each galaxy using the periods and mean magnitudes of individual RR Lyrae as determined in Section 3, the PWZ from Section 4.1, and a GMM formalism, which provides a way to account for the effect of outlier RR Lyrae stars on the galaxy distance without resorting to hard cuts (e.g., sigma-clipping). We followed the formalism for a GMM as laid out in Hogg et al. (2010) and Foreman-Mackey et al. (2014).

For each RR Lyrae in a given galaxy, we used the measured mean magnitudes and period to calculate a distance modulus μ_k by applying Equations (1), (2), and (3). We assumed the measured μ_k to be sampled from a normal distribution ${ \mathcal N }(\mu ,{\sigma }_{k}^{\mu })$, where μ is the true distance modulus of our galaxy and ${\sigma }_{k}^{\mu }$ is determined, for each star, through propagation of the measurement uncertainties (P and W) and of the uncertainties in the PWZ coefficients. The metallicity of each RR Lyrae, needed to use Equations (2) and (3), is not known (i.e., they are too faint for metallicity determinations), so we left it as a free parameter, with the assumption that the value is the same for all RR Lyrae stars in the galaxy.

Despite our extensive efforts to remove contaminants from our RR Lyrae sample, our sample may not be 100% pure. For example, stars with few epochs could have incorrect periods due to aliasing or could be non–RR Lyrae stars that are misidentified, both of which would introduce contamination into the sample, possibly affecting the distance determinations.

To account for contamination, we adopted a GMM to model a second population that is drawn from ${ \mathcal N }({\mu }_{\mathrm{false}},{\sigma }_{\mathrm{false}})$, where μ_false and σ_false are both free nuisance parameters of our model. We did not enforce a binary decision on whether a given measurement belonged to either population, but rather assigned to each star a continuous probability, Q_k , of being a contaminant. We modeled the probability function as a sigmoid-like function:

$$\begin{eqnarray}&&{Q}_{k}=\displaystyle \frac{1}{1+\exp (-s({R}_{k}-2))},\end{eqnarray} \tag{ 4 }$$

where ${R}_{k}=| {\mu }_{k}-\overline{\mu }| /{\sigma }_{k}^{\mu }$ is the absolute deviation, expressed in units of ${\sigma }_{k}^{\mu }$, from the weighted average of the μ_k distribution. We chose a sigmoid so that measurements in the vicinity of the observed period–Wesenheit sequence contribute fully to the fit, while scattered outliers have significantly reduced constraining power. The probability Q_k is defined to be equal to 0.5 at R_k = 2 and to approach 1 as R_k tends to infinity. The parameter s sets the slope of the sigmoid cutoff, and it was left free to accommodate galaxy-by-galaxy differences in data quality and contamination level. Within this framework, the likelihood of a given parameter combination can be written as:

$$\begin{eqnarray}&&\begin{array}{l}p({\mu }_{k}| \mu ,[\mathrm{Fe}/{\rm{H}}],{\mu }_{\mathrm{false}},{\sigma }_{\mathrm{false}},s)\\ \,\,=\,(1-{Q}_{k})\exp \left(-{\left({\mu }_{k}-\mu \right)}^{2}/2{\left({\sigma }_{k}^{\mu }\right)}^{2}\right)\\ \,\,\,\,+\,{Q}_{k}\,\exp \left(-{\left({\mu }_{k}-{\mu }_{\mathrm{false}}\right)}^{2}/2{\left({\sigma }_{\mathrm{false}}\right)}^{2}\right).\end{array}\end{eqnarray} \tag{ 5 }$$

We explored this parameter space with emcee, using the same convergence criterion defined in Section 3.1, the sum of the log-likelihoods of all RRab stars as our merit function, and flat, uninformative priors on our parameters. Specifically, we defined the priors as a top-hat function, whose limits are listed in Table 4. For μ, μ_false, σ_false, and s, we chose the width of the top-hat to be large enough to accommodate extensive exploration of our MCMC walkers. The prior on [Fe/H] is flat for −2 < [Fe/H] < −1 and zero outside of this interval. The reason for this choice is that our data have little to no informative power on [Fe/H], meaning that, in virtually every fit, the recovered PPD on [Fe/H] closely tracks the prior. Choosing a top-hat prior means that we effectively marginalized over a flat metallicity posterior, while ensuring that [Fe/H] remains within a physically motivated interval. This essentially captures the effect of our ignorance on the metallicity in the final distance uncertainties. We chose the limits of [Fe/H] to be representative of the expected metallicity range predicted by stellar evolution theory (e.g., Savino et al. 2020) and observed in local dwarf galaxies (e.g., Clementini et al. 2005; Bernard et al. 2008).

Table 4. Parameters and Priors Used in Our PWZ Modeling

Parameter	Prior	Description
μ	${ \mathcal U }(23,26)$	Distance modulus of the galaxy
[Fe/H]	${ \mathcal U }(-2,-1)$	Metallicity of the RR Lyrae population
μfalse	${ \mathcal U }(20,30)$	Mean distance modulus of the contaminants
σfalse	${ \mathcal U }(0,10)$	Scatter of the contaminant population
s	${ \mathcal U }(1,10)$	Steepness of the Q sigmoid

Download table as:  ASCIITypeset image

Figure 7 shows examples from our PWZ models, illustrating different regimes of data quality within our sample: a clean, populated sample (Figure 7(a), And i); a populated, noisy sample (Figure 7(b), And vii); and a sparsely populated sample (Figure 7(c), And xix). However, the large numbers of RR Lyrae stars mean that we are able to achieve robust PWZ fits even in the presence of a modest population of contaminants or uncertain pulsation properties on individual stars, consistent with expectations for simple model constraints from noisy data (e.g., Hogg et al. 2010). The distance moduli measured through our PWZ models (50th percentile of the PPD), and the respective uncertainties (15.9th and 84.1th percentiles of the PPD), are summarized in Table 5.

Zoom In Zoom Out Reset image size

Table 5. Our Measured Distance Moduli and the Resulting Distances Relative to Both the Sun and M31.

ID	Nab	Nc	μ	D⊙	DM31	QPlane	MV	rh	References
				kpc	kpc			pc
M31	28	21	24.45 ± 0.06	776.2${}_{-21}^{+22}$	0	⋯	⋯	⋯	⋯
GSS	16	5	24.58 ± 0.07	${824.1}_{-26}^{+27}$	${53.4}_{-26}^{+30}$	⋯	⋯	⋯	⋯
M32	145	49	24.44 ± 0.06	${772.7}_{-21}^{+22}$	${20.6}_{-13}^{+21}$	0.953	−16.8 ± 0.1	106 ± 12	(1)
M33	47	10	24.67 ± 0.06	${859.0}_{-23}^{+24}$	${226.7}_{-11}^{+15}$	∼0	⋯	⋯	⋯
NGC147	90	34	24.33 ± 0.06	${734.5}_{-20}^{+21}$	${107.0}_{-8.0}^{+15}$	0.856	−16.6 ± 0.07	${1431}_{-43}^{+44}$	(2)
NGC185	387	192	24.06 ± 0.06	648.6 ± 18	${154.1}_{-21}^{+23}$	0.931	−15.6 ± 0.07	555 ± 17	(2)
NGC205	261	145	24.61 ± 0.06	835.6 ± 23	${58.0}_{-29}^{+30}$	0.919	−16.7 ± 0.1	${598}_{-29}^{+30}$	(1)
And i	181	53	24.45 ± 0.05	776.2 ± 18	${48.0}_{-3.2}^{+10}$	0.980	−11.4 ± 0.2	${880}_{-30}^{+31}$	(3)
And ii	152	51	24.12 ± 0.05	${666.8}_{-15}^{+16}$	${168.9}_{-16}^{+19}$	0.017	−11.7 ± 0.2	${1028}_{-30}^{+31}$	(3)
And iii	82	23	24.29 ± 0.05	${721.1}_{-16}^{+17}$	${84.9}_{-14}^{+19}$	0.188	−9.5 ± 0.2	420 ± 43	(3)
And V	122	39	24.40 ± 0.06	758.6 ± 21	${110.5}_{-3.5}^{+7.0}$	0.002	−9.3 ± 0.2	${353}_{-24}^{+35}$	(3)
And vi	203	55	24.60 ± 0.06	831.8 ± 23	${281.6}_{-7.1}^{+8.6}$	∼0	−11.6 ± 0.2	489 ± 22	(4)
And vii	418	145	24.58 ± 0.06	824.1 ± 23	${230.8}_{-6.5}^{+8.4}$	∼0	−13.3 ± 0.3	815 ± 28	(4)
And ix	41	11	24.23 ± 0.06	${701.5}_{-19}^{+20}$	${82.0}_{-24}^{+26}$	0.339	−8.6 ± 0.3	${408}_{-42}^{+62}$	(3)
And x	48	16	24.00 ± 0.06	${631.0}_{-17}^{+18}$	${162.2}_{-24}^{+25}$	0.057	−7.3 ± 0.3	${202}_{-37}^{+74}$	(3)
And xi	15	9	24.38 ± 0.07	${751.6}_{-24}^{+25}$	${104.2}_{-4.2}^{+11}$	0.669	−6.4 ± 0.4	131 ± 44	(3)
And xii	7	5	${24.28}_{-0.07}^{+0.08}$	${717.8}_{-23}^{+27}$	${107.7}_{-13}^{+20}$	0.762	−6.6 ± 0.5	${376}_{-147}^{+251}$	(3)
And xiii	11	3	24.57 ± 0.07	${820.4}_{-26}^{+27}$	${126.4}_{-8.0}^{+16}$	0.803	−6.8 ± 0.4	${191}_{-72}^{+96}$	(3)
And XIV	90	22	24.44 ± 0.06	${772.7}_{-21}^{+22}$	${160.8}_{-4.2}^{+3.8}$	0.835	−8.6 ± 0.3	337 ± 46	(3)
And xv	75	34	24.37 ± 0.05	748.2 ± 17	${95.8}_{-4.8}^{+12}$	0.004	−8.4 ± 0.3	283 ± 23	(3)
And xvi	3	4	23.57 ± 0.08	517.6 ± 19	${280.0}_{-27}^{+26}$	0.874	−7.5 ± 0.3	239 ± 25	(3)
And XVII	27	20	24.40 ± 0.07	${758.6}_{-24}^{+25}$	${49.9}_{-5.8}^{+17}$	0.798	−7.8 ± 0.3	315 ± 68	(3)
And xix	4	0	${24.55}_{-0.08}^{+0.09}$	${812.8}_{-29}^{+34}$	${113.3}_{-6.9}^{+18}$	0.001	−10.1 ± 0.3	${3357}_{-465}^{+816}$	(3)
And XX	12	3	24.35 ± 0.08	${741.3}_{-27}^{+28}$	${128.4}_{-5.5}^{+12}$	∼0	−6.4 ± 0.4	${86}_{-22}^{+43}$	(3)
And xxi	21	2	${24.44}_{-0.07}^{+0.06}$	${772.7}_{-25}^{+22}$	${124.4}_{-3.8}^{+5.1}$	∼0	−8.9 ± 0.3	${922}_{-95}^{+182}$	(3)
And XXII	16	4	24.39 ± 0.07	${755.1}_{-24}^{+25}$	${216.8}_{-5.6}^{+5.7}$	∼0	−6.4 ± 0.4	${198}_{-44}^{+66}$	(3)
And XXIII	15	3	24.36 ± 0.07	744.7 ± 24	${128.1}_{-4.9}^{+10}$	∼0	−9.8 ± 0.2	${1170}_{-94}^{+95}$	(3)
And xxiv	16	5	23.92 ± 0.07	${608.1}_{-19}^{+20}$	${194.5}_{-24}^{+25}$	0.006	−7.6 ± 0.3	${460}_{-90}^{+178}$	(3)
And xxv	18	9	${24.38}_{-0.06}^{+0.07}$	${751.6}_{-21}^{+25}$	${85.2}_{-4.4}^{+12}$	0.553	−${9.1}_{-0.2}^{+0.3}$	${590}_{-47}^{+90}$	(3)
And XXVI	21	9	${24.48}_{-0.07}^{+0.06}$	${787.0}_{-25}^{+22}$	${104.6}_{-3.5}^{+6.8}$	0.312	−${6.0}_{-0.5}^{+0.7}$	${229}_{-115}^{+138}$	(3)
And XXVIII	36	43	24.36 ± 0.05	744.7 ± 17	${368.8}_{-7.3}^{+7.8}$	∼0	−${8.8}_{-1.0}^{+0.4}$	240 ± 46	(5)
And XXIX	45	10	24.26 ± 0.06	${711.2}_{-19}^{+20}$	${189.1}_{-8.8}^{+12}$	∼0	−8.2 ± 0.4	${352}_{-42}^{+43}$	(6)
And xxx	37	15	23.74 ± 0.06	${559.8}_{-15}^{+16}$	${238.6}_{-24}^{+24}$	0.855	−${7.7}_{-0.2}^{+0.3}$	245 ± 33	(3)
And XXXI	42	15	24.36 ± 0.05	744.7 ± 17	${261.4}_{-5.9}^{+6.9}$	∼0	−11.6 ± 0.7	${910}_{-110}^{+89}$	(7)
And XXXII	71	14	24.52 ± 0.06	${801.7}_{-22}^{+23}$	${146.8}_{-4.2}^{+7.8}$	0.582	−12.3 ± 0.7	${1516}_{-237}^{+283}$	(7)
And XXXIII	35	6	24.24 ± 0.06	${704.7}_{-19}^{+20}$	${340.3}_{-8.7}^{+10}$	∼0	−10.1 ± 0.7	348 ± 83	(8)
Pisces	52	8	23.91 ± 0.05	605.3 ± 14	${292.1}_{-16}^{+17}$	0.900	−9.8 ± 0.1	370 ± 36	(9)
Peg DIG	530	182	24.74 ± 0.05	${887.2}_{-20}^{+21}$	${458.2}_{-9.4}^{+11}$	⋯	−12.3 ± 0.2	⋯	(10)
IC 1613	58	23	24.32 ± 0.05	731.1 ± 17	${511.1}_{-9.8}^{+10}$	⋯	−15.2 ± 0.1	⋯	(10)

Note. The values reported are the best fits to the RRab variables. We list the number of RRab and RRc we identify in each system. For the satellite galaxies, we indicate the probability Q_Plane of belonging to the planar structure identified in Section 6.2. We also report updated absolute luminosities and physical half-light radii for the dwarf galaxies, using the new distances. The references for the apparent V magnitudes and apparent sizes are also listed. For M32, NGC 147, NGC 185, and NGC 205, the reported radius is not the half-light radius, but rather the effective radius of a Sérsic profile with index 4 (M32 and NGC 205), 1.69 (NGC 147), and 1.78 (NGC 185).

References. (1) Choi et al. (2002); (2) Crnojević et al. (2014); (3) Martin et al. (2016); (4) McConnachie & Irwin (2006a); (5) Slater et al. (2011); (6) Bell et al. (2011); (7) Martin et al. (2013b); (8) Martin et al. (2013c); (9) Lee (1995); (10) McConnachie (2012).

Download table as:  ASCIITypeset image

For most of our targets, the dominant source of random uncertainties is the unknown metallicity of the RR Lyrae, which limits precision to roughly 0.05 mag. Uncertainties in the variability parameters and small number of RR Lyrae stars are only a significant contribution to the distance uncertainties of the least-populated (N_ab ≲ 15) galaxies, such as And xix. This is due to our choice of incorporating a broad [Fe/H] prior in our modeling, so that the random uncertainties resulting from our fits take into account the precision of our pulsation parameters, the PWZ uncertainties reported in Table 3, as well as the poorly constrained metallicity of the RR Lyrae stars.

At the heart of determining the 3D structure of the M31 satellite system is a self-consistent distance to M31 itself. From our re-reduction of archival M31 HST imaging, we apply our PWZ and distance modeling to 28 RRab variables to find μ = 24.45 ± 0.06, which corresponds to a physical distance of ${776.2}_{-21}^{+22}$ kpc. Previous analysis of the RR Lyrae stars in the same field resulted in μ = 24.48 ± 0.15 (Jeffery et al. 2011). Though this field is located in the inner stellar halo at a projected distance of 11 kpc from the photometric center of M31, under the approximation of reasonable spherical symmetry for metal-poor halo RR Lyrae, our distance value is virtually identical to the distance to M31's center.

Recent literature values of the M31 distance obtained from Cepheids are μ = 24.41 ± 0.03 (Li et al. 2021) and μ = 24.46 ± 0.20 (Bhardwaj et al. 2016), while eclipsing binary studies yield μ = 24.44 ± 0.12 (Ribas et al. 2005) and μ = 24.36 ± 0.08 (Vilardell et al. 2010). From a meta-analysis of 34 literature measurements, de Grijs & Bono (2014) derived a M31 distance modulus of 24.46 ± 0.1. Our measured value is in good agreement with these results and has been derived consistently with the distance to the other target in our sample. Therefore we adopted the distance of μ = 24.45 ± 0.06 as our anchor point for the M31 system and use it to derive relative distances to M31 (listed in Table 5) and the 3D structure of our sample. We detail the procedure for deriving 3D physical distances in Section 6.1.

The structure and morphology of the M31 satellites have been studied in great detail over past decades (e.g., Choi et al. 2002; McConnachie & Irwin 2006b; McConnachie 2012). As more satellites have been discovered, a large effort to uniformly derive structural parameters was undertaken by the PAndAS team (e.g., Crnojević et al. 2014; Martin et al. 2016). At the level of the full M31 satellite system, however, absolute luminosities and physical sizes existing in the literature are still somewhat heterogeneous, with some of the pre-PAndAS determinations being based on a range of distances and extinctions. This makes population-wide comparisons challenging.

As one step toward homogenizing the physical properties of M31 satellites, we use our new distance catalog to update some of the structural parameters (i.e., sizes and luminosities) for our galaxy sample. Though in general, distances to most galaxies do not change dramatically, there are a handful of cases in which we do see large changes compared to values used for computing sizes and luminosities (e.g., And ix, And xv, and And xxiv).

To update sizes and luminosities, we took apparent magnitudes and angular sizes from a set of literature studies chosen to be as homogeneous as possible (references given in Table 5) and converted them into absolute V luminosities and physical half-light radii using our RR Lyrae distances and extinction values from Schlafly & Finkbeiner (2011). We report the updated luminosities and sizes, the first homogeneously derived for virtually the full M31 system, in Table 5. For most galaxies, the luminosities and sizes change by a modest amount, comparable to or smaller than the measurement errors, with respect to literature values. The only galaxy for which we find substantial differences is And xxiv, for which we derive M_V = −7.6 ± 0.3 and ${r}_{h}={460}_{-90}^{+178}$ pc, compared to literature values of M_V −8.4 ± 0.4 and ${r}_{h}={680}_{-140}^{+250}$ pc (Martin et al. 2016).

We use these values throughout the remainder of this paper. A full reanalysis of the structural fits, particularly for the faintest galaxies, is the subject of future investigation as part of our Treasury program.

Our RR Lyrae sample consists of both RRab and RRc variables. However, the distances of Table 5 are based on RRab pulsators alone. We made this choice to avoid uncertainties that can arise from the more challenging theoretical modeling of RRc pulsators and to their less robust observational characterization (due to light-curve shape, shorter periods, and smaller amplitudes). Verifying the impact of this choice and quantifying how the inclusion of RRc pulsation properties affect our distances provides insight into the level of systematics arising from the stellar pulsation models. To explore this, we have run our PWZ models on the full RR Lyrae sample and fit the PWZ of RRc stars using the appropriate coefficients from Table 3, which are calibrated on the first-overtone pulsation period, P_FO . We evaluate the correction term of Equation (3) using the fundamentalized period ($\mathrm{log}{P}_{F}=\mathrm{log}{P}_{{FO}}+0.127$; e.g., Braga et al. 2015).

Figure 8 shows the difference between the distances (heliocentric and relative to M31) determined from RRab only (Table 5) and from the full RR Lyrae sample (i.e., RRab and RRc). The inclusion of RRc stars in our fit has a modest effect. Over the whole sample, the distance modulus variation is symmetric around zero, showing no significant systematic bias. The median absolute deviation in μ (top panel) is 0.013 mag (0.6% change in distance from the MW), which is generally consistent with statistical uncertainties. The largest effect is observed on AND xiii (0.082 mag; 3.8%). This galaxy represents one of the lowest signal-to-noise cases in our sample, having both the lowest number of photometric epochs (nine per filter) and a limited number of RR Lyrae (14) used for distance determination. The effect on the relative distances to M31 (bottom panel) is also generally small, with no indication of strong biases. The median absolute deviation in relative distance is 0.5 kpc (0.4% variation). The largest difference is observed for AND xxiv (13 kpc; 6.7%).

Zoom In Zoom Out Reset image size

In spite of these small differences, we found that the increased sample resulting from the inclusion of the RRc variables does not generally increase our distance precision (due to the dominant contribution of the metallicity term on the error budget). Therefore, we made the conservative choice of only measuring distances from the RRab stars and avoiding the added uncertainties of first-overtone models.

Another choice in our modeling is the application of the empirical correction term of Equation (3). The purpose of this is to bring our measurements onto the Gaia eDR3 distance scale (N21) rather than leaving them on the M15 distance scale.

Figure 9 shows differences in the distance moduli and the relative distances to M31 due to the exclusion of Equation (3) from our framework. The exclusion of the eDR3-based correction term has a nonnegligible systematic effect on the distance moduli (top panel). Retaining the M15 original distance scale systematically increases the measured distance by a median of 0.073 mag (3.4% or ∼26 kpc at the distance of M31). The effect is relatively constant among our sample, and the largest variation, observed in NGC 185, is 0.088 mag (4.1%). This is in line with expectations from Figure 6, consistent with the findings of N21. We discuss this further in Section 5.7.

Zoom In Zoom Out Reset image size

Because we self-consistently use the distance to M31 as an anchor point, the relative distances are much less affected by our choice of distance scale (bottom panel). The relative distances to M31 increase by a median amount of 4.1 kpc (3.4% variation), with a maximum change of 17 kpc in the case of the Pegasus DIG (3.3% of the 511 kpc distance from M31). As a comparison, the average random uncertainty in the relative distances is ∼10 kpc or 8%.

Our distance uncertainties are primarily driven by the unknown RR Lyrae metallicity. It is therefore useful to examine the choices we made for treatment of metallicity and its effect on distances.

Our first decision is to adopt a single [Fe/H] value for the entire RR Lyrae population. This assumption is not critical, since the RR Lyrae metallicity spreads that are quantified by observational (e.g., Clementini et al. 2005) and theoretical studies (e.g., Savino et al. 2020) are a subdominant contributor to scatter in the PWZ, compared to the measurement uncertainties on just W. To further verify this, we repeated our analysis with a 0.3 dex spread in the assumed RR Lyrae metallicity and found it affected our distances by <0.01 mag. N21 (see their Appendix E) came to a similar conclusion.

The choice of the metallicity prior has a more significant impact. We could have chosen to adopt a fixed value of [Fe/H] (e.g., Martínez-Vázquez et al. 2017; Oakes et al. 2022), or to use a tighter prior when sampling the posterior. While this would have lowered our precision floor, it would have increased the contribution of metallicity to the systematic error budget.

For example, it seems intuitive to use a prior based on the [Fe/H] that was informed by RGB star spectroscopy of a given galaxy (e.g., Collins et al. 2013) or by the luminosity–metallicity (L-Z) correlation observed in local galaxies (e.g., Kirby et al. 2013). However, such a choice is less justified than it might appear at first glance. This is because the characteristic metallicity of the RR Lyrae population is not guaranteed to be representative of the average galaxy metallicity. In fact, due to the specific effective temperature range required to trigger radial pulsation in HB stars, RR Lyrae are produced by a very specific subset of a galaxy's stellar population, whereas RGB star metallicities come from a broader range of populations. As a result, the typical metallicity of the RR Lyrae population is only weakly dependent on the galaxy metallicity distribution function.

Instead, RR Lyrae metallicities are much more sensitive to a galaxy's SFH. Counterintuitively, RR Lyrae metallicities are higher for older stellar populations. This concept was first shown by Lee (1992), and further quantified by Savino et al. (2020). More concretely, Figure 10 of Savino et al. (2020) clearly shows how the metallicity of the RR Lyrae can differ by as much as 1.5 dex from the average RGB metallicity.

This effect is unlikely to be severe for dwarfs of intermediate luminosity, for which the L-Z relation predicts average metallicities close to those inferred through the models of Savino et al. (2020). However, it can have a significant impact on the brightest (faintest) galaxies. In these regimes, the RR Lyrae production efficiency is expected to peak at much lower (higher) metallicities than the average RGB metallicity. We quantify this effect by re-running our models and swapping our flat [Fe/H] prior for a Gaussian prior, with standard deviation of 0.5 dex, and centered on the [Fe/H] predicted by the L-Z relation of Kirby et al. (2013). Over the whole sample, the median absolute change in distance modulus, compared to the values of Table 5, is 0.06 mag or ∼3% in heliocentric distance. For galaxies in the luminosity range −16.5 < M_V < −7 (∼70% of the sample), the inferred distances change by less than 0.1 mag from the values of Table 5, while for galaxies outside this interval, the change is larger. The strongest effect is observed in M31 itself, for which a prior centered on [Fe/H] = −0.4 (consistent with RGB metallicities in the inner halo, e.g., Gilbert et al. 2014), results in a distance of ∼716 kpc (Δμ = 0.18 mag), much lower than the canonically accepted range of 760–780 kpc (Section 4.3). This discrepancy is understood by considering that the metal-rich helium-burning stars in the M31 halo mostly occupy the red clump and therefore do not contribute to the RR Lyrae population, which is generated by more metal-poor stars. A prior of [Fe/H] = −0.4 for the RR Lyrae is therefore not well motivated.

The argument for M31 applies, to a lesser degree, to all galaxies for which the L-Z-based metallicity results in an inefficient RR Lyrae production. Given these difficulties in adopting a physically consistent prior on [Fe/H], we have chosen, in Section 4.2, to be agnostic about the RR Lyrae metallicity, except for the boundaries of the prior, which were chosen to limit [Fe/H] to values that are expected to produce RR Lyrae efficiently.

The TRGB is another well-established Population II distance indicator (see Beaton et al. 2018, and references therein). Due to its brightness, the TRGB can be used to much larger distances than RR Lyrae, making it excellent for mapping the Local Volume and anchoring H₀ (e.g., Tully et al. 2016; McQuinn et al. 2017; Freedman et al. 2020). However, unlike RR Lyrae, a large number of stars are necessary to clearly identify the TRGB and measure precise and/or accurate distances (Madore & Freedman 1995). In this section, we measure TRGB distances to select galaxies in our sample to better quantify the consistencies and limitations of these Population II distance indicators.

We selected a TRGB validation sample consisting of all of the galaxies that, in our ACS field of view, have a sufficiently high number of stars to allow for a robust measurement. We quantified this number by using total absolute magnitudes, sizes, ellipticities, and position angles of each galaxy (from Table 5 and references therein) to build light-profile models for our target galaxies. We used exponential light profiles for all galaxies except M32, NGC 147, NGC 185, and NGC 205, for which we used Sérsic profiles (indexes provided in the caption of Table 5). We calculated the fraction of the total galactic light captured in the ACS field of view and derived an effective absolute magnitude, ${M}_{V}^{\mathrm{eff}}$, for our photometric catalogs (e.g., an ACS field that contains 10% of the galaxy light would result in an ${M}_{V}^{\mathrm{eff}}$ 2.5 magnitudes fainter than the galaxy's M_V ). We selected galaxies that have ${M}_{V}^{\mathrm{eff}}\lt -9.5$. This corresponds to a stellar mass of roughly 10⁶ M_⊙, comparable to globular clusters such as ω Cen and 47 Tuc. This is generally regarded as the lower limit for reliable TRGB measurements (e.g., Bellazzini et al. 2001, 2004; Bono et al. 2008; Soltis et al. 2021). We excluded M32, as its projected proximity to the disk of M31 could introduce significant uncertainties in the amount of foreground extinction and contamination from M31's red giant population. We also excluded galaxies with F475W/F814W data, as there is no TRGB calibration for that filter combination in the literature. Our validation sample is therefore composed of: NGC 147, NGC 185, NGC 205, And vi, and And vii. While these galaxies are the only ones that can provide a meaningful comparison with our RR Lyrae–based distances, we also attempted to measure the TRGB in galaxies with fainter ${M}_{V}^{\mathrm{eff}}$, to explore the robustness of TRGB measurements as a function of observed stellar population size.

We measured the observed TRGB magnitude using the procedure described in McQuinn et al. (2016a, 2016b), which uses the TRGB calibration of Rizzi et al. (2007). We dereddened the photometry using the extinction map of Schlafly & Finkbeiner (2011). Then we identified the observed TRGB magnitude using a maximum likelihood modeling of the RGB luminosity function that takes into account photometric uncertainties and completeness. We define the uncertainty in the TRGB distance as the quadrature sum of the uncertainty in the apparent TRGB magnitude, the uncertainty on the Rizzi et al. (2007) calibration (as prescribed in McQuinn et al. 2017), and the uncertainty in the dust correction. We took the dust extinction uncertainty as 10% of the Schlafly & Finkbeiner (2011) extinction (as motivated in the original calibration of Schlegel et al. 1998) with a floor of 0.02 in E(B − V) (to account for the systematic uncertainties reported in the re-calibration of Schlafly & Finkbeiner 2011).

Figure 10 presents a comparison between our RR Lyrae–based distances and our TRGB distances. For all of the galaxies in the validation sample (filled dots), the two distances agree within 0.1 mag and are well within the measurement uncertainties. In this luminosity regime, we found that TRGB distance moduli are systematically smaller than RR Lyrae distance moduli, with a median difference of 0.049 mag (2.3% distance difference). Again, this effect is comparable with the random uncertainties of the PWZ distances.

Zoom In Zoom Out Reset image size

For the 16 galaxies with ${M}_{V}^{\mathrm{eff}}\,\gt$ −9.5 for which the TRGB modeling converged (empty symbols), the TRGB distance is systematically larger than that inferred through the RR Lyrae, with the difference strongly increasing with fainter ${M}_{V}^{\mathrm{eff}}$ (up to 0.7 mag, or ∼35% distance difference, for ${M}_{V}^{\mathrm{eff}}\,\approx$ −6), i.e., the magnitude of TRGB is measured to be too faint with respect to what was predicted from the RR Lyrae distance. This trend is a result of the stochastic RGB sampling in low-luminosity galaxies and it is a clear demonstration of the need of a large stellar sample for robust TRGB measurements, as first demonstrated by Madore & Freedman (1995) and further discussed by, e.g., Makarov et al. (2006), Conn et al. (2011), and McQuinn et al. (2013).

The comparison of Figure 10 provides an empirical quantification of the luminosity regime below which measured TRGB magnitudes lose fidelity as a distance indicator. In principle, well-calibrated HB magnitudes could provide a robust distance scale that extends below this luminosity threshold. However, the HB suffers from similar stochasticity problems as the TRGB and reliable HB magnitudes become hard to measure for galaxies with M_V ≳ −7. RR Lyrae stars, on the other hand, can provide robust distances even in a small-sample regime, and have been identified in stellar systems as faint as M_V = −2.7 (Martínez-Vázquez et al. 2019). Therefore, outside the MW halo, they remain the sole option to attain a precise distance scale that is consistent over a large dynamic range of luminosities and SFHs. This is a particularly important takeaway in view of the many faint stellar systems that should be uncovered by wide-area searches for galaxies in and beyond the Local Group (e.g., Simon 2019; Sand et al. 2022; searches of Roman/Rubin data).

A subset of M31 satellites have RR Lyrae distances in the literature. These are generally derived for single galaxies or small sets using a range of observations and techniques. The largest set of M31 satellite RR Lyrae distances (six galaxies; And i, And ii, And iii, And xv, And xvi, and And XXVIII) were published by Martínez-Vázquez et al. (2017) as part of the HST ISLAndS program (Skillman et al. 2017). This study derived the variable star properties using the same HST/ACS observations of this paper and obtained distances using PWZ relations based on M15 models. However, their results are based on different photometric reductions, light-curve analyses, and PWZ formalism. The distances to the Pisces and IC 1613 dwarf galaxies (Bernard et al. 2010; Hidalgo et al. 2011) were studied as part of the LCID program (Gallart et al. 2015). The observations used are also the same used in this paper (although in Pisces, we complement the F814W data with additional observations from GO-13738), but the distances were derived using a period–luminosity–metallicity (PLZ) relation. The RR Lyrae of IC 1613 were also studied in the context of the Carnegie-Chicago Hubble Program (Beaton et al. 2016), and a distance was derived using a combination of PWZ and PLZ from a completely independent HST data set (Hatt et al. 2017). Finally the mean V magnitude of the RR Lyrae population has been used as a distance proxy for And xix, And xxi, and And xxv using ground-based data (Cusano et al. 2013, 2015, 2016); for And i, And ii, And iii, And vi, And xi, and And xiii using HST/WFPC2 observations (Pritzl et al. 2002, 2004, 2005; Yang & Sarajedini 2012); and for the Peg DIG using HST/ACS data (Cole et al. 2017).

Figure 11 shows differences in RR Lyrae distances from this paper and the literature. The distances of the largest sample, that of Martínez-Vázquez et al. (2017), are systematically larger by a median of 0.06 mag. This effect is consistent with the eDR3 correction (Equation (3)), which is known to produce distances that are ∼0.06 mag closer than M15 (N21).

Zoom In Zoom Out Reset image size

The comparison with other studies shows a somewhat larger scatter, which may reflect the variety of techniques used. We find a median Δμ = −0.028 mag for the ensemble. On an individual galaxy basis, most differences are at the level of 0.1 mag and are consistent within the reported uncertainties. The only notable exception to this trend is And xxv, for which the existing literature distance is 0.25 mag larger than that determined herein. The cause of this discrepancy is not clear, but existing literature measurements for this galaxy show a similar degree of disagreement (e.g., Conn et al. 2012 measured a distance to And xxv of $\mu ={24.33}_{-0.21}^{+0.07}$, as opposed to the value of μ = 24.63 ± 0.17 derived by Cusano et al. 2016). A detailed exploration of this discrepancy is not within the scope of this paper. Overall, the scatter observed in Figure 11 exemplifies the effect of methodological differences on available literature distances, and highlights the importance of large self-consistent compilations. Such a lack of homogeneity in the literature was recently highlighted by the review of Monelli & Fiorentino (2022), who estimated systematic effects that are in good quantitative agreement with the comparison shown here.

The distances presented in this paper represent the first homogeneous analysis of nearly the full M31 system as we know it today. However, there are already a few notable compilations that comprise a significant sample of M31 satellites. Using TRGB measurements from ground-based observations, McConnachie et al. (2005) derived distances for 14 M31 satellites and M31 itself. The TRGB from ground-based observations was also used by Conn et al. (2012, hereafter C12) on 26 M31 satellites, extending the TRGB measurements to a number of fainter galaxies discovered by the PAndAS survey (Martin et al. 2006, 2009; Ibata et al. 2007; McConnachie et al. 2008, 2009, 2018; Richardson et al. 2011), while Weisz et al. (2019b) derived distances to 16 galaxies using the observed magnitude of their HB.

Figure 12 presents a comparison between our RR Lyrae–based distances and the results from these three compilations. Galaxies are plotted in order of decreasing luminosity. To guide the reader, we have indicated with a vertical line M_V = −9.5, which is the luminosity limit we have used to define our TRGB validation sample (Section 5.4). Visually, this plot reveals significant scatter among the literature compilations themselves. For the same galaxy, differences on the order of 0.2 mag are not uncommon and, in some cases, the distances differ by up to 0.5 mag. In the following sections, we discuss the trends related to this figure in more detail.

Zoom In Zoom Out Reset image size

5.6.1. Literature TRGB Distances

5.6.2. Literature HB Distances

The extensive set of internal and external comparisons presented in this section reveals the high accuracy of our distance measurements. Changes to our distance anchor (Sections 5.1, 5.2) and cross-checks with the TRGB (Section 5.4) resulted in a maximum systematic variation of 3.4% in the distances. This corresponds to a difference ∼26 kpc (0.07 mag) from the MW and of ∼4 kpc from M31. Similarly, examination of our metallicity assumptions (Section 5.3) and comparisons with literature measurements (Sections 5.5, 5.6) indicates that systematic differences are generally contained within 0.1 mag.

Compared to previous compilations, our RR Lyrae distances are a significant improvement in homogeneity and precision, enabling more robust characterizations of the M31 satellite system global properties. The largest literature M31 distance catalog, prior to this work, was that of C12, which represented a cornerstone for determinations of the spatial and structural properties of the M31 system (e.g., Conn et al. 2013; Lewis et al. 2013; Watkins et al. 2013; Martin et al. 2016; Salomon et al. 2016; McConnachie et al. 2018). Compared to C12, our work increases the galaxy sample by ∼50% and by using RR Lyrae as a distance indicator, we improve on distance precision at all galaxy luminosities. The median precision of our distances is 21 kpc (2.9%). Our largest distance uncertainty is 32 kpc (3.9%; And xix). These are substantially lower than the median error of C12 (36 kpc) and their largest uncertainty for faint galaxies, which is 190 kpc (And XXVI).

Following C12, Figure 13 shows a sky-projected view of the stellar systems in our sample as it would be observed from the center of M31, comparing the previous map from C12 (Figure 13(a)) and our new data (Figure 13(b)). The sky positions are expressed in galactocentric coordinates, with b = 90° pointing at M31's north galactic pole and l = 0° in the direction of the MW.

Zoom In Zoom Out Reset image size

This figure nicely illustrates the clear improvement that our new data set brings in terms of sample size and sky localization. Many galaxies can now be located in M31's sky within a few degrees, and only some of the satellites in the innermost halo (e.g., M32 not present in the C12 compilation) retain an uncertain sky localization. This is due to their small distance to M31, which translates small errors in distance into large errors on the polar coordinates. In the following section, we leverage our improved distances to examine select structural and spatial properties of the M31 satellite population.

One application of our precise RR Lyrae distances is precise 3D cartography of the M31 satellite system. To do this, we sample the MCMC chains of each dwarf galaxy (from Section 4.2), in parallel with the chains for M31, to combine the distance posterior probabilities and obtain a probability distribution for the coordinates of the satellite in a Cartesian reference frame centered on M31. For ease of comparison with past work, we oriented the reference frame to be identical to that of McConnachie & Irwin (2006b) and C12. In this reference frame, the XY-plane lies on the M31 disk, with the positive Z-axis pointing to M31's north galactic pole and the positive X-axis having the same azimuth, in the M31 disk, as the M31-Earth vector.

Figure 14 presents the projection of the satellite positions onto the three orthogonal Cartesian planes. The most striking feature that can be appreciated from this view is the stark asymmetry in the satellite distribution around M31. Approximately 80% of the dwarf galaxies around M31 are located in one hemisphere of the M31 halo. This hemisphere also appears to be roughly aligned with the direction of the MW.

Zoom In Zoom Out Reset image size

The asymmetry in the M31 satellites was previously identified by McConnachie & Irwin (2006b), and later by Conn et al. (2013), using TRGB distances, but it is even more pronounced with our new updated distances, as a handful of dwarfs previously beyond M31 are now on its near side (e.g., And xii, And XXII, and And xxiv). The reason for this asymmetric distribution in M31 is not known.

The alignment of this asymmetry with the MW direction raises the possibility that observational effects could play a role. Recently, Doliva-Dolinsky et al. (2022) calculated detailed dwarf galaxy detection limits for the PAndAS survey, which is responsible for finding most faint M31 satellites. They found that the limiting magnitude between the near and far side of M31's virial radius should change by no more than 1 magnitude. This means that, apart from a small sky region directly behind M31's disk, galaxies brighter than M_V ≲ −7.5 should be detectable throughout M31's halo. The persistence of strong asymmetry (at the same ∼80% level) among galaxies above this brightness indicates that, if observational effects are playing a role in the observed dwarf galaxy distribution, they are more nuanced than photometric completeness alone. A more detailed model of the spatial properties of M31 satellites, including a quantitative treatment of the spatial completeness maps, is already ongoing (A. Doliva-Dolinsky 2022, in preparation).

In the past few years, there has been some observational (e.g., Libeskind et al. 2016; Brainerd & Samuels 2020) and theoretical (e.g., Pawlowski et al. 2017; Forero-Romero & Arias 2018; Wan et al. 2020; Wang et al. 2021) evidence that at least some lopsidedness can occur in the satellite system of massive hosts, especially if a nearby massive companion (such as the MW) is present. Even so, the M31 system seems to represent a rather extreme configuration. On the other hand, M31 is essentially the only system, apart from the MW, where an accurate 3D reconstruction of the satellite-galaxy system can be obtained, while for more distant galaxy groups, only projected spatial distributions can be examined. In the original discovery paper, McConnachie & Irwin (2006b) provided a detailed discussion of possible dynamical and cosmological explanations for the M31 lopsided satellite distribution; we refer the reader to this paper for more details. Among several of these explanations, one is that much of the satellite population was accreted recently and has not had time to fully phase mix. This scenario is broadly consistent with various pictures of M31's recent accretion history (e.g., Ferguson & Mackey 2016; D'Souza & Bell 2018; McConnachie et al. 2018; Mackey et al. 2019) and could reconcile the extreme properties of M31 with respect to the broader galaxy population. Clearly, further work is needed to clarify the (lack of) consistency of the M31 dwarf galaxy system with the properties of observed and simulated satellite populations (e.g., Boylan-Kolchin 2021; Pawlowski 2021).

Another large-scale feature that has sparked discussion and controversy in recent years (see, e.g., Boylan-Kolchin 2021; Pawlowski 2021) is a degree of planar alignment among ∼50% of the dwarf satellites around Andromeda, emerging from both spatial (McConnachie & Irwin 2006b; Koch & Grebel 2006; Conn et al. 2013) and kinematic (Ibata et al. 2013; Sohn et al. 2020) data. While the lack of 6D phase-space information prevents us from firmly establishing whether the "great plane of Andromeda" (GPoA) is a long-lived feature or a transitory configuration, the reports of similar features in other local galaxies (e.g., Pawlowski et al. 2012; Müller et al. 2017, 2018) have prompted a debate on the significance of these planar structures for our current understanding of structure assembly (e.g., Kroupa et al. 2005; Libeskind et al. 2005; Zentner et al. 2005; Sawala et al. 2016, 2022; Fernando et al. 2017, 2018; Pawlowski 2018; Samuel et al. 2021).

We leveraged our larger sample of precise, homogeneous distances to reassess the properties of the GPoA. We modeled the 3D positions of our galaxy sample using a GMM analogous to that described in Section 4. For this model, we excluded IC 1613 and the Peg DIG, because their large distances from M31 could place a high weight on the recovered orientation of the GPoA, but also decrease the likeliness of physical association with the rest of the satellite system. We also excluded the GSS, since a meaningful modeling of this system in the context of the GPoA would require knowledge of the orbital plane within the GSS itself, which our single line of sight does not provide. For the remaining galaxy sample, we followed the approach of Conn et al. (2013) and used the Cartesian coordinates of Figure 14 to derive the distance of each satellite from a given plane, passing through M31, as:

$$\begin{eqnarray}&&{D}_{\mathrm{Plane}}={ax}+{by}+{cz},\end{eqnarray} \tag{ 6 }$$

where a, b, and c are the orthonormal components of the plane perpendicular vector. We assumed the plane to have a Gaussian density profile along its perpendicular axis, so that D_Plane can be described as being drawn from a population ${ \mathcal N }(0,\sqrt{{\mathrm{rms}}_{\mathrm{Plane}}^{2}+{\sigma }_{{DP}}^{2}})$, with σ_DP being the measurement uncertainty on D_Plane, and rms_Plane being the intrinsic dispersion of the planar distribution. We assumed a second, spherical isotropic population to be superimposed onto the planar component. Finally, we assumed that the radial distance D₃₁ of dwarfs that belong to this component is drawn from the distribution ${ \mathcal N }(0,\sqrt{{\mathrm{rms}}_{\mathrm{Sphere}}^{2}+{\sigma }_{D31}^{2}})$, with σ_D31 and rms_Sphere being the observational error and intrinsic dispersion of D31.

Because we use a GMM, we did not make binary decisions on which galaxies belonged to which component. Rather, we assigned each galaxy a probability of belonging to the planar component, defined as:

$$\begin{eqnarray}&&{Q}_{\mathrm{Plane}}=2(1-\mathrm{CDF}({D}_{\mathrm{Plane}})),\end{eqnarray} \tag{ 7 }$$

where CDF(D_Plane) is the cumulative distribution function, evaluated at the point D_Plane, of the Gaussian ${ \mathcal N }(0,\sqrt{{\mathrm{rms}}_{\mathrm{Plane}}^{2}+{\sigma }_{{DP}}^{2}})$. We used this probability to scale the likelihoods arising from the planar and spherical components, in a similar manner as Equation (5).

We use emcee to explore the parameter space of our model (rms_Plane, rms_Sphere), and two parameters for the plane orientation), using uniform priors, with 150 independent walkers randomly distributed in plane orientation. All of them recovered the same predominantly planar component. Viewed from the MW, this is a nearly edge-on feature comprising roughly half of the M31 satellites. Both the orientation and the probable members are in agreement with previous characterizations of the GPoA (e.g., Conn et al. 2013; Ibata et al. 2013; Pawlowski et al. 2013; Santos-Santos et al. 2020).

Given its close-to-edge-on nature, it is not surprising that the GPoA is preserved with our updated distances. However, with our measurements, we update the structural properties of the plane. We find that the inclination of the plane with respect to the MW is roughly 14°. This is larger than the original ∼ 1° value reported in Ibata et al. (2013) and more in line with the 12° figure determined by Santos-Santos et al. (2020).

The left panel of Figure 15 shows a projection of the M31 system along the direction of minimum scatter for this planar structure. Galaxies are color-coded by their probability of being plane members, Q_Plane (also provided in Table 5). Our model results in a best-fit ${\mathrm{rms}}_{\mathrm{Plane}}={23.4}_{-7.1}^{+11}$ kpc. This value, however, is somewhat dependent on the details of our model implementation. Given our choice of GMM formalism, and the dependence of Q_Plane on rms_Plane, an accurate measurement of the plane rms thickness depends on the fidelity of the isotropic contamination model. Considering the satellite anisotropy discussed in Section 6.1, and possible deviations of the radial density profile from the assumed one, a poor fit of the spherically symmetric component could result in an artificial inflation of the best-fit rms_Plane value, which should be considered an upper limit on the real dispersion of this planar structure. Indeed, if we select only the 15 satellites with Q_Plane > 0.5, we measure a dispersion of only ∼7 kpc around the best-fitting plane. This value, on the other hand, should be interpreted as a lower limit, because a selection of the highest-probability members will, by construction, return those galaxies that have the smallest distance to the plane, therefore biasing the measured thickness to lower values. Given these considerations, we find that our limits on the plane thickness are compatible with the values of 12–14 kpc reported in the literature (Conn et al. 2013; Ibata et al. 2013; Pawlowski et al.2013).

Zoom In Zoom Out Reset image size

One aspect that should be noted is that, while in this section we have referred to the GPoA as a "planar" structure, its members do not distribute homogeneously throughout the best-fit plane. This can be appreciated in the right panel of Figure 15, which shows the view of the M31 system from a direction normal to the GPoA. Here it can be seen that the high-probability GPoA members participate in the general satellite asymmetry described in Section 6.1. In fact, we find that a similar fraction of galaxies brighter than M_V = −7.5 is found on the near side of M31 in the subsamples with Q_Plane > 0.5 (82%) and Q_Plane < 0.5 (77%). In addition to the nomenclature implications on whether this structure should be best referred to as an "arc," rather than a "plane," this degree of asymmetry in the GPoA has relevance for some of its formation models, such as the tidal dwarf scenario (e.g., Metz & Kroupa 2007; Metz et al. 2007; Hammer et al. 2013; Banik et al. 2022), which requires emergence of this feature at early times and, therefore, should result in an advanced degree of phase-mixing in its constituent members.

While the long-term stability of the GPoA and its significance in the context of Λ cold dark matter have been the subject of discussion for more than a decade, major insights into the nature of this structure will ultimately be enabled by a full orbital analysis of its candidate members. Our program is laying the foundation for this orbital characterization, by providing precise distances now and by establishing an astrometric first epoch for many of the M31 galaxies that have no available proper motion. While proper motions have already been measured for a handful of M31 satellites (e.g., Sohn et al. 2020), members of our team are already in the process of obtaining second epoch imaging for additional galaxies (HST-GO-16273, PI: T. Sohn; JWST GTO-1305, PI: van der Marel). Only with full phase-space information can the plane be fully understood in a cosmological context.

In addition to a virtually complete sample of known M31 satellites, our data set includes RR Lyrae that are along a line of sight to the GSS, which is situated at a projected distance of ∼20 kpc from M31. For this specific field of the GSS, we find μ = 24.58 ± 0.07, corresponding to a distance of ${D}_{M31}={53.4}_{-26}^{+30}$ kpc behind M31 as seen from the MW. Analyzing the same field, Jeffery et al. (2011) reported a distance modulus of μ = 24.52 ± 0.19.

The GSS HST field also overlaps a Canada France Hawaii Telescope (CFHT) field of McConnachie et al. (2003, Field 7), for which they found a distance modulus of μ = 24.59 ± 0.05. The McConnachie et al. (2003) field is, in turn, contained in the GSS2 field of Conn et al. (2016), for which they reported a TRGB distance modulus of $\mu ={24.40}_{-0.02}^{+0.03}$.

It is unclear whether the difference with the distance estimated by Conn et al. (2016) represents a concern. The Conn et al. (2016) field of view covers a much larger area and averages over a much larger portion of the GSS compared to the GSS HST field. Conn et al. (2016) noted that this region of the sky is particularly difficult to model, with projected overlap between distinct features and relatively high contamination from the smooth inner halo of M31. Ultimately, this highlights the problem of interpreting the geometric properties of tidal debris even in relatively close galaxies, as low surface brightness and complex morphology can hamper the deprojection of 2D structures to a sufficient level of accuracy. In the case of M31, a more extensive survey of the RR Lyrae populations across the GSS and its neighboring structures could certainly provide precious information for a reliable 3D reconstruction, which in turn will yield deeper insight into the accretion events that produced these tidal features (e.g., Hammer et al. 2018).

NGC 147 and NGC 185 are two of the brightest satellite galaxies around M31 and in the Local Group. Due to their close proximity on the sky (∼12 kpc, projected), it has long been suspected that NGC 147 and NGC 185 formed a bound galaxy pair (e.g., van den Bergh 1998; Arias et al. 2016). Using information from radial velocities, and adopting a separation of ∼60 kpc, Geha et al. (2010) suggested that NGC 147 and NGC 185 are indeed gravitationally bound. More recently, Sohn et al. (2020) measured HST-based proper motions for these galaxies and, through orbital history modeling, were able to rule out NGC 147/185 as a bound galaxy pair. This conclusion was drawn assuming a separation of ∼89 kpc between the two galaxies.

Using our distances, we find a separation between NGC 147 and NGC 185 of ${83.8}_{-24.6}^{+25.4}$ kpc. This is only marginally smaller than previous determinations, and the result of Sohn et al. (2020) is unaffected.

As the third largest galaxy in the Local Group, M33 is expected to host its own system of satellite galaxies (e.g., Dooley et al. 2017; Bose et al. 2018; Patel et al. 2018). This prediction compares with a paucity of observed dwarf galaxies associated with M33, with the limiting magnitude of existing observations and modest coverage of the M33 virial region likely playing a significant role in the difference (Patel et al. 2018). The most likely known dwarf galaxy associated with M33 is And XXII as it has a close projected separation (∼43 kpc). However, the physical association between M33 and And XXII ultimately depends on their deprojected distance. Using CFHT photometry, Chapman et al. (2013) estimated a separation of ${59}_{-14}^{+21}$ kpc between And XXII and M33, and suggested that the two systems are currently gravitationally bound.

Using our distances, we find a separation of ${113.4}_{-29.7}^{+31.2}$ kpc between these two galaxies. While this is still within the estimated virial radius of M33 (161 kpc; Patel et al. 2017), this new measurement makes the status of And XXII as a bound satellite of M33 ultimately much more sensitive to the orbital history of the latter. Orbital solutions that exclude recent close encounters with M31 (e.g., Patel et al. 2017; van der Marel et al. 2019) would be compatible with And XXII being a current satellite of M33. However, models such as those put forward by McConnachie et al. (2009) and Putman et al. (2009) advocate for a recent <100 kpc pericenter passage of M33 around M31. They imply that And XXII, if it ever was associated with M33, is now likely unbound from its former host. Conversely, the assumption that And XXII is currently associated with M33 would disfavor the hypothesis of a close M31–M33 interaction. Ultimately, the status of And XXII, and its relevance in the context of M33's orbital history, will become much clearer when measurements of the proper motion for this dwarf galaxy become available.

With regards to the orbital history of M33, it should be noted that, while our updated distance to M31 is consistent with the value of 770 ± 40 kpc used in recent models of the M31–M33 interaction (e.g., Patel et al. 2017; van der Marel et al. 2019), the distance to M33 has changed more significantly from the value of 794 ± 23 kpc used in those models. The resulting M31–M33 distance is therefore increased by ∼25 kpc. Such change is unlikely to have profound consequences on the orbit of M33, but work has already been undertaken to include these new distance determinations into an updated dynamical modeling of the M31–M33 pair (E. Patel et al. 2022, in preparation).

Recently, a new tentative satellite of M33 was identified with the discovery of Pisces VII (Martínez-Delgado et al. 2022). With an estimated TRGB distance of ${1.0}_{-0.2}^{+0.3}$ Mpc, it is currently not clear whether Pisces VII is an isolated galaxy (like And xviii; Makarova et al. 2017) or is indeed associated with the M31–M33 system. Unfortunately, due to its very recent discovery, this galaxy is not present in our sample. The low luminosity of Pisces VII (M_V = −6.8 ± 0.2), especially in light of the results of Section 5.4, prevents a more robust localization on the basis of the TRGB alone. The only way to measure a secure distance to this galaxy is through cadenced observations like those presented in this paper.