<figure><img src="image_1.png" /></figure>

HAL
open science

# The PAndAS View of the Andromeda Satellite System.
## IV Global properties

Amandine Doliva-Dolinsky, Nicolas F Martin, Zhen Yuan, Alessandro Savino, Daniel R Weisz, Annette M.N Ferguson, Rodrigo A Ibata, Stacy Y Kim, Geraint F Lewis, Alan W Mcconnachie, et al.

<figure><img src="image_2.png" /></figure>

► To cite this version:

Amandine Doliva-Dolinsky, Nicolas F Martin, Zhen Yuan, Alessandro Savino, Daniel R Weisz, et al. The PAndAS View of the Andromeda Satellite System. IV Global properties. Astrophys.J., 2023, 952 (1), pp.72. 10.3847/1538-4357/acdcf6. hal-04042764

HAL Id: hal-04042764

https://hal.science/hal-04042764

Submitted on 9 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# The PAndAS View of the Andromeda Satellite System. IV. Global Properties

Amandine Doliva-Dolinsky<fnref n="1" /> ID, Nicolas F. Martin<fnref n="1" />,<fnref n="2" /> ID, Zhen Yuan (袁珍)<fnref n="1" /> ID, Alessandro Savino<fnref n="3" /> ID, Daniel R. Weisz<fnref n="3" /> ID, Annette M. N. Ferguson<fnref n="4" /> ID, Rodrigo A. Ibata<fnref n="1" /> ID, Stacy Y. Kim<fnref n="5" />, Geraint F. Lewis<fnref n="6" /> ID, Alan W. McConnachie<fnref n="7" /> ID, and Guillaume F. Thomas<fnref n="8" />,<fnref n="9" /> ID

Received 2023 March 1; revised 2023 June 7; accepted 2023 June 7; published 2023 July 19

## Abstract

We build a statistical framework to infer the global properties of the satellite system of the Andromeda galaxy (M31) from the properties of individual dwarf galaxies located in the Pan-Andromeda Archaelogical Survey (PAndAS) and the previously determined completeness of the survey. Using forward modeling, we infer the slope of the luminosity function of the satellite system, the slope of its spatial density distribution, and the size-luminosity relation followed by the dwarf galaxies. We find that the slope of the luminosity function is $\beta = -1.5 \pm 0.1$. Combined with the spatial density profile, it implies that, when accounting for survey incompleteness, M31 hosts $92_{-26}^{+19}$ dwarf galaxies with $M_V < -5.5$ and a sky-projected distance from M31 between 30 and 300 kpc. We conclude that many faint or distant dwarf galaxies remain to be discovered around Andromeda, especially outside the PAndAS footprint. Finally, we use our model to test if the higher number of satellites situated in the hemisphere facing the Milky Way could be explained simply by the detection limits of dwarf galaxy searches. We rule this out at >99.9% confidence and conclude that this anisotropy is an intrinsic feature of the M31 satellite system. The statistical framework we present here is a powerful tool to robustly constrain the properties of a satellite system and compare those across hosts, especially considering the upcoming start of the Euclid or Rubin large photometric surveys that are expected to uncover a large number of dwarf galaxies in the Local Volume.

*Unified Astronomy Thesaurus concepts: Andromeda Galaxy (39); Dwarf galaxies (416); Local Group (929)*

## 1. Introduction

During the last decades, faint dwarf galaxies ($L < 10^6 L_\odot$) have proven to be powerful test beds for cosmological and galaxy formation models. The majority of these constraints are obtained from the dwarf galaxy satellite system of the Milky Way (MW; e.g., Koposov et al. 2009; Kim et al. 2018; Nadler et al. 2021) because of the difficulty in detecting those faint objects beyond our immediate surroundings with current panoptic photometric surveys (e.g., Koposov et al. 2008; Drlica-Wagner et al. 2020). However, as the Milky Way satellites and past satellite accretion may not be typical (e.g., Martin et al. 2017; Weisz et al. 2019; Evans et al. 2020), it is important to explore the satellite systems of other similar hosts, the most accessible of which is the Andromeda galaxy (M31).

M31 and the cohort of dwarf galaxies that inhabit its halo are close enough (~800 kpc; Savino et al. 2022) to be resolvable into stars with modern observing capabilities. At the turn of the century, systematic efforts were undertaken to survey the surroundings of our cosmic neighbor, with the Sloan Digital Sky Survey (SDSS; Abazajian et al. 2003) and, more particularly, with the Pan-Andromeda Archaeological Survey (PAndAS;

McConnachie et al. 2018), the sample of known dwarf galaxies near M31 has increased significantly. From about 10 dwarf galaxies known at the end of the 20th century (Herschel 1789; van den Bergh 1972; Karachentsev & Karachentseva 1999), we now know of ~40 dwarf galaxies that are likely satellites of Andromeda. Among them, 4 were discovered from SDSS photometry (Zucker et al. 2004, 2007; Bell et al. 2011; Slater et al. 2011), 4 from more fortuitous efforts (Majewski et al. 2004; Irwin et al. 2008; Collins et al. 2022; Martínez-Delgado et al. 2022), 3 from searches based on the Panoramic Survey Telescope and Rapid Response System survey (Martin et al. 2013b, 2013c), and 19 from the exploration of the deeper PAndAS data (Martin et al. 2006, 2009; Ibata et al. 2007; McConnachie et al. 2008; Richardson et al. 2011). Because they are significantly more distant than their MW counterparts, these newly discovered dwarf galaxies are also somewhat brighter than the faintest MW dwarf galaxy satellites, but they nevertheless reach total luminosities as faint as $10^{4.2 \pm 0.4} L_\odot$ ($M_V = -6.0_{-0.5}^{+0.7}$ for And XXVI; Savino et al. 2022).

Beyond the mere discovery of satellites, it is essential to also quantify the completeness of those large surveys (Koposov et al. 2008; Drlica-Wagner et al. 2020) as these are key ingredients to properly fold in observational biases when comparing the known dwarf galaxies (or dwarf galaxy systems) between themselves or with simulations of the faint end of galaxy formation in a given cosmological model. This step in

<figure><img src="image_5.png" /></figure>

<footnote n="1">¹ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000, France; amandine.doliva-dolinsky@astro.unistra.fr</footnote>

<footnote n="2">² Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany</footnote>

<footnote n="3">³ Department of Astronomy, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="4">⁴ Institute for Astronomy, University of Edinburgh Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK</footnote>

<footnote n="5">⁵ Department of Physics, University of Surrey, Guildford, GU2 7XH, UK</footnote>

<footnote n="6">⁶ Sydney Institute for Astronomy, School of Physics, A28, The University of Sydney, NSW 2006, Australia</footnote>

<footnote n="7">⁷ NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada</footnote>

<footnote n="8">⁸ Instituto de Astrofísica de Canarias, Calle Vía Láctea, s/n, E-38205. La Laguna, Santa Cruz de Tenerife, Spain</footnote>

<footnote n="9">⁹ Departamento de Astrofísica, Universidad de La Laguna, E-38206, La Laguna, Tenerife, Spain</footnote>
turn requires building well-understood search algorithms that
can be run on artificial dwarf galaxies injected in the survey
data. In the case of the PAndAS survey, Martin et al. (2013a)
developed a likelihood-based algorithm that runs on the
survey's photometric catalog and determines the probability
of there being a dwarf galaxy at any location of the survey
based on the distribution of local stars on the sky and in the
color-magnitude space. Doliva-Dolinsky et al. (2022) then
used this algorithm to characterize the detection limits of the
survey, injecting half a million artificial dwarf galaxies with
varying sizes, luminosities, and positions. The resulting
detection limits show significant variations that are driven, as
expected, by the size and luminosity of the systems (i.e., their
surface brightness), but also by the location in the survey. With
PAndAS spanning more than 20° on the sky and M31 being
located fairly close to the Milky Way plane (b = -22°), the
strongly varying MW foreground contamination between the
southern edge of the survey (b ~ -35°) and its northern edge
(b ~ -11°) leads to significant variations of the surface
brightness limits (from ~30.5 mag arcsec⁻² far from the MW
plane to ~29 mag arcsec⁻² closest to the plane, respectively).
As shown in Doliva-Dolinsky et al. (2022), variations in the
foreground contamination are at least as important as changes
to the heliocentric distance of a satellite. Both lead to variations
of ~1.5 mag arcsec⁻² in the surface brightness detection limit
over the M31 halo.

With this knowledge in mind, it is possible to reliably infer
the global properties of the dwarf galaxy satellite system of
M31, taking detection limit biases into account. Among those
global properties, the shape of its luminosity function is an
important observational probe as it is sensitive to cosmology, to
feedback processes and to reionization. Indeed, the normal-
ization, shape and/or the existence of a break in the faint end of
the luminosity function is an imprint of the properties of dark
matter (Spergel & Steinhardt 2000; Bode et al. 2001) and of the
suppression of star formation from stellar outflows and
reionization (Bullock et al. 2000; Somerville 2002; Mash-
chenko et al. 2008; Koposov et al. 2009; Boylan-Kolchin et al.
2015; Wheeler et al. 2015; Weisz & Boylan-Kolchin 2017).
The radial distribution of dwarf galaxy satellites around their
host can also be shaped by the physics of reionization (Ocirk
& Aubert 2011; Dooley et al. 2017) and the disruption of
subhalos by the central disk (D'Onghia et al. 2010; Kelley et al.
2019; Samuel et al. 2020).

Another challenge comes from the distribution of satellites
that does not appear as isotropic as expected from ΛCDM
(Pawlowski 2018). Disk-like distributions of satellite dwarf
galaxies have been found around the MW, M31, and Centaurus
A (Lynden-Bell 1976; Kroupa et al. 2005; Metz et al. 2007;
Conn et al. 2013; Ibata et al. 2013; Müller et al. 2018). In
addition, when looking at the position of M31 dwarf galaxies,
most of them appear to lie closer to the MW than on the
opposite hemisphere (McConnachie & Irwin 2006; Conn et al.
2012; Wan et al. 2020). With updated RR Lyrae-based
distances, Savino et al. (2022) reenforces those conclusions,
further highlighting the anisotropy in the M31 satellite
distribution. The detection limits of dwarf galaxy searches
could lead to an anisotropy between the close and far
hemispheres of M31 as the farther the dwarf galaxy the harder
it is to detect. It is therefore essential to fold in these detection
limits when inferring the global properties of the M31 satellite

system to check if this anisotropy could simply be the results of
observational biases.

Faced with the issue of comparing biased, incomplete
observations with models, it may seem more convenient and
straightforward to simply correct observed properties, for
instance a binned luminosity function, with correction factors
calculated from the detection limits. This technique is however
plagued by noise in the case of small samples, as is the case for
dwarf galaxy systems. Therefore, while it is computationally
more expensive, it is much more reliable to forward model the
limitations of the data (detection limits, irregular shape of the
survey) directly into the model; this is the approach that we
follow here.

We use a forward-modeling approach to infer the combined
properties of the luminosity function, the radial distribution,
and the size-luminosity relation of the dwarf galaxy system of
M31. In Section 2, we detail the sample of observed satellites
and the dwarf galaxy completeness of PAndAS. Sections 3.1
and 3.2 describe the framework and the model used to obtain
the results presented in Section 4. Finally, we summarize and
discuss the main properties of the dwarf galaxy satellite system
of M31 in Section 5.

## 2. Sample

PAndAS (McConnachie et al. 2009, 2018) was conducted
from 2008 to 2011 with the 1 square degree MegaCam wide
field image at the Canada-France-Hawaii Telescope (CFHT).
Combined with previous observations (Ibata et al. 2007;
McConnachie et al. 2008), this Large Program resulted in a
survey of over 400 square degrees surrounding M31 and M33,
reaching out to ~150 kpc and ~50 kpc in projected distance
from these galaxies, respectively. For the details of the survey
and the creation of the catalogs, we refer the reader to
McConnachie et al. (2018), but it is worth mentioning that the g
and i band photometry is obtained for all fields with a median
depth of 26.0 and 24.8 for 5σ detections, respectively (Ibata
et al. 2014).

The 24 dwarf galaxies known within this footprint are listed
in Table 1. As the search algorithm struggles to separate dwarf
galaxies from stellar structures near M31 (Martin et al. 2013a)
and as the completeness in this region is not well constrained,
we choose to mask the inner 30 kpc. The luminosity and size of
each dwarf galaxy are taken from Martin et al. (2016) and
Savino et al. (2022). Where needed, distance-related properties
(physical half-light radii, absolute magnitudes) are updated
using the distances from Savino et al. (2022). Given the
uncertain nature of And XXVII that may well be a disrupted
system (Preston et al. 2019) and has large uncertainties in its
structural properties (Richardson et al. 2011; Martin et al.
2016), we choose not to add it to our sample.

The search for dwarf galaxies suffers from spatial and
photometric incompleteness. The former arises from the
complex PAndAS coverage on the sky and its correction is
quite straightforward, while the latter stems from the complex
detection process and is very sensitive to the characteristics of a
given dwarf galaxy but also to its location within the survey,
mainly because of the varying MW and M31 stellar
contamination. The detection limits were derived by Doliva-
Dolinsky et al. (2022) via the injection of nearly half a
million artificial dwarf galaxies in the PAndAS catalog to
obtain the recovery fraction for each MegaCam field on a M<sub>V</sub>
and log(r<sub>h pc</sub>) grid defined by -8.5 ≤ M<sub>V</sub> ≤ -4.5 and
<table><caption>**Table 1**
Sample of the Dwarf Galaxies Present in the PAndAS Survey</caption><thead><tr><th>Name</th><th>α(J2000)</th><th>δ(J2000)</th><th>r<sub>h</sub>(arcmin)</th><th>m<sub>V</sub></th><th>D<sub>MW</sub>(kpc)</th><th>M<sub>V</sub></th><th>r<sub>h</sub> (pc)</th><th>D<sub>M31</sub> (kpc)</th><th>Recovery Fraction</th></tr></thead><tbody><tr><td>And I</td><td>00<sup>h</sup>45<sup>m</sup>39<sup>s</sup>7</td><td>+38°02′15″</td><td>3.9<sup>+0.1</sup><sub>-0.1</sub></td><td>13.1<sup>+0.2</sup><sub>-0.2</sub></td><td>775<sup>+19</sup><sub>-17</sub></td><td>-11.4 ± 0.2</td><td>880<sup>+31</sup><sub>-30</sub></td><td>48.0<sup>+10</sup><sub>-3.2</sub></td><td>1.00</td></tr><tr><td>And II</td><td>01<sup>h</sup>16<sup>m</sup>26<sup>s</sup>8</td><td>+33°26′07″</td><td>5.3<sup>+0.1</sup><sub>-0.1</sub></td><td>12.4<sup>+0.2</sup><sub>-0.2</sub></td><td>667<sup>+16</sup><sub>-15</sub></td><td>-11.7 ± 0.2</td><td>1028<sup>+31</sup><sub>-30</sub></td><td>168.9<sup>+19</sup><sub>-16</sub></td><td>1.00</td></tr><tr><td>And III</td><td>00<sup>h</sup>35<sup>m</sup>30<sup>s</sup>9</td><td>+36°29′56″</td><td>2.0<sup>+0.2</sup><sub>-0.2</sub></td><td>14.8<sup>+0.2</sup><sub>-0.2</sub></td><td>721<sup>+17</sup><sub>-16</sub></td><td>-9.5 ± 0.2</td><td>420 ± 43</td><td>84.9<sup>+19</sup><sub>-14</sub></td><td>1.00</td></tr><tr><td>And V</td><td>01<sup>h</sup>10<sup>m</sup>17<sup>s</sup>5</td><td>+47°37′42″</td><td>1.6<sup>+0.2</sup><sub>-0.1</sub></td><td>15.1<sup>+0.2</sup><sub>-0.2</sub></td><td>759<sup>+21</sup><sub>-20</sub></td><td>-9.3 ± 0.2</td><td>353<sup>+35</sup><sub>-24</sub></td><td>110.5<sup>+7</sup><sub>-3.5</sub></td><td>1.00</td></tr><tr><td>And IX</td><td>00<sup>h</sup>52<sup>m</sup>53<sup>s</sup>4</td><td>+43°11′57″</td><td>2.0<sup>+0.2</sup><sub>-0.2</sub></td><td>15.6<sup>+0.3</sup><sub>-0.3</sub></td><td>702<sup>+19</sup><sub>-20</sub></td><td>-8.6 ± 0.3</td><td>408<sup>+62</sup><sub>-42</sub></td><td>82.0<sup>+26</sup><sub>-24</sub></td><td>1.00</td></tr><tr><td>And X</td><td>01<sup>h</sup>06<sup>m</sup>35<sup>s</sup>4</td><td>+44°48′27″</td><td>1.1<sup>+0.4</sup><sub>-0.2</sub></td><td>16.7<sup>+0.3</sup><sub>-0.3</sub></td><td>630<sup>+18</sup><sub>-18</sub></td><td>-7.3 ± 0.3</td><td>202<sup>+74</sup><sub>-37</sub></td><td>162.2+25<sub>-24</sub></td><td>1.00</td></tr><tr><td>And XI</td><td>00<sup>h</sup>46<sup>m</sup>19<sup>s</sup>7</td><td>+33°48′10″</td><td>0.6<sup>+0.2</sup><sub>-0.2</sub></td><td>18.0<sup>+0.4</sup><sub>-0.4</sub></td><td>751<sup>+23</sup><sub>-22</sub></td><td>-6.4 ± 0.4</td><td>131 ± 44</td><td>104.2<sup>+11</sup><sub>-4.2</sub></td><td>0.97</td></tr><tr><td>And XII</td><td>00<sup>h</sup>47<sup>m</sup>28<sup>s</sup>3</td><td>+34°22′38″</td><td>1.8<sup>+0.2</sup><sub>-0.7</sub></td><td>17.7<sup>+0.5</sup><sub>-0.5</sub></td><td>718<sup>+25</sup><sub>-26</sub></td><td>-6.6 ± 0.5</td><td>376<sup>+251</sup><sub>-147</sub></td><td>107.7<sup>+20</sup><sub>-13</sub></td><td>0.80</td></tr><tr><td>And XIII</td><td>00<sup>h</sup>51<sup>m</sup>51<sup>s</sup>0</td><td>+33°00′16″</td><td>0.8<sup>+0.4</sup><sub>-0.3</sub></td><td>17.8<sup>+0.4</sup><sub>-0.4</sub></td><td>821<sup>+28</sup><sub>-26</sub></td><td>-6.8 ± 0.4</td><td>191<sup>+96</sup><sub>-72</sub></td><td>126.4<sup>+16</sup><sub>-8.0</sub></td><td>0.99</td></tr><tr><td>And XIV</td><td>00<sup>h</sup>51<sup>m</sup>35<sup>s</sup>0</td><td>+29°41′23″</td><td>1.5<sup>+0.2</sup><sub>-0.2</sub></td><td>15.8<sup>+0.3</sup><sub>-0.3</sub></td><td>773<sup>+21</sup><sub>-21</sub></td><td>-8.6 ± 0.3</td><td>337 ± 46</td><td>160.8<sup>+3.8</sup><sub>-4.2</sub></td><td>1.00</td></tr><tr><td>And XV</td><td>01<sup>h</sup>14<sup>m</sup>18<sup>s</sup>3</td><td>+38°07′11″</td><td>1.3<sup>+0.1</sup><sub>-0.1</sub></td><td>16.0<sup>+0.3</sup><sub>-0.3</sub></td><td>746<sup>+17</sup><sub>-18</sub></td><td>-8.4 ± 0.3</td><td>283 ± 23</td><td>95.8<sup>+12</sup><sub>-4.8</sub></td><td>1.00</td></tr><tr><td>And XVI</td><td>00<sup>h</sup>59<sup>m</sup>30<sup>s</sup>3</td><td>+32°22′34″</td><td>1.0<sup>+0.1</sup><sub>-0.1</sub></td><td>16.1<sup>+0.3</sup><sub>-0.3</sub></td><td>517<sup>+18</sup><sub>-19</sub></td><td>-7.5 ± 0.3</td><td>239 ± 25</td><td>280.0<sup>+26</sup><sub>-27</sub></td><td>1.00</td></tr><tr><td>And XVII</td><td>00<sup>h</sup>37<sup>m</sup>06<sup>s</sup>3</td><td>+44°19′23″</td><td>1.4<sup>+0.3</sup><sub>-0.3</sub></td><td>16.6<sup>+0.3</sup><sub>-0.3</sub></td><td>757<sup>+24</sup><sub>-23</sub></td><td>-7.8 ± 0.3</td><td>315 ± 68</td><td>49.9<sup>+17</sup><sub>-5.8</sub></td><td>1.00</td></tr><tr><td>And XIX</td><td>00<sup>h</sup>19<sup>m</sup>34<sup>s</sup>5</td><td>+35°02′41″</td><td>14.2<sup>+3.4</sup><sub>-1.9</sub></td><td>14.5<sup>+0.3</sup><sub>-0.3</sub></td><td>813<sup>+31</sup><sub>-31</sub></td><td>-10.1 ± 0.3</td><td>3357<sup>+816</sup><sub>-465</sub></td><td>113.3<sup>+18</sup><sub>-6.9</sub></td><td>1.00</td></tr><tr><td>And XX</td><td>00<sup>h</sup>07<sup>m</sup>30<sup>s</sup>6</td><td>+35°07′37″</td><td>0.4<sup>+0.2</sup><sub>-0.1</sub></td><td>18.0<sup>+0.4</sup><sub>-0.4</sub></td><td>741<sup>+27</sup><sub>-27</sub></td><td>-6.4 ± 0.4</td><td>86<sup>+43</sup><sub>-22</sub></td><td>128.4<sup>+12</sup><sub>-5.5</sub></td><td>0.98</td></tr><tr><td>And XXI</td><td>23<sup>h</sup>54<sup>m</sup>47<sup>s</sup>9</td><td>+42°28′14″</td><td>4.1<sup>+0.8</sup><sub>-0.4</sub></td><td>15.5<sup>+0.3</sup><sub>-0.3</sub></td><td>770<sup>+23</sup><sub>-22</sub></td><td>-8.9 ± 0.3</td><td>922<sup>+182</sup><sub>-95</sub></td><td>124.4<sup>+5.1</sup><sub>-3.8</sub></td><td>1.00</td></tr><tr><td>And XXII</td><td>01<sup>h</sup>27<sup>m</sup>40<sup>s</sup>4</td><td>+28°05′25″</td><td>0.9<sup>+0.3</sup><sub>-0.2</sub></td><td>18.0<sup>+0.4</sup><sub>-0.4</sub></td><td>754<sup>+24</sup><sub>-23</sub></td><td>-6.4 ± 0.4</td><td>198<sup>+66</sup><sub>-44</sub></td><td>216.8<sup>+5.7</sup><sub>-5.6</sub></td><td>0.90</td></tr><tr><td>And XXIII</td><td>01<sup>h</sup>29<sup>m</sup>21<sup>s</sup>0</td><td>+38°43′26″</td><td>5.4<sup>+0.4</sup><sub>-0.4</sub></td><td>14.6<sup>+0.2</sup><sub>-0.2</sub></td><td>745<sup>+24</sup><sub>-25</sub></td><td>-9.8 ± 0.2</td><td>1170<sup>+95</sup><sub>-94</sub></td><td>128.1<sup>+10</sup><sub>-4.9</sub></td><td>1.00</td></tr><tr><td>And XXIV</td><td>01<sup>h</sup>18<sup>m</sup>32<sup>s</sup>7</td><td>+46°22′13″</td><td>2.6<sup>+0.5</sup><sub>-0.5</sub></td><td>16.3<sup>+0.3</sup><sub>-0.3</sub></td><td>609<sup>+19</sup><sub>-20</sub></td><td>-7.6 ± 0.3</td><td>460<sup>+178</sup><sub>-90</sub></td><td>194.5<sup>+35</sup><sub>-24</sub></td><td>0.92</td></tr><tr><td>And XXV</td><td>00<sup>h</sup>30<sup>m</sup>09<sup>s</sup>9</td><td>+46°51′41″</td><td>2.7<sup>+0.4</sup><sub>-0.2</sub></td><td>15.3<sup>+0.3</sup><sub>-0.2</sub></td><td>752<sup>+23</sup><sub>-23</sub></td><td>-9.1<sup>+0.3</sup><sub>-0.2</sub></td><td>590<sup>+90</sup><sub>-47</sub></td><td>85.2<sup>+12</sup><sub>-4.4</sub></td><td>1.00</td></tr><tr><td>And XXVI</td><td>00<sup>h</sup>23<sup>m</sup>46<sup>s</sup>3</td><td>+47°54′43″</td><td>1.0<sup>+0.6</sup><sub>-0.5</sub></td><td>18.5<sup>+0.7</sup><sub>-0.5</sub></td><td>786<sup>+24</sup><sub>-23</sub></td><td>-6.0<sup>+0.7</sup><sub>-0.5</sub></td><td>229<sup>+138</sup><sub>-115</sub></td><td>104.6<sup>+6.8</sup><sub>-3.5</sub></td><td>7.00 × 10<sup>-5</sup></td></tr><tr><td>And XXX</td><td>00<sup>h</sup>36<sup>m</sup>34<sup>s</sup>6</td><td>+49°38′49″</td><td>1.5<sup>+0.2</sup><sub>-0.2</sub></td><td>16.0<sup>+0.3</sup><sub>-0.3</sub></td><td>558<sup>+17</sup><sub>-16</sub></td><td>-7.7<sup>+0.3</sup><sub>-0.2</sub></td><td>245 ± 33</td><td>238.6<sup>+74</sup><sub>-24</sub></td><td>1.00</td></tr><tr><td>NGC 147</td><td>00<sup>h</sup>47<sup>m</sup>27<sup>s</sup>0</td><td>+34°22′29″</td><td>6.70 ± 0.09</td><td>7.76 ± 0.06</td><td>773<sup>+21</sup><sub>-20</sub></td><td>-16.6 ± 0.07</td><td>1431<sup>+44</sup><sub>-43</sub></td><td>107.0<sup>+15</sup><sub>-8</sub></td><td>1.00</td></tr><tr><td>NGC 185</td><td>00<sup>h</sup>38<sup>m</sup>58<sup>s</sup>0</td><td>+48°20′15″</td><td>2.94 ± 0.04</td><td>8.46 ± 0.06</td><td>650<sup>+18</sup><sub>-18</sub></td><td>-15.6 ± 0.07</td><td>555 ± 17</td><td>154.1<sup>+23</sup><sub>-21</sub></td><td>1.00</td></tr></tbody></table>

Note. The apparent magnitude and apparent size values are taken from Martin et al. (2016), except for those of NGC 147 and NGC 185 that are taken from Crnojević et al. (2014). All absolute magnitudes, physical sizes, and distances are from Savino et al. (2022). While being in the PAndAS footprint, some galaxies are not part of this sample because they are in a region where the completeness was not determined (M32, NGC 205; Doliva-Dolinsky et al. 2022), because their structural parameters are too uncertain (And XXVII; Richardson et al. 2011), or because their distances from M31 is beyond 300 kpc (And XVIII; Savino et al. 2022).

1.8 ≤ log(r<sub>h(PC)</sub>) ≤ 3.0, with a step size of 0.25 and 0.10, respectively. The recovery (or lack thereof) of a dwarf galaxy is performed with the search algorithm developed by Martin et al. (2013a) and that looks for overdensities of stars both spatially and along a red-giant-branch feature in the color-magnitude diagram. An analytical model is fitted to the resulting M<sub>V</sub>-log(r<sub>h(PC)</sub>) recovery fraction grid so that the recovery fraction of any galaxy at a given location, with a given magnitude and size, can easily be calculated. We have also built an analytical model to account for the impact of the heliocentric distance to a dwarf galaxy on recovery fractions. Although the impact of distance to a dwarf galaxy is less important than other parameters, the effect is still not negligible and needs to be taken into account (Doliva-Dolinsky et al. 2022). From these, the resulting efficiency of detection for all 24 dwarf galaxies in the sample are listed in Table 1. For the recovery fraction of dwarf galaxies with M<sub>V</sub> < -8.5 and/or log(r<sub>h(PC)</sub>) > 3, we extrapolate the analytical model for larger/brighter dwarf galaxies.

### 3. Model

Here, we discuss our methodology to infer the global properties of the M31 dwarf galaxy system from the observed properties of the individual M31 dwarf galaxies. Accounting for the PAndAS survey detection limits, we infer the underlying M31 luminosity function, size-luminosity relation, and spatial distribution via forward modeling.

#### 3.1. Dwarf Galaxy Probabilistic Model

Consider a dwarf galaxy whose observed properties, listed in Table 1, are: its coordinates on the sky, (α, δ), its apparent magnitudes in both the g and i PAndAS bands, m<sub>g</sub> and m<sub>i</sub>, its angular half-light radius, r<sub>h</sub><sup>ang</sup>, and its heliocentric distance, D<sub>MW</sub>. As detailed in McConnachie et al. (2018), the apparent magnitude are corrected for extinction following Schlegel et al. (1998) and Schlafly & Finkbeiner (2011), and m<sub>V</sub> is obtained from m<sub>g</sub> and m<sub>i</sub> using the transformation equations derived in Ibata et al. (2014). Then, using D<sub>MW</sub>, it is straightforward to transform these observed properties into the intrinsic properties of the systems: the absolute magnitude, M<sub>V</sub>, and the physical half-light radius, r<sub>h</sub>. We also use the observed properties of the dwarf galaxy D<sub>MW</sub> and (α, δ) to calculate the spherical coordinates of a dwarf galaxy in the M31-centric referential (r<sub>M31</sub>, θ, φ). The dwarf galaxy properties considered here are therefore D = {M<sub>V</sub>, log r<sub>h</sub>, r<sub>M31</sub>, θ, φ}.

We chose to define the dwarf galaxy probabilistic model that depends on a set of parameters P as the combination of three independent components: the probability of a dwarf galaxy to have a given absolute magnitude, i.e., the shape of the luminosity function of the satellite system, P<sub>LF</sub>(M<sub>V</sub>|P); the probability for a dwarf galaxy to have a given size knowing its magnitude, i.e., the size-luminosity relation of the satellite system, P<sub>log r<sub>h</sub>|M<sub>V</sub></sub>(log r<sub>h</sub>|M<sub>V</sub>, P); and the probability of a dwarf galaxy to be at a given (sky-projected or three-dimensional) location around M31, P<sub>sp</sub>(r<sub>M31</sub>, θ, φ|P). We assume that those components are independent of each other (see discussion in
<figure><img src="image_6.png" /><figcaption>Figure 1. Resulting correlation graphs and marginalized PDFs (full line) for each parameter of the 3D model. Black crosses represent the highest likelihood value for each couple of parameters. The marginalized PDFs for the 2D model are represented by the dashed lines.</figcaption></figure>

Section 4.3.1), which allows us to simply write the probabilistic model as

$$P(D|\mathcal{P}) \propto P_{\text{LF}}(M_V|\mathcal{P}) P_{\log r_h|M_V}(\log r_h|M_V, \mathcal{P}) \times P_{\text{sp}}(r_{\text{M31}}|\mathcal{P}). \quad (1)$$

Following Tollerud et al. (2008), we model the shape of the luminosity function of the satellite system as a power law with exponent $\beta$ over the magnitude range that we consider here for M31 dwarf galaxies ($M_V < -5.5$; this choice is discussed further in Section 4.3.1).

$$P_{\text{LF}}(M_V|\beta) \propto \frac{\log 10}{2.5} 10^{-(\beta+1)(M_V-4.83)/2.5}. \quad (2)$$

Following Shen et al. (2003) and Brasseur et al. (2011), we assume a linear relation between $M_V$ and the mean $\log(r_h)$,

$\langle \log r_h \rangle$, such that

$$\langle \log r_h \rangle = z_p + s(M_V + 6.0), \quad (3)$$

with $s$ being the slope and $z_p$ being the value of the relation for $M_V = -6.0$. The intrinsic dispersion, $\sigma$, around the relation is modeled as a Gaussian distribution along the $\log r_h$ direction and yields

$$P_{\log r_h|M_V}(\log r_h|M_V, z_p, s, \sigma) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{1}{2}\left(\frac{\log r_h - \langle\log r_h\rangle}{\sigma}\right)^2\right). \quad (4)$$

Finally, as the distribution of M31 satellite dwarf galaxies appears to be circularly but not spherically isotropic
(Savino et al. 2022), we consider two cases for the spatial distribution part of the model: a sky-projected (2D) and a volumetric (3D) distribution model. In both cases, we choose a simple and agnostic shape for the radial density distribution function, a power law, with parameters $\alpha_{2D}$ and $\alpha_{3D}$, respectively.<fnref n="10" /> At this stage, we introduce the assumption of an isotropic distribution of the dwarf galaxies around M31 (an assumption we will revisit later) to simplify the problem at hand. This allows us to remove the impact of the spherical coordinate angles on any model we define and, for the 3D case, we have

$$P_{\text{sp}}(r_{\text{M31}}, \theta, \phi|\alpha_{\text{3D}}) = P_{\text{sp}}(r_{\text{M31}}|\alpha_{\text{3D}}) \\ \propto \int_{0}^{2\pi} \int_{0}^{\pi} r^2 r^{\alpha_{\text{3D}}} \sin(\theta) d\theta d\phi \\ \propto 4\pi r^{2+\alpha_{\text{3D}}}. \qquad (5)$$

Similarly, for the sky-projected model,

$$P_{\text{sp}}(r_{\text{M31}}|\alpha_{\text{2D}}) \propto 2\pi r^{1+\alpha_{\text{2D}}}. \qquad (6)$$

In summary, the probabilistic model has 5 parameters $\mathcal{P} = \{\beta, z_p, s, \sigma, \alpha\}$, with $\alpha = \alpha_{2D}$ or $\alpha = \alpha_{3D}$ in the 2D and 3D cases, respectively. Folding everything together and introducing the normalization constant $A(\mathcal{P})$ to ensure that $P(\mathcal{D}|\mathcal{P})$ integrates to unity over the full space (i.e., $-17 < M_V < -5.5$, $1.8 < \log r_h < 4$ and $30 < r_{\text{M31}} < 300$ kpc), Equation (1) becomes

$$P(\mathcal{D}|\mathcal{P}) = A(\mathcal{P}) P_{\text{sp}}(r_{\text{M31}}|\alpha_{\text{3D}}) P_{\text{log}}(r_{\text{h}}|M_{\text{V}}, z_{\text{p}}, s, \sigma) \\ \times P_{\text{LF}}(M_{\text{V}}|\beta), \qquad (7)$$

with $A(\mathcal{P})$ such that

$$\iiint \mathcal{P}(\mathcal{D}|\mathcal{P}) dr_{\text{M31}} d \log r_h d M_V = 1. \qquad (8)$$

## 3.2. Final Likelihood Function

The probabilistic model presented above describes the distribution of dwarf galaxies in the data space but does not provide any constraint on the theoretical number of dwarf galaxies that inhabit the M31 satellite system, $N_{\text{true}}$, over the chosen ranges of observed properties.<fnref n="11" /> At this stage we also introduce the data variable $N_{\text{obs}}$ that is the number of observed dwarf galaxies in the considered magnitude range and volume.<fnref n="12" /> For simplicity, we define $\mathcal{D}' = \mathcal{D} \cup \{N_{\text{obs}}\}$ and $\mathcal{P}' = \mathcal{P} \cup \{N_{\text{true}}\}$.

To constrain $N_{\text{true}}$ using, in particular, $N_{\text{obs}}$, we add another layer to the statistical framework and now consider the theoretical density function, unaffected by the survey footprint and detection limits, $\rho_{\text{true}}(\mathcal{D}|\mathcal{P}') = N_{\text{true}}P(\mathcal{D}|\mathcal{P})$. Folding in the detections limits yields the observed version of this function, $\rho_{\text{obs,th}}$, simply defined as

$$\begin{aligned} \rho_{\text{obs,th}}(\mathcal{D}|\mathcal{P}') &= \tau(\mathcal{D}) \rho_{\text{true}}(\mathcal{D}|\mathcal{P}') \\ &= \tau(\mathcal{D}) N_{\text{true}} P(\mathcal{D}|\mathcal{P}), \end{aligned} \qquad (9)$$

with $\tau(\mathcal{D})$ the probability of detecting a dwarf galaxy depending on its properties $\mathcal{D}$ (Doliva-Dolinsky et al. 2022).

Using the formalism of Kepner et al. (1999), Rykoff et al. (2012) and Drlica-Wagner et al. (2020), $N_{\text{true}}$ can be constrained by first virtually binning the data space. In any bin $i$, the likelihood $\ell_i(\mathcal{D}'|\mathcal{P}')$ of generating a sample of $N_{\text{obs},i}$ dwarf galaxy in bin $i$ can be described by the Poisson likelihood $\mathfrak{P}(N_{\text{obs},i}|N_{\text{obs,th},i})$. Here, the expectation $N_{\text{obs,th},i}$ is the theoretically observed number of dwarf galaxies in bin $i$, or

$$N_{\text{obs,th},i} = \tau(\mathcal{D}_i) N_{\text{true}} P(\mathcal{D}_i|\mathcal{P}) d\mathcal{D}. \qquad (10)$$

The total likelihood of the dwarf galaxy system can therefore be expressed as

$$\begin{aligned} \mathcal{L}(\mathcal{D}'|\mathcal{P}') &= \prod_{i \in \text{bins}} \ell_i(\mathcal{D}'|\mathcal{P}') \\ &= \prod_{i \in \text{bins}} \mathfrak{P}(N_{\text{obs},i}|N_{\text{obs,th},i}) \\ &= \prod_{i \in \text{bins}} (N_{\text{obs,th},i})^{N_{\text{obs},i}} \exp(-N_{\text{obs,th},i})/N_{\text{obs},i}!, \end{aligned} \qquad (11)$$

which further translates into

$$\log(\mathcal{L}(\mathcal{D}'|\mathcal{P}')) = - \sum_{i \in \text{bins}} N_{\text{obs,th},i} + \sum_{i \in \text{bins}} N_{\text{obs},i} \log(N_{\text{obs,th},i}) \\ + \text{const.} \qquad (12)$$

The first term of this equation is simply the integral of $\rho_{\text{obs,th}}$ over the data space (i.e., the number of observable satellites predicted by the model). In addition, if we consider bins that are small enough to contain either one or no galaxy, the second part of the equation then becomes a sum over the $N_{\text{obs}}$ bins that contain a galaxy. Equation (12) therefore becomes

$$\begin{aligned} \log(\mathcal{L}(\mathcal{D}'|\mathcal{P}')) &= - \int_{\mathcal{D}} \rho_{\text{obs,th}} d\mathcal{D} + \sum_{i=1}^{N_{\text{obs}}} \log(N_{\text{obs,th},i}) + \text{const} \\ &= - \int_{\mathcal{D}} \tau(\mathcal{D}) N_{\text{true}} P(\mathcal{D}|\mathcal{P}) d\mathcal{D} \\ &\quad + \sum_{i=1}^{N_{\text{obs}}} \log(\tau(\mathcal{D}_i) N_{\text{true}} P(\mathcal{D}_i|\mathcal{P})) + \text{const.} \end{aligned} \qquad (13)$$

Here, $\mathcal{D}_i$ are the Data Values of Dwarf Galaxy $i$.

With the assumed isotropic distribution of the satellites, the integral of Equation (13) can be marginalized over $\theta$ (in the 2D case) or $\theta$ and $\phi$ (in the 3D case), which introduces the mean fraction of detected dwarf galaxies at radius $r_{\text{M31}}$, $\langle \tau(M_V, \log r_h, r_{\text{M31}}) \rangle$. This allows us to drop the dependence

<footnote n="10">¹⁰ It may be tempting to assume models informed by the distribution of dark matter subhalos in simulations, such as, for example, a Navarro–Frenk–White profile (Navarro et al. 1996). However, considering that the region with $r_{\text{M31}} < 30$ kpc is masked in our study, the concentration of the profile would be difficult to constrain. We will further explore a more complex modeling of the radial density distribution function in a future contribution.</footnote>

<footnote n="11">¹¹ In all that follows, we choose the magnitude range $-17 < M_V < -5.5$ that, at the bright end, includes the brightest M31 dwarf that is in the survey footprint (NGC 147) and, at the faint end, is fainter than the faintest dwarf galaxy known around M31 (And XXVI; $M_V = -6^{+0.7}_{-0.5}$). The volume we consider is delimited by 30 kpc $< r_{\text{M31}} < 300$ kpc, bound by a rough estimate of the virial radius of M31 and an inner boundary that corresponds to a region in which the search for dwarf galaxy is made extremely difficult by the presence of the galaxy’s disk (Doliva-Dolinsky et al. 2022). We also choose to explore a size range of $1.8 < \log r_h < 4$ which encompass the size of all known M31 dwarf galaxies (Table 1).</footnote>

<footnote n="12">¹² While we consider a sample constructed for the 24 dwarf galaxies listed in Table 1, $N_{\text{obs}}$ may not always be 24 as our drawing from the uncertainties on the parameters of the dwarf galaxies and, in particular, their distance may, in a small number of cases, push a sample dwarf galaxy outside of the studied volume.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 2.** PDFs of the inferred number of satellites with $M_V < -5.5$ within 300 kpc of M31 and within the PAndAS footprint for the 3D and 2D cases, with and without AndXXVI. We can note that the 3D models predict ~10 more satellites to be found in PAndAS with $M_V < -5.5$, while the 2D model predicts ~50.</figcaption></figure>

<table><caption>**Table 2**
Values for the Model Parameters in the Case of a 2D and 3D Spatial Distribution of Dwarf Galaxies</caption><thead><tr><th></th><th>β</th><th>z<sub>p</sub></th><th>s</th><th>σ</th><th>α</th><th>N<sub>true</sub></th></tr></thead><tbody><tr><td>2D</td><td>-1.5 ± 0.1</td><td>2.5<sup>+0.2</sup><sub>-0.1</sub></td><td>-0.05<sup>+0.03</sup><sub>-0.02</sub></td><td>0.32<sup>+0.07</sup><sub>-0.05</sub></td><td>-0.1<sup>+0.3</sup><sub>-0.5</sub></td><td>136<sup>+65</sup><sub>-35</sub></td></tr><tr><td>3D</td><td>-1.5 ± 0.1</td><td>2.5<sup>+0.2</sup><sub>-0.1</sub></td><td>-0.05<sup>+0.03</sup><sub>-0.02</sub></td><td>0.33 ± 0.06</td><td>-1.7<sup>+0.4</sup><sub>-0.3</sub></td><td>92<sup>+19</sup><sub>-26</sub></td></tr><tr><td>2D<sub>med.-68%</sub></td><td>-1.5 ± 0.1</td><td>2.6<sup>+0.2</sup><sub>-0.1</sub></td><td>-0.05 ± 0.03</td><td>0.36<sup>+0.1</sup><sub>-0.07</sub></td><td>-0.2 ± 0.5</td><td>172<sup>+90</sup><sub>-55</sub></td></tr><tr><td>3D<sub>med.-68%</sub></td><td>-1.5 ± 0.1</td><td>2.6<sup>+0.2</sup><sub>-0.1</sub></td><td>-0.05 ± 0.03</td><td>0.36<sup>+0.1</sup><sub>-0.07</sub></td><td>-1.6 ± 0.4</td><td>95<sup>+35</sup><sub>-25</sub></td></tr></tbody></table>

**Notes.** This table presents the inferred global properties of the M31 satellites system in the case of a 2D and 3D spatial distribution. As defined in Section 3, $\beta$ is the slope of the luminosity function, $z_p$ is the zero-point, $s$ is the slope, and $\sigma$ is the dispersion of the linear relation between the size and the luminosity of dwarf galaxies, $\alpha$ is the slope of the spatial distribution and $N_{\text{true}}$ is the expected number of M31 satellites. 2D<sub>med.-68%</sub> and 3D<sub>med.-68%</sub> are the results obtained if we use the median and quantiles.

on $\theta$ and $\phi$ and the likelihood finally becomes

$$ \log(\mathcal{L}(D'|P')) = -\int_{M_V} \int_{\log r_h} \int_{r_{M31}} \langle \tau(M_V, \log r_h, r_{M31}) \rangle \times N_{\text{true}} P(D|P) dr_{M31} d\log r_h dM_V + \sum_{i=1}^{N_{\text{obs}}} \log(\tau(D_i) N_{\text{true}} P(D_i|P)) + \text{const.} \quad (14) $$

### 3.3. Implementation

We sample the likelihood with our own Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970). In order to obtain the probability distribution function (PDF) for each parameter while taking into account the uncertainties on the observed properties of the dwarf galaxies, we fold in the PDFs on the observed parameters instead of using a single value for each property. Following Conn et al. (2012), the likelihood becomes the convolution of the likelihood for a single value (Equation (14)) with the PDF of each observed property of the satellite system. With $\Omega$ the sample of all possible sets of values $D'$ and $g(D')$ being the probability of a given set, the likelihood function becomes

$$ \mathcal{L}(D'_{\Omega}|P') = \int_{\Omega} \mathcal{L}(D'|P')g(D')dD'. \quad (15) $$

This integral is calculated numerically via a Monte Carlo method and the random drawing of 50 satellite systems generated from the PDFs of $m_V$, $r_h$, $D_{\text{MW}}$ and $D_{\text{M31}}$ for all galaxies in the sample. The final distribution is the sum of the resulting chains.

Finally, we have only considered likelihoods to this point, but we seek to determine the probability of the model given the data $P(P'|D'_{\Omega})$. It is linked to the probability of the data given the model $\mathcal{L}(D'_{\Omega}|P')$ via the prior $P(P')$ such that

$$ P(P'|D'_{\Omega}) \propto \mathcal{L}(D'_{\Omega}|P')P(P'). \quad (16) $$

For simplicity, we choose uniform priors on all parameters and we impose that $0 < \beta < -4$, $1 < s < -1$, $1 < \alpha < -4$, $0 < N_{\text{true}} < 1000$, $0 < \sigma < 1$, and $0 < z_p < 3$.

## 4. Results

### 4.1. Inferred Global Properties of the M31 Dwarf Galaxy System

The constraints on the global dwarf galaxy satellite system of M31 are listed in Table 2 for both models with a sky-projected (2D, $\alpha = \alpha_{2D}$) and a volumetric (3D, $\alpha = \alpha_{3D}$) radial distribution component. The marginalized, posterior PDF of the different parameters of $P'$ are presented in Figure 2 for the case of the 3D radial distribution model. The favored parameters listed in Table 2 correspond to the peak of a parameter's marginalized one-dimensional PDF and the
associated credible intervals are bound by the parameter values that have a PDF values of 0.61 of the maximum (equivalent to a ±1σ confidence interval in the case of a Gaussian PDF and that we prefer over the 68% central confidence interval in the case of skewed PDFs). In order to provide an extensive view of our results, we also provide in Table 2 the median and 68% credible interval for each of the parameters. In any case, the MCMC chains are available at https://github.com/dolivadolinsky.

We first note that all six parameters of the model are well constrained and that the posterior PDFs are rarely perfect Gaussians or Poisson distribution in the case of $N_{\text{true}}$. This is likely a consequence of the complexity of the model and the nontrivial impact of the detection limits on the model. In addition to constraints on the individual parameters, the statistical framework we have developed makes it very easy to study the correlations (or lack thereof) between different parameters. For instance, we note the strong correlation between $\beta$, the slope of the luminosity function, and $N_{\text{true}}$, the number of M31 dwarf galaxies in the considered volume and magnitude range. This correlation is expected as changes to the slope of the luminosity function will directly lead to a change in the number of dwarf galaxies predicted by the model. The correlation between the slope of the size–luminosity relation, $s$, and the value of the slope at $M_V = -6.0$, $z_p$, is expected as there is a tradeoff between making the relation flatter and higher to ensure it goes through the cloud of data points. Similarly, a relation with a higher zero-point will result in a larger scatter. However, it is interesting to note the weak correlation between the slope and the scatter that might result from the symmetry assumption on $\sigma$. This choice is further discussed in Section 4.3.1.

Focusing on this part of the model, the size–luminosity relation for the M31 dwarf galaxies is shown in Figure 3, overlaid on the data of the 24 dwarf galaxies of the sample (teal points with error bars) and the average detection limits determined by Doliva-Dolinsky et al. (2022) in this space (gray-scale background). From the marginalized one-dimensional posterior PDFs, we derive $z_p = 2.5_{-0.1}^{+0.2}$, $s = -0.05_{-0.02}^{+0.03}$, and $\sigma = 0.33 \pm 0.06$. This relation is compatible with the one obtained by Manwadkar & Kravtsov (2022) for the Milky Way satellites, even if theirs is slightly steeper. It is also similar to but shallower than the one determined by Brasseur et al. (2011), the light orange model in the figure; $z_p = 2.34 \pm 0.1$, $s = -0.09 \pm 0.02$, and $\sigma = 0.23_{-0.07}^{+0.02}$, also determined through forward modeling, but with binary detection limits (recovery fractions of 0 or 1) that follow the dotted line in Figure 3. The differences between the two favored models likely stems from these distinct detection limits and our model infers higher overall values for $r_h$ at the fainter end as it compensates for the large and faint dwarf galaxies that are yet undiscovered because of their low surface brightness limits. Even with this slight difference, the discussion of Brasseur et al. (2011) is still relevant as our result is compatible with the relation found for late type galaxies by Shen et al. (2003). This similarity may hint that the relations found for dwarf spheroidal galaxies are a continuation of the one obtained for more massive low-concentration galaxies that includes dwarf ellipticals.

The radial distribution of the assumed isotropic distribution of dwarf galaxies around M31 is well-constrained with $\alpha_{2D} = -0.1_{-0.5}^{+0.3}$ for the 2D case. It implies an almost flat surface density of dwarf galaxies on the sky and confirms

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Relation between the size and the luminosity of M31's dwarf galaxies as inferred through our modeling. The best relation and corresponding width are represented by the purple full line and dashed lines, with the uncertainties on the mean model shown as the high opacity purple band. The best model derived by Brasseur et al. (2011) is represented by the orange lines, and the binary completeness limits they used by the dotted black line. The average detection limits folded in our analysis are represented by the gray background scale (100% recovery in white and 0% recovery in dark gray). Given those, the inferred model compensates for the undiscovered large and faint dwarf galaxies and therefore is slightly shifted from what we would naïvely expect from the cloud of known dwarf galaxies (teal dots).</figcaption></figure>

previous hints that this is the case (McConnachie et al. 2009). For the 3D case, we determine a slope $\alpha_{3D} = -1.7_{-0.3}^{+0.4}$. A previous constraint on the slope of the 3D radial distribution function was provided by Conn et al. (2012), who determined $\alpha_{3D} = -1.52_{-0.35}^{+0.32}$ from the forward modeling of a very similar sample of dwarf galaxies but different distance values for the dwarf galaxies<fnref n="13" /> and also without taking the detection limits into account. The two constraints are nevertheless compatible within their uncertainties and yield a fairly steep density profile, even though it is not as steep as the Navarro–Frenk–White profile in the external part of the halo ($\alpha \sim -3$; Navarro et al. 1996).

One of the parameters that is most affected by the detection limits is certainly the slope of the luminosity function as dwarf galaxies that are missed because they are too faint to be discovered in PAndAS will cause the observed luminosity function to drop significantly at faint magnitudes. Our analysis yields a strong constraint on this slope and we infer that the intrinsic luminosity function of M31 dwarf galaxies has a slope $\beta = -1.5 \pm 0.1$ (in the 3D case, similar in the 2D case). It is steeper than the one derived by Crnojević et al. (2019,

<footnote n="13">¹³ We now use the updated RR Lyrae distances determined by Savino et al. (2022) instead of distances of the tip of the red giant branch determined by Conn et al. (2012).</footnote>
$\beta = -1.22^{+0.11}_{-0.10}$), who did not model the impact of the detection limits. This difference arises from the consideration of the detection limits of the survey and therefore highlights their importance in deriving the faint end luminosity function.

Finally, we infer the total number of M31 dwarf galaxies within the considered magnitude range ($-17.0 < M_V < -5.5$) and volume ($30 < r_{M31/\,kpc} < 300$), $N_{\text{true}} = 136^{+65}_{-35}$ in the 2D case, or $N_{\text{true}} = 92^{+19}_{-26}$ in the 3D case. Combined with our framework that only considers galaxies brighter than $M_V = -5.5$ ($\sim 10^4 L_\odot$), and even though the realm of dwarf galaxies extends to much fainter systems (e.g., around the MW; McConnachie 2012) that are not detectable in PAndAS (Martin et al. 2013a), these values are in line with the expectation that a galaxy like M31 is surrounded by hundreds of dwarf galaxies, most of them faint (e.g., Garrison-Kimmel et al. 2019).

To check the quality of the model inference, a comparison of the favored model with the cumulative distribution of observed dwarf galaxies is shown in Figure 4 for the luminosity function and in Figure 5 for the 3D radial distribution. The inferred model is shown in black, with the gray band tracking the corresponding uncertainties, while the dashed lines represent the model, corrected to include only the PAndAS footprint, and the dashed-dotted line further adds the impact of the detection limits. This final line is directly comparable with the observations (teal points) and shows a good agreement for both cumulative distributions. This is the sign that, despite its complexity, the favored model is a good representation of the known population of M31 dwarf galaxies. These figures also make it evident that the majority of still undiscovered M31 dwarf galaxies brighter than $M_V = -5.5$ are located outside of the PAndAS footprint, mainly beyond 150 kpc (the difference between the full and dashed line) but that about half of the dwarf galaxies in the magnitude range $-5.5 > M_V > -7.0$ remain to be discovered in the PAndAS footprint (difference between the dashed and dashed-dotted line). Some of these are likely to be among the list of candidate satellites already published (Martin et al. 2013a; Mackey et al. 2019).

## 4.2. Anisotropy in the Satellite Distribution

With the inferred properties of the isotropic model we have constructed to represent the M31 dwarf galaxy system, we can now explore the perceived anisotropy of the satellite system and, in particular, whether it could be an artifact produced by the detection limits of the PAndAS survey. Looking at the RR Lyrae distances obtained by Savino et al. (2022), overlaid on the average detection limits in Figure 6, it is clear that the distribution of Andromeda's dwarf galaxies is not isotropic. Among the 24 dwarf galaxies present in the PAndAS footprint and that contribute to our sample, 21 systems are located on the MW side of M31 and produce a strong anisotropy. The (in) completeness impacts the distribution of known satellites in two different ways: it is easier to detect a dwarf galaxy that has a smaller heliocentric distance but, because of the increasing foreground contamination, it is harder to detect a dwarf galaxy closer to the MW plane. Therefore, we aim to test if the favored inferred model, observed through the detection limits (the contours in Figure 6) could naturally produce this observed anisotropy.

To quantify the significance of the anisotropy, we use a simple Monte Carlo procedure to generate 10,000 satellite systems drawn from the favored isotropic model, folding in the PAndAS recovery fractions. In practice, we start by drawing

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Cumulative number of dwarf galaxies as a function of their magnitude. The inferred model is represented by the black line and gray band. The dashed line shows the inference for the PAndAS footprint and the dotted-dashed line the favored model once the average detection limits are applied. This line is directly comparable to, and shows good agreement with, the cumulative distribution of known dwarf galaxies (teal points).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Same as Figure 4, but for the cumulative number of dwarf galaxies as a function of their distance to M31. The favored model, once the spatial and detection limits are accounted for, is compatible with the cumulative distribution of known dwarf galaxies.</figcaption></figure>

the distance to M31 from the PDF obtained by Savino et al. (2022). Then, we randomly locate this satellite around M31 using the favored density model and, finally, we test them against the detection limits of Doliva-Dolinsky et al. (2022). We reject dwarf galaxies that do not pass this test and repeat this procedure until the sample of “observed” dwarf galaxies contains 24 satellites.

Of these 10,000 systems drawn from the favored isotropic model, we find that only 7 systems have a distribution that is at least as anisotropic as M31's (at least 21 dwarf galaxies on the MW side of M31). Therefore, we conclude that the asymmetric dwarf galaxy completeness limits of the survey are very unlikely to explain, on their own, the observed anisotropy of the M31 dwarf galaxy satellite system.

## 4.3. Caveats

### 4.3.1. Limit of the Model

Given the small number of data points (i.e., the 24 known dwarf galaxies is our sample) and our decision to design an diagnostic and empirical framework, we chose to restrict the different components of the model to rather simple functional forms in order to avoid an increase in the number of parameters that would only be poorly constrained. Here, we propose to
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Right panel: Map of the completeness in the y–z plane where z is the distance to M31 along the line of sight to the MW. The contamination due to MW and M31 stars lead to the variation of the completeness along the y-axis, while the impact of the distance is visible along the z-axis. Left panel: Map of the observable satellite density predicted by the best-fit (isotropic) model in the y–z plane. The color scale includes the masked region within 30 kpc of M31. The red dots represent the known dwarf galaxies. The model predicts a slightly higher density of observed satellites on the side closer to the MW (+z) but this impact is compensated by the increase of the observed space on the negative end.</figcaption></figure>

review the limitations of those choices and the different biases
that they may potentially produce.

First, reionization is expected to affect galaxies as bright as $M_V = -7.0$ (Koposov et al. 2009; Brown et al. 2014; Weisz et al. 2014). As we choose to model the luminosity function as a simple power law and we extrapolate our model until $M_V = -5.5$, it might result in a slightly overestimation of the number of dwarf galaxies expected around M31.

Then, simulations predict a change in the slope of the radial
density profile between the most central regions, where it is
shallower, and the outer parts of the halo, where it is steeper
(Springel et al. 2008; D’Onghia et al. 2010; Garrison-Kimmel
et al. 2017; Sawala et al. 2017; Samuel et al. 2020). This
change is driven by the stronger impact of tides closer to the
host. For the model, we used a simple, single power law
throughout the considered range between 30 and 300 kpc from
M31. By masking the inner 30 kpc, we should not be biased by
the strongest tidal effects but the simplicity of our model could
be responsible for the small overprediction of dwarf galaxies in
the region within 100 kpc displayed in Figure 5 (although we
note that this apparent discrepancy is not very significant).
Including the impact of tidal effects is, at this stage, very
difficult as the full orbital information is not available for the
vast majority of M31 satellites and tidal stripping is related to
both the orbit and the mass/density of the satellites (Springel
et al. 2008). In the outer parts of the halo, while imposing a
hard cutoff at 300 kpc could lead to an overestimation of the
number of satellites between 200 and 300 kpc, it may be
compensated for by the exclusion of dwarf galaxies that are
slightly further away.

Finally, we could design a model where the dispersion of the
size–luminosity relation changes as a function of the magnitude
but, given the small number of known dwarf galaxies, it would
be difficult to constrain such a parameter without adding strong
priors build, for instance, from simulations. We have decided to
be keep our analysis data-driven so our results can be more
appropriately and independently compared with the results of
simulations.

Although the simplicity of our model does not muster all the
subtleties of the reality for the faint end of dwarf galaxies, we

do not anticipate our choices will have a major effect on our
results given the small sample of known M31 dwarf galaxies.
However, searching for and discovering the expected 92<sup>+19</sup><sub>-26</sub>
satellites would open the possibility to complicate the model
and take into account the effects mentioned above.

4.3.2. AndXXVI

While AndXXVI is part of our sample, its significance
of detection is just below the threshold used to obtain the
recovery fraction (i.e., it has a significance of 5.9 with a
threshold of 6). Considering that a threshold of 5.9 or 6 should
not lead to a significant change in the detection limits, given the
uncertainties of the recovery fraction, we decide to include
AndXXVI in our sample. However, retrospectively, it may
have been preferable to derive the detection limits with a
threshold ≪5.9. Nevertheless, it is important to note that its
low significance of detection is reflected in τ(*D*), the
probability of detecting a dwarf galaxy with properties *D*. To
quantify the impact of its inclusion on our results, we also
inferred the parameters when excluding this satellite for the
3D case. We find *β* = −1.5 ± 0.1, *z*<sub>p</sub> = 2.5<sup>+0.3</sup><sub>−0.2</sub>, *s* = −0.1<sup>+0.07</sup><sub>−0.06</sub>,
*σ* = 0.34 ± 0.08, *α* = −1.8<sup>+0.5</sup><sub>−0.4</sub>, and *N*<sub>true</sub> = 75<sup>+23</sup><sub>−26</sub>. Although
there are some variations in the parameter values, they are
statistically compatible with the parameters obtained by
including AndXXVI.

4.3.3. Impact of the Anisotropy

It is interesting to note that α<sub>3D</sub> is not quite α<sub>2D</sub> − 1,
suggesting that the surface density model, according to our
definitions of these two parameters, is not simply the
integration of the 3D density model, and therefore the chosen
models are not a perfect representation of the data. A plausible
cause could be that the assumption of isotropy may not be
entirely valid in the 3D case. Indeed, as seen in Section 4.2, the
asymmetry in the spatial distribution of M31 satellites is a real
feature, as it is unlikely to be explained by detection limits
alone. Therefore, imposing the 3D spatial distribution to be
symmetric could introduce a bias in the determination of α<sub>3D</sub>.
## 5. Discussion and Conclusion

In this paper we inferred the global properties of the M31 dwarf galaxy satellite system within the PAndAS survey through forward modeling, carefully folding in the detection limits of dwarf galaxies in this survey. We conclude that M31 hosts $92_{-26}^{+19}$ or $136_{-35}^{+65}$ satellites, depending on whether the spatial distribution of satellites is modeled in 2D or 3D, over the magnitude range $-17 < M_V < -5.5$ and $30 < r_{M31/kpc} < 300$ kpc, even though only 24 are known within the footprint and 33 overall. Both results would be slightly overestimated if the true luminosity function of M31 satellite dwarf galaxies were to differ from a power law at fainter magnitudes than $M_V = -7$, for instance because of reionization (Brown et al. 2014; Weisz et al. 2014). In order to better constrain the luminosity function at the faint end and to model the possible turn over, it is essential to increase the size of the faint dwarf galaxy sample. As shown in Figure 4, about half of the dwarf galaxies remain undiscovered in the PAndAS footprint and for the chosen magnitude limit. Yet, the most promising regions to search for these still fairly bright but elusive galaxies is to search for them outside the survey footprint: about half of the expected tally of dwarf galaxies brighter than $M_V = -5.5$ reside beyond the edges of the survey, out to the projected virial radius of M31. This also shows the importance of probing a wider area of the M31 surroundings to better constrain its global satellite properties. For the MW, model predictions range from ~20 to ~40 satellites with $M_V > -5.5$ and $r_{MW} < 300$ kpc (Jethwa et al. 2018; Newton et al. 2018; Nadler et al. 2020; Manwadkar & Kravtsov 2022) and Drlica-Wagner et al. (2020) derived an empirical estimates of ~30 satellites. For M31 we found $92_{-26}^{+19}$, which means that it host ~2–3 times more satellites than the MW for the same luminosity range. Those results are in agreement with the larger number of M31 globular clusters compared to the MW (Huxor et al. 2014) and could be easily explained if M31 is significantly more massive than the MW.

Assuming a linear relation between $\log r_h$ and $M_V$, we infer a slope $s = -0.05_{-0.02}^{+0.03}$, a zero-point $z_p = 2.5_{-0.1}^{+0.2}$ at magnitude $M_V = -6.0$, and a scatter $\sigma = 0.33 \pm 0.06$. These are consistent with the values derived by Brasseur et al. (2011; $z_p = 2.38_{-0.13}^{+0.16}$, $s = -0.03 \pm 0.03$, and $\sigma = 0.2_{-0.02}^{+0.08}$) for a similar model for the dwarf galaxies of M31 and, also of the MW. We therefore conclude that the satellite population of M31 and the MW do not show any significant difference in their size–luminosity relation.

We assumed that the slope of the luminosity function could be represented by a power law, which we constrained to be $\beta = -1.5 \pm 0.1$. This result is compatible with previous studies and for other satellite systems. It is, for instance, coherent with the results based on the Galaxy and Mass Assembly survey: Loveday et al. (2012) derive a slope of $-1.50 \pm 0.03$ for the faint end of the luminosity function of blue galaxies. For the MW, Tollerud et al. (2008) found $\beta = -1.9 \pm 0.2$ while Koposov et al. (2008) derived $\beta \sim -1.25$ and, for M81, Chiboucas et al. (2013) derived a $\beta = -1.27 \pm 0.04$. While some of those differences may be due to the merger history of the host galaxy (about 10 Gyr for the MW,<fnref n="14" /> 2–4 Gyr for M31, and ongoing for M81; e.g., D’Souza & Bell 2018; Helmi et al.

2018; Okamoto et al. 2019), it is important to note that these constraints on $\beta$ were obtained from very different sample sizes and techniques and without always taking a survey’s detection limits into account. In the case of M81, for instance, Chiboucas et al. (2013), were limited to only five galaxies with $M_V > -10$, among which only one was fainter than $M_V = -8.0$. The same authors also derived $\beta = -1.13 \pm 0.06$ for M31 Chiboucas et al. (2009), which is significantly different from our inference, albeit from an earlier sample of known M31 dwarf galaxies and without taking the survey footprint and the detection limits into account.

This comparison underlines the importance of the impact of the detection limits on constraining the global properties of a satellite system. Determining these limits can admittedly be quite a tedious and computationally taxing task (Koposov et al. 2008; Drlica-Wagner et al. 2020; Doliva-Dolinsky et al. 2022), but it is absolutely essential to fold them into the analysis. The simplicity of using binary detection limits with magnitude cuts (e.g., Brasseur et al. 2011; Bennet et al. 2019; Crnojević et al. 2019) can be appealing but requires a very conservative approach (and throwing away some of the data) to not induce a bias on the slope of the luminosity function. We explore this effect in Figure 7 as we simplify the Doliva-Dolinsky et al. (2022) detection limits to binary limits with a magnitude cut at $M_V = -6.0$. This leads to a clear bias on the slope of the luminosity function as undetected galaxies at the faint end are not compensated for. In fact, the value of $\beta$ inferred in this case is similar to those determined by Chiboucas et al. (2009) and Crnojević et al. (2019) without taking the detection limits into account. The importance of determining accurate detection limits is complementary to the use of Bayesian inference. Forward modeling has been used to go from the luminosity function to the stellar-halo mass relation and vice versa (Danieli et al. 2022), but the most common method to derive the luminosity function from observations is through their noisy correction from the incompleteness and, in some case, the unadvised fitting of the cumulative function (Chiboucas et al. 2009; Geha et al. 2017; Bennet et al. 2019; Crnojević et al. 2019). While this avoids the complex work of forward modeling, biases may arise from the loss of information implied by the use of a correction factor (e.g., an average of the detection limits, a magnitude cut). Additionally, while fitting the cumulative function may seem like a straightforward process, those data points are intrinsically correlated, which makes it difficult to properly handle statistics and, consequently, uncertainties. The possible correlation between the slope of the luminosity function and the age of the last merger suffered by the host should therefore be consolidated by redetermining the slopes in a homogeneous analysis that would take the completeness of the different surveys into, for example by following the path traced in this paper.

Finally, we show that the radial distribution of the satellites can be modeled by a power law of exponent $\alpha_{2D} = -0.1_{-0.5}^{+0.3}$ and $\alpha_{3D} = -1.7_{-0.3}^{+0.4}$ in the sky-projected and 3D radial distributions, respectively. However, we determined that the observed anisotropy of the satellite dwarf galaxies (McConnachie & Irwin 2006; Savino et al. 2022) is unlikely to be the consequence of the detection limits that make it more difficult to detect dwarf galaxies on the more distant side of M31. The observed anisotropy is rarely reproduced by drawing random satellite systems from the favored satellite system model that fold in the limits: only 0.07% of these have an anisotropy that

<footnote n="14">¹⁴ While a merger with the LMC satellites and the MW is ongoing (Battaglia et al. 2022), those dwarf galaxies were not taken into account by Tollerud et al. (2008) who studied the satellites in the northern hemisphere.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Marginalized probability distribution function of the parameter $\beta$, the slope of the luminosity function, in a case where the dwarf galaxy detection limits are taken into account (black PDF) or they are approximated by a simple magnitude cut (teal PDF).</figcaption></figure>

is at least as strong as the observations. As the incompleteness of the sample is ruled out as a major factor responsible for the observed anisotropy, another observational effect could be a bias in the distances. However, the fact that the anisotropy manifests using both RR Lyrae and tip of the red giant branch distances disfavor such biases (Conn et al. 2012; Savino et al. 2022). In other words, we have, by using robust distances, combined with our detailed modeling of the (in)completeness of the data, built a strong case for the anisotropy to be a real physical configuration. Different solutions are proposed to explain such a distribution (Pawlowski et al. 2017; Thomas et al. 2018; Wan et al. 2020) and complexity could be added to the model presented here to explore these possibilities.

More generically, the statistical framework developed here is very flexible and can easily be expanded. One can, for instance, imagine replacing the current isotropic spatial distribution component of the model to make it more intricate and parameterize the observed anisotropy of the system. It would also be straightforward to add more components to Equation (1), for instance to use the M31 dwarf galaxies to place constraints on the dark matter particle (e.g., Kim et al. 2018) or the faint end of galaxy formation (e.g., Koposov et al. 2009). Considering that, according to our inference, two thirds of dwarf galaxies with $M_V < -5.5$ are still undiscovered around M31, applying this framework to a wider and deeper M31 satellite sample (e.g., thanks to the surveys that will be conducted with the Euclid, Rubin, or Roman telescopes), would yield important and tighter constraints on the global properties of this system and, possibly, on the dark matter properties. Finally, the method we describe could similarly be applied to other galaxies, either the MW itself (Drlica-Wagner et al. 2020), or galaxies outside the Local Group (Mao et al. 2021; Mutlu-Pakdil et al. 2021; Carlsten et al. 2022; Nashimoto et al. 2022). Reproducing the current analysis on the dwarf galaxy systems of other hosts will be essential to place robust constraints on cosmological and galaxy formation physics and to continue refining the use of dwarf galaxies as cosmological probes. This framework will be particularly useful considering the future arrival of the Large Synoptic Survey Telescope and Euclid that are going to revolutionize the study of satellite systems beyond the Local Group.

## Acknowledgments

We thank Denija Crnojević and Erik Tollerud for stimulating discussions about the work presented here. We are grateful for the referee for his useful comments and suggestions that allowed us to make a more comprehensive and complete paper.

Based on observations obtained with MegaPrime/Mega-Cam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii.

The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data.

D.W. and A.S. were supported for this work by NASA through grants GO-13768, GO-15746, GO-15902, AR-16159, and GO-16273 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555.

G.F.T. acknowledge support from the Agencia Estatal de Investigación (AEI) under grant Ayudas a centros de excelencia Severo Ochoa convocatoria 2019 with reference CEX2019-000920-S, and from the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (AEI-MCINN) under grant En la frontera de la arqueología galacáctica: evolución de la materia luminoso y obscura de la vía Láctea y las galaxias enenas del Grupo Local with reference PID2020-118778GB-I00.

R.I. and N.M. acknowledge funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 834148).

## ORCID iDs

Amandine Doliva-Dolinsky https://orcid.org/0000-0001-9775-9029

Nicolas F. Martin https://orcid.org/0000-0002-1349-202X

Zhen Yuan (袁珍) https://orcid.org/0000-0002-8129-5415

Alessandro Savino https://orcid.org/0000-0002-1445-4877

Daniel R. Weisz https://orcid.org/0000-0002-6442-6030

Annette M. N. Ferguson https://orcid.org/0000-0001-7934-1278

Rodrigo A. Ibata https://orcid.org/0000-0002-3292-9709

Geraint F. Lewis https://orcid.org/0000-0003-3081-9319

Alan W. McConnachie https://orcid.org/0000-0003-4666-6564

Guillaume F. Thomas https://orcid.org/0000-0002-2468-5521

## References

Abazajian, K., Adelman-McCarthy, J., Agueros, M., et al. 2003, AJ, 126, 2081

Battaglia, G., Taibi, S., Thomas, G. F., & Fritz, T. K. 2022, A&A, 657, A54

Bell, E. F., Slater, C. T., & Martin, N. F. 2011, ApJ, 742, L15

Bennet, P., Sand, D. J., Crnojević, D., et al. 2019, ApJ, 885, 153

Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93

Boylan-Kolchin, M., Weisz, D. R., Johnson, B. D., et al. 2015, MNRAS, 453, 1503
Brasseur, C. M., Martin, N. F., Macciò, A. V., Rix, H.-W., & Kang, X. 2011, ApJ, 743, 179

Brown, T. M., Tumlinson, J., Geha, M., et al. 2014, ApJ, 796, 91

Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2000, ApJ, 539, 517

Carlsten, S. G., Greene, J. E., Beaton, R. L., Danieli, S., & Greco, J. P. 2022, ApJ, 933, 47

Chiboucas, K., Jacobs, B. A., Tully, R. B., & Karachentsev, I. D. 2013, AJ, 146, 126

Chiboucas, K., Karachentsev, I. D., & Tully, R. B. 2009, AJ, 137, 3009

Collins, M. L. M., Charles, E. J. E., Martínez-Delgado, D., et al. 2022, MNRAS, 515, L72

Conn, A. R., Ibata, R. A., Lewis, G. F., et al. 2012, ApJ, 758, 11

Conn, A. R., Lewis, G. F., Ibata, R. A., et al. 2013, ApJ, 766, 120

Crnojević, D., Ferguson, A., Irwin, M., et al. 2014, MNRAS, 445, 3862

Crnojević, D., Sand, D. J., Bennet, P., et al. 2019, ApJ, 872, 80

Danieli, S., Greene, J. E., Carlsten, S., et al. 2022, arXiv:2210.14233

D'Onghia, E., Springel, V., Hernquist, L., & Keres, D. 2010, ApJ, 709, 1138

D'Souza, R., & Bell, E. F. 2018, NatAs, 2, 737

Doliva-Dolinsky, A., Martin, N. F., Thomas, G. F., et al. 2022, ApJ, 933, 135

Dooley, G. A., Peter, A. H. G., Yang, T., et al. 2017, MNRAS, 471, 4894

Drlica-Wagner, A., Bechtol, K., Mau, S., et al. 2020, ApJ, 893, 47

Evans, T. A., Fattahi, A., Deason, A. J., & Frenk, C. S. 2020, MNRAS, 497, 4311

Garrison-Kimmel, S., Wetzel, A., Bullock, J. S., et al. 2017, MNRAS, 471, 1709

Garrison-Kimmel, S., Hopkins, P., Wetzel, A., et al. 2019, MNRAS, 487, 1380

Geha, M., Wechsler, R., Mao, Y.-Y., et al. 2017, ApJ, 847, 4

Hastings, W. K. 1970, Biometrika, 57, 97

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85

Herschel, W. 1789, RSPT, 79, 212

Huxor, A. P., Mackey, A. D., Ferguson, A. M. N., et al. 2014, MNRAS, 442, 2165

Ibata, R., Martin, N. F., Irwin, M., et al. 2007, ApJ, 671, 1591

Ibata, R. A., Lewis, G., Conn, A., et al. 2013, Natur, 493, 62

Ibata, R. A., Lewis, G., McConnachie, A., et al. 2014, ApJ, 780, 128

Irwin, M. J., Ferguson, A. M. N., Huxor, A. P., et al. 2008, ApJL, 676, L17

Jethwa, P., Erkal, D., & Belokurov, V. 2018, MNRAS, 473, 2060

Karachentsev, I. D., & Karachentseva, V. E. 1999, A&A, 341, 355

Kelley, T., Bullock, J. S., Garrison-Kimmel, S., et al. 2019, MNRAS, 487, 4409

Kepner, J., Fan, X., Bahcall, N., et al. 1999, ApJ, 517, 78

Kim, S. Y., Peter, A. H. G., & Hargis, J. R. 2018, PhRvL, 121, 211302

Koposov, S. E., Yoo, J., Rix, H.-W., et al. 2009, ApJ, 696, 2179

Koposov, S., et al. 2008, ApJ, 686, 279

Kroupa, P., Theis, C., & Boily, C. M. 2005, A&A, 431, 517

Loveday, J., Baldry, I. K., Driver, S. P., et al. 2012, MNRAS, 420, 1239

Lynden-Bell, D. 1976, MNRAS, 174, 695

Mackey, A. D., Ferguson, A., Huxor, A., et al. 2019, MNRAS, 484, 1756

Majewski, S. R., Kalirai, J., Gilbert, K., et al. 2004, AAS Meeting Abstracts, 205, 1411.13

Manwadkar, V., & Kravtsov, A. V. 2022, MNRAS, 516, 3944

Mao, Y.-Y., Geha, M., Wechsler, R. H., et al. 2021, ApJ, 907, 85

Martin, N. F., Ibata, R. A., Irwin, M. J., et al. 2006, MNRAS, 371, 1983

Martin, N. F., Ibata, R. A., McConnachie, A. W., et al. 2013a, ApJ, 776, 80

Martin, N. F., McConnachie, A., Irwin, M., et al. 2009, ApJ, 705, 758

Martin, N. F., Slater, C., Schlafly, E., et al. 2013b, ApJL, 779, L10

Martin, N. F., Schlafly, E., Slater, C., et al. 2013c, ApJ, 772, 15

Martin, N. F., Ibata, R., Lewis, G., et al. 2016, ApJ, 833, 167

Martin, N. F., Weisz, D., Albers, S., et al. 2017, ApJ, 850, 16

Martínez-Delgado, D., Karim, N., Charles, E. J. E., et al. 2022, MNRAS, 509, 16

Mashchenko, S., Wadsley, J., & Couchman, H. M. P. 2008, Sci, 319, 174

McConnachie, A. W. 2012, AJ, 144, 4

McConnachie, A. W., & Irwin, M. J. 2006, MNRAS, 365, 902

McConnachie, A. W., Huxor, A., Martin, N. F., et al. 2008, ApJ, 688, 1009

McConnachie, A. W., Irwin, M., Ibata, R., et al. 2009, Natur, 461, 66

McConnachie, A. W., Ibata, R., Martin, N., et al. 2018, ApJ, 868, 55

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, JChPh, 21, 1087

Metz, M., Kroupa, P., & Jerjen, H. 2007, MNRAS, 374, 1125

Müller, O., Pawlowski, M. S., Jerjen, H., & Lelli, F. 2018, Sci, 359, 534

Mutlu-Pakdil, B., Sand, D., Crnojević, D., et al. 2021, ApJ, 918, 88

Nadler, E. O., Wechsler, R. H., Bechtol, K., et al. 2020, ApJ, 893, 48

Nadler, E. O., Drlica-Wagner, A., Bechtol, K., et al. 2021, PhRvL, 126, 091101

Nashimoto, M., Tanaka, M., Chiba, M., et al. 2022, ApJ, 936, 38

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Newton, O., Cautun, M., Jenkins, A., Frenk, C. S., & Helly, J. C. 2018, MNRAS, 479, 2853

Ocvirk, P., & Aubert, D. 2011, MNRAS, 417, L93

Okamoto, S., Arimoto, N., Ferguson, A. M. N., et al. 2019, ApJ, 884, 128

Pawlowski, M. S. 2018, MPLA, 33, 1830004

Pawlowski, M. S., Ibata, R. A., & Bullock, J. S. 2017, ApJ, 850, 132

Preston, J., Collins, M., Ibata, R., et al. 2019, MNRAS, 490, 2905

Richardson, J. C., Irwin, M., McConnachie, A., et al. 2011, ApJ, 732, 76

Rykoff, E. S., Koester, B. P., Rozo, E., et al. 2012, ApJ, 746, 178

Samuel, J., Wetzel, A., Tollerud, E., et al. 2020, MNRAS, 491, 1471

Savino, A., Weisz, D., Skillman, E., et al. 2022, ApJ, 938, 101

Sawala, T., Pihajoki, P., Johansson, P. H., et al. 2017, MNRAS, 467, 4383

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978

Slater, C. T., Bell, E. F., & Martin, N. F. 2011, ApJL, 742, L14

Somerville, R. S. 2002, ApJL, 572, L23

Spergel, D. N., & Steinhardt, P. J. 2000, PhRvL, 84, 3760

Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685

Thomas, G. F., Famaey, B., Ibata, R., et al. 2018, A&A, 609, A44

Tollerud, E. J., Bullock, J. S., Strigari, L. E., & Willman, B. 2008, ApJ, 688, 277

van den Bergh, S. 1972, ApJL, 171, L31

Wan, Z., Oliver, W. H., Lewis, G. F., Read, J. I., & Collins, M. L. M. 2020, MNRAS, 492, 456

Weisz, D. R., & Boylan-Kolchin, M. 2017, MNRAS, 469, L83

Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014, ApJ, 789, 148

Weisz, D. R., Martin, N., Dolphin, A., et al. 2019, ApJL, 885, L8

Wheeler, C., Oñorbe, J., Bullock, J. S., et al. 2015, MNRAS, 453, 1305

Zucker, D. B., Kniazev, A., Bell, E., et al. 2004, ApJL, 612, L121

Zucker, D. B., Kniazev, A., Martínez-Delgado, D., et al. 2007, ApJL, 659, L21