# Forward-modelling the Luminosity, Distance, and Size distributions of the
Milky Way Satellites

Viraj Manwadkar1 and Andrey V. Kravtsov1,2,3
1Department of Astronomy & Astrophysics, The University of Chicago, Chicago,
IL 60637 USA
2Kavli Institute for Cosmological Physics, The University of Chicago, Chicago,
IL 60637 USA
3Enrico Fermi Institute, The University of Chicago, Chicago, IL 60637 USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We use GRUMPY, a simple regulator-type model for dwarf galaxy formation and
evolution, to forward model the dwarf galaxy satellite population of the Milky
Way (MW) using the Caterpillar zoom-in simulation suite. We show that
luminosity and distance distributions of the model satellites are consistent
with the distributions measured in the DES, PS1 and SDSS surveys, even without
including a model for the orphan galaxies. We also show that our model for
dwarf galaxy sizes can simultaneously reproduce the observed distribution of
stellar half-mass radii, $r_{\rm 1/2}$, of the MW satellites and the overall
$r_{\rm 1/2}-M_{\star}$ relation exhibited by observed dwarf galaxies. The
model predicts that some of the observed faint stellar systems with $r_{\rm
1/2}<10$ pc are ultra-faint dwarf galaxies. Scaling of the stellar mass
$M_{\star}$ and peak halo mass $M_{\rm peak}$ for the model satellites is not
described by a power law, but has a clear flattening of $M_{\star}-M_{\rm
peak}$ scaling at $M_{\rm peak}<10^{8}\,\,M_{\odot}$ imprinted by
reionization. As a result, the fraction of low mass haloes ($M_{\rm
peak}<10^{8}\,M_{\odot}$) hosting galaxies with $M_{V}<0$ is predicted to be
50% at $M_{\rm peak}\sim 3.6\times 10^{7}\,\,M_{\odot}$. We find that such
high fraction at that halo mass helps to reproduce the number of dwarf
galaxies discovered recently in the HSC-SSP survey. Using the model we
forecast that there should be the total of $440^{+201}_{-147}$ (68% confidence
interval) MW satellites with $M_{\rm V}<0$ and $r_{\rm 1/2}>10$ pc within 300
kpc and make specific predictions for the HSC-SSP, DELVE-WIDE and LSST
surveys.

###### keywords:

galaxies:dwarf, galaxies: haloes, galaxy: formation, galaxy: structure,
galaxies: luminosity function

††pubyear: 2021††pagerange: Forward-modelling the Luminosity, Distance, and
Size distributions of the Milky Way Satellites–5

##  1 Introduction

Most galaxies in the Universe are dwarf galaxies with luminosities of
$L\lesssim 10^{9}\,L_{\odot}$. Due to their shallow potential wells, dwarf
galaxies are particularly sensitive to feedback processes (e.g., Dekel & Silk,
1986; Efstathiou, 2000; Pontzen & Governato, 2012), such as heating of
intergalactic gas during reionization (e.g., Efstathiou, 1992; Thoul &
Weinberg, 1996; Barkana & Loeb, 1999; Bullock et al., 2000; Gnedin, 2000;
Somerville, 2002; Benson et al., 2002a) and galactic outflows (e.g., Hopkins
et al., 2014; Muratov et al., 2015; Pandya et al., 2021). Dwarf galaxies are
thus a key laboratory for probing and understanding these feedback processes.

In addition, the faintest of dwarf galaxies—the ultra-faint dwarf galaxies
(UFDs, $M_{\rm V}\gtrsim-8$ or $L\lesssim 10^{5}\,L_{\odot}$, see Simon, 2019,
for a review)–predicted by theoretical models (Ricotti & Gnedin, 2005; Gnedin
& Kravtsov, 2006; Bovill & Ricotti, 2009; Salvadori & Ferrara, 2009), showed
that such galaxies are expected to be the fossils of the earliest stages of
galaxy formation. These predictions were born out by observations so far,
which showed that UFDs have formed most of their stars very early (Tolstoy et
al., 2009; Brown et al., 2012, 2014; Weisz et al., 2014; Sacchi et al., 2021)
and their stellar populations are very metal poor (see Simon, 2019, for a
review). Due to the broad cosmological significance of dwarf galaxies,
observational searches for faintest and lowest surface brightness galaxies
have been the subject of much observational effort over the past two decades
(e.g., Belokurov et al., 2010; Martin et al., 2015; Longeard et al., 2018; Mau
et al., 2020; Cerny et al., 2021b; Mutlu-Pakdil et al., 2022).

Prior to 2000, discoveries of dwarf galaxies were mainly based on searches
using photographic plates (e.g., Harrington & Wilson, 1950; Wilson, 1955;
Cannon et al., 1977; Irwin et al., 1990; Ibata et al., 1994, cf. Fig. 1 in
Simon 2019). However, in the last two decades wide area digital sky surveys,
such as the Sloan Digital Sky Survey (SDSS), Dark Enery Survey (DES),
PanSTARRS (PS1), Hyper Suprime-Cam Strategic Survey Program (HSC-SSP), ATLAS,
MagLiteS, etc. have revolutionized studies of the faintest dwarf galaxies. The
number of known dwarf galaxies in the Local Group volume was tripled and
dozens of the UFDs were discovered (e.g., Willman et al., 2005; Belokurov et
al., 2007; Koposov et al., 2015; Drlica-Wagner et al., 2015a, see Simon 2019
and Drlica-Wagner et al. 2020 for recent summaries). The advent of the Vera C.
Rubin Observatory (Ivezić et al., 2019) is expected to greatly expand the
discovery space of dwarf galaxies again (e.g., Mutlu-Pakdil et al., 2021).

The low surface brightnesses and faint luminosities of UFDs make them
challenging to detect even in the most advanced modern surveys and so far
observations only probed the tip of the iceberg (e.g., Figure 5,6 in Drlica-
Wagner et al., 2020) of the UFD galaxy population expected in the Local Group
volume (e.g., Gnedin & Kravtsov, 2006). Although ongoing and future surveys
are expected to discover many more UFDs, most of the expected UFD discoveries,
especially $M_{\rm V}\gtrsim-4$, are expected to be confined to within
$\sim$300 kpc from the MW (e.g., Mutlu-Pakdil et al., 2021). Thus, for the
foreseeable future, the faintest satellites of the MW will continue to be the
window into formation of the faintest galaxies (e.g., Koposov et al., 2009;
Jethwa et al., 2018; Newton et al., 2018; Nadler et al., 2019; Nadler et al.,
2020).

The luminosity function (LF) is a fundamental property of galaxy population in
general, with the faint-end of LF in the dwarf regime of particular importance
as it is expected to be shaped by feedback processes accompanying galaxy
formation (e.g., Bullock et al., 2000; Benson et al., 2002a; Somerville, 2002;
Busha et al., 2010; Macciò et al., 2010; Font et al., 2011; Bose et al., 2018;
Rey et al., 2019; Katz et al., 2020). The UFD tail of the LF, which is
expected to be particularly sensitive to both outflows and suppression of gas
accretion due to UV heating after reionization, is particularly interesting.
In the next decade, the luminosity function of the MW satellites will remain
the most powerful observational probe of these processes in the UFD regime.

Likewise, distribution of the MW satellites around Milky Way can be used to
probe processes shaping evolution of dwarf satellite galaxies around $\sim
L_{\star}$ galaxies, such as disruption of subhaloes by the central disk
(e.g., Nadler et al., 2018; Kelley et al., 2019; Richings et al., 2020; Webb &
Bovy, 2020; Green et al., 2021), merger history of their host halo (e.g., Bose
et al., 2020) or details of reionization (e.g., Ocvirk & Aubert, 2011). In
fact, so far models have struggled to reproduce the concentration of the
radial distribution of the luminous satellites of the Milky Way and a number
of other nearby galaxies (e.g., Carlsten et al., 2020, although see Mao et al.
2021), leading to conclusions that a significant fraction of subhaloes is
missing in simulations for numerical reasons and that special treatment of the
“orphan” galaxies of these missing subhaloes is needed (Bose et al., 2020).
Indications that a surprisingly large fraction of UFDs are close to the
pericentre of their orbit (e.g., Simon, 2018; Fritz et al., 2018; Li et al.,
2021) is another puzzle to be solved.

Cosmological simulations that include all relevant physical and environmental
processes should ultimately provide the basis for our understanding of dwarf
galaxy formation (e.g., see Wheeler et al., 2015; Wheeler et al., 2019; Rey et
al., 2019, 2020; Agertz et al., 2020; Katz et al., 2020; Applebaum et al.,
2021; Jeon et al., 2021; Rey et al., 2022, for recent examples). However,
current simulations face many challenges, such as the dynamic range challenge
in resolving relevant small-scale structure both in dark matter and baryon
distribution and associated computational expense. At the same time, key
physical processes affecting galaxy evolution, such as radiative and stellar
feedback or the process of star formation in low-metallicity diffuse
environments of dwarf galaxies, are not yet well understood. All of this
demands convergence studies and parameter exploration to understand the full
range of scenarios, however, the complexity and computational expense of
simulations make this difficult in practice.

Semi-analytic models have thus been used as an alternative modelling approach
for dwarf satellite galaxies (e.g., Kauffmann et al., 1993; Somerville, 2002;
Benson et al., 2002a, b). Such approaches had good success in modelling the
observed dwarf satellites (e.g., Salvadori & Ferrara, 2009; Starkenburg et
al., 2013; Bose et al., 2018). During the past decade, approaches based on
simple phenomenological ansatzes for how galaxy properties are related to the
properties of their dark matter haloes have become widely used (e.g.,
Kravtsov, 2010; Jethwa et al., 2018; Kim et al., 2018; Nadler et al., 2020;
Jiang et al., 2021). However, the simplicity of this approach has limitation
as it does not include modelling of the full range of observed dwarf galaxy
properties, which include sizes, gas masses and metallicities.

“Regulator” type models have recently been used successfully to model a
variety of galaxy statistics and scaling relations across wide range of epochs
(see, e.g., §4.1.2 of Somerville & Davé, 2015, and references therein). These
models are simpler than classical semi-analytic models, as they aim to model
only the most crucial aspects of galaxy formation. At the same time, they
retain sufficient complexity to predict evolution of key observable properties
of galaxies, such as their gas mass, stellar and gas phase metallicities, etc.
The relative simplicity and flexibility of such models allows straightforward
exploration of the effect different parts of the model have on galaxy
properties. At the same time, the regulator models proved to be quite capable
in reproducing a wide range of observations.

Indeed, in Kravtsov & Manwadkar (2022) we presented a regulator-type model
focusing on modelling galaxies in the dwarf regime and demonstrated that the
model can simultaneously match a wide range of properties of observed nearby
dwarf galaxies and their observed correlations, including star formation
histories, stellar and gas phase metallicities, gas masses, sizes, star
formation rates. In this paper, we use results of this successful model along
with realistic survey detection probabilities to forward model observed
satellite population. Such forward modelling is particularly important for
modelling UFD satellite population and was used in a number of previous
theoretical studies of the MW satellite population based on phenomenological
models of galaxy-halo connection (e.g., Jethwa et al., 2018; Nadler et al.,
2020).

In this study we focus on forward-modelling the satellite luminosity function
down to the faintest luminosities, distribution of dwarf satellits in their
radial distance to the Milky Way, and dwarf satellite sizes probed in various
current and upcoming observational surveys. This paper is organized as
follows. In Section 2, we describe the main features of the galaxy formation
model and the forward modelling steps to produce predicted satellite
populations. In Section 3, we compare model predictions to the luminosity
function, radial distribution, and size distribution of the observed
population of MW satellites. In Section 4, we discuss implications of our
results for the census of satellite dwarf galaxies in the MW environment and
their sizes and surface brightness and then summarize our results in Section
5. In the Appendices, we discuss the impact of simulation resolution on
predicted satellite statistics (Appendix A), parametrization of relations
between peak halo mass, stellar mass and V-band luminosity (Appendix B),
sensitivity of model predictions to the current uncertainty in the MW halo
mass (Appendix C), and list the relevant properties of observed dwarf galaxies
and ultra faint star clusters used in this study (Appendix D).

##  2 The Model

To model the satellite population of the Milky Way (MW) we use the GRUMPY
dwarf galaxy model presented in Kravtsov & Manwadkar (2022). As a baseline for
modelling different galaxy components it uses halo mass accretion histories
extracted from simulations. In this section we describe halo catalogs used to
make model runs and a general overview of the model and details particularly
relevant to this study.

Figure 1: The average subhalo mass functions of subhaloes within i) 50 kpc
(dotted lines), ii) 100 kpc (dashed lines) and iii) virial radius (solid
lines) for the Caterpillar suites at the LX13 and LX14 resolution level.

###  2.1 Halo catalogs

Figure 2: (Left) The subhalo mass functions within virial radius for the two
MW-sized hosts from Caterpillar LX14 suites (denoted by colored lines) that
contain a LMC analog satellite according to criterion in Section 2.3. For
comparison, we also show the subhalo mass functions from other hosts in
Caterpillar LX14 that do not have an LMC analog (thin grey lines). (Right) The
mass assembly histories (MAHs) of the two MW-sized hosts from Caterpillar LX14
suites (denoted by colored lines) that contain a LMC analog satellite
according to criterion in Section 2.3. For comparison, we also show the MAHs
from other hosts in Caterpillar LX14 that do not have an LMC analog (thin grey
lines).

As a baseline for the galaxy formation model we use halo catalogs and halo
mass accretion histories (MAHs) derived from the Caterpillar suite (Griffen et
al., 2016) of zoom-in simulations of MW-sized
haloes.111https://www.caterpillarproject.org The Caterpillar suite consists of
simulations of 32 MW-like host haloes222A total of 35 MW-like haloes were
simulated, but a few of them are not used here as they are contaminated by low
resolution particles. Refer to Griffen et al. (2016) for more details on the
contamination study.. Each of these host haloes has been re-simulated at three
different resolution levels, namely LX12, LX13 and LX14, with LX14 being the
highest resolution of the three (see Griffen et al., 2016, for details). The
existence of several resolution levels allows checks for convergence of the
results and LX14 resolution is higher than the resolution of ELVIS simulations
used in our previous study (Kravtsov & Manwadkar, 2022). We therefore adopt
the Caterpillar suite for our modelling of dwarf satellite population.

The haloes were identified using the modified version of the Rockstar halo
finder and the Consistent Trees Code (Behroozi et al., 2013a). As discussed by
Griffen et al. (2016, see their Section 2.5), the modification is related to
the misclassification of subhaloes with large initial fraction of unbound
particles by the Rockstar, which results in a loss of $\approx 30\%$ of
massive subhaloes.

Figure 1 shows the average cumulative subhalo mass functions for the
Caterpillar suite at the LX13 and LX14 resolution levels. The mass resolution
limit of a simulation can be approximately deduced as the mass at which the
cumulative subhalo mass function begins to flatten. The figure shows that the
LX14 suite has a resolution limit of $\mu=M_{\rm peak}/M_{\rm host}\approx
4\times 10^{-6}$ ($M_{\rm peak}\approx 4\times 10^{6}\,M_{\odot}$) even in the
innermost regions of the host ($r<50$ kpc). Thus, the LX14 suite allows us to
model MW satellites with a full range of subhalo masses that are expected to
host faintest UDF galaxies probed by current observations. A detailed
convergence test via comparisons of the LX13 and LX14 simulations is presented
in the Appendix A.

We run our model using haloes with $M_{\rm peak}>10^{6}\,M_{\odot}$ and
$t_{\rm start}<4$ Gyr, where $t_{\rm start}$ is the earliest time when a given
halo is detected in the simulation. This cut is used to avoid spurious haloes
and ill-defined MAH. Only $\lesssim$1% of haloes with $M_{\rm
peak}>10^{7}\,M_{\odot}$ are removed with this cut and it will thus will not
affect our results significantly.

In our model, we do not track and model for ‘orphan haloes’, that is, dark
matter haloes which have disrupted below the detection resolution limit of
Rockstar halo finder. We discuss the rationale of this choice detail in
Section 4.4.4.

Throughout this paper, we use halo masses defined within the radius enclosing
either the “virial” density contrast (Bryan & Norman, 1998), $M_{\rm vir}$, or
the contrast of 200 times the critical density of the universe at the redshift
of analysis, $M_{\rm 200c}$. $M_{\rm vir}$ "virial" is used when computing
subhalo disruption probabilities with the Nadler et al. (2018) model and in
the studies of numerical effects on subhalo population in Section 4.4.4, while
we use $M_{\rm 200c}$ elsewhere to compute evolution of halo mass in the
galaxy formation model, $M_{\rm peak}$ for subhaloes, etc. because it is much
less affected by the pseudo-evolution effect (Diemer et al., 2013).

###  2.2 Subhalo disruption model

To account for subhalo disruption due to the extra tidal forces exerted by the
central baryonic disk, we use the random forest subhalo disruption model from
Nadler et al. (2018). This model computes a disruption/survival probability of
a subhalo given its internal and orbital properties, namely, the
proper/physical333physical/proper distance, i.e., the comoving distance
multiplied by the scale factor. distance and scale factor at the first
pericentric passage after accretion, and the virial mass, maximum circular
velocity and scale factor at time of accretion. All this information is
extracted from the evolution tracks of the simulations described in Section
2.1.

###  2.3 Selecting haloes with Large Magellanic Cloud Analog

The satellite population of the MW has been shown to have significantly
impacted by the presence of the Large Magellanic Cloud (LMC, e.g., Lynden-
Bell, 1976; Lu et al., 2016; Dooley et al., 2017; Nadler et al., 2020). The
LMC appears to be on its first infall into the MW (e.g., Besla et al., 2007;
Kallivayalil et al., 2013) and is thought to have brought with it a
significant fraction of the MW satellites (e.g., D’Onghia & Lake, 2008; Nadler
et al., 2020). It is thus important to account for the presence of the LMC-
size object along with the satellites it brought when modelling the population
of the MW satellites, especially in the surveys with the sky footprint close
to the LMC, such as the Dark Energy Survey (DES).

Thus, we focus on the two Caterpillar LX14 hosts that have a realistic LMC
analog defined as a subhalo with i) mass of $M_{\rm peak}>1\times
10^{11}\,M_{\odot}$ consistent with recent estimates of the LMC mass (e.g.,
Erkal et al., 2019; Erkal & Belokurov, 2020; Vasiliev et al., 2021; Shipp et
al., 2021), and ii) is on its first infall into the MW and iii) is at a
physical host-centric distance close to 50 kpc.

Among all the Caterpillar hosts, we did not find a single host that had an
LMC-mass subhalo on its first pericentric passage and at a host-centric
distance of $\sim 50$ kpc. However, in Caterpillar 4 and Caterpillar 53 hosts,
an LMC-mass subhalo went through a recent first pericentric passage around at
$z\sim 0.05$. For these hosts, we extract all halo track data at a previous
epoch ($z_{\rm extract}$) when LMC was on its first infall and at a host-
centric distance closest to 50 kpc. The halo properties, subhalo distribution,
and our galaxy model predictions at $z_{\rm extract}$ are then used.

Although the subhalo disruption prescription described in Section 2.2 was
calibrated at $z=0$, in the subsequent analysis we assume that disk properties
and halo structure to be similar at $z\sim 0.05$ and that the same disruption
model can be used in such hosts. Information for the two hosts with LMC
analogs is listed in Table 1. The MAHs and subhalo mass functions of these two
hosts are shown in Figure 2. Caterpillar 4 mass function is pretty close to
the median, while Caterpillar 53 mass function is on the high edge of the
distribution in amplitude, as could be expected if the LMC counterpart brought
in a significant number of subhaloes with it.

We find 4 other hosts, Caterpillar $1$, $6$, $27$, and $40$, that have an LMC
mass subhaloes on the first infall. However, in these haloes LMC-sized
subhaloes are still at distances of $201$ kpc, $131$ kpc, $235$ kpc and $225$
kpc, respectively. Many of the faint UFD galaxies brought by these subhaloes
would not be detectable at these distances and we therefore do not consider
them as having an LMC analog. We do find differences in the observed radial
distribution and luminosity function of satellites of such hosts, which we
discuss in Section 3.3.2.

Table 1: Summary of properties of the MW-sized host haloes in the Caterpillar suite that contain a realistic LMC analog (see Section 2.3 for details). For each host we list the epoch closest to $z=0$ at which LMC analog is on its first approach to the pericentre and is at a distance closet to 50 kpc. At this epoch, namely $z_{\rm extract}$, we extract the host $M_{\rm 200c}$ mass, host mass enclosed within 100 kpc ($M(\rm<100kpc)$), host virial radius ($r_{\rm vir}$), host halo concentration ($c_{-2,\rm host}$), redshift at which 50% of the current host halo mass was assembled ($z_{\rm form,50\%}$, the LMC analog $M_{\rm peak}$ mass and its distance from the host halo centre. Note that all the distances below ($r_{\rm vir}$ and $d_{\rm host}$) $d_{\rm host}$ is the proper distance, that is, comoving distance times the scale factor. Simulation Name | $z_{\rm extract}$ | $M_{\rm 200c}$ | $M(\rm<100kpc)$ |  $r_{\rm vir}$(kpc) | c-2 | $z_{\rm form,50\%}$ | LMC $M_{\rm peak}$ ($\,M_{\odot}$) | LMC $d_{\rm host}$(kpc)
---|---|---|---|---|---|---|---|---
Caterpillar 4 | 0.057 | $1.15\times 10^{12}$ | $6.41\times 10^{11}$ | 299.6 | 9.97 | 0.95 | $1.3\times 10^{11}$ | 55.5
Caterpillar 53 | 0.057 | $1.09\times 10^{12}$ | $6.29\times 10^{11}$ | 293.7 | 10.91 | 0.88 | $1.16\times 10^{11}$ | 47.7

###  2.4 Correcting for the difference of host and MW halo masses

In addition to accounting for the presence of the LMC, it is important to make
sure that simulated host haloes have masses consistent with current estimates
of the MW halo mass, because the number of satellites scales roughly linearly
with host halo mass (e.g., Kravtsov et al., 2004).

In our analyses we adopt recent constraints for the MW halo mass enclosed
within 100 kpc (denoted as $M(<100)$) by Vasiliev et al. (2021), who deduced
$M_{\rm MW}(<100)=5.6\pm 0.4\times 10^{11}\,M_{\odot}$. To approximately
correct for the differences in halo mass between MW and simulated haloes, we
multiply the number of satellites in each host by the fraction $M_{\rm
MW,r}(<100)/M_{\rm host}(<100)$, where $M_{\rm MW,r}(<100)$ is a random number
drawn from a Gaussian distribution with the mean of $5.6\times
10^{11}\,M_{\odot}$ and rms of $0.4\times 10^{11}\,M_{\odot}$ and $M_{\rm
host}(<100)$ is the mass enclosed within 100 kpc in the simulated hosts. We
compute the latter by using the virial mass $M_{\rm vir}$, virial radius
$r_{\rm vir}$ and scale radius $r_{\rm s}$ from the halo catalog444Note that
$r_{\rm vir}$ and $r_{\rm s}$ should be in units of proper/physical distance.
Thus, whenever we analyze an epoch with $z_{\rm extract}>0$, we do the
appropriate conversion by multiplying by scale factor at that epoch. and
assuming an NFW profile:

$M(<100\,\text{kpc})=M_{\rm vir}\frac{f(cx)}{f(x)},$ (1)

where $c=r_{\rm vir}/r_{\rm s}$ and $x=100/r_{\rm vir}$ and

$f(x)=\ln(1+x)-\frac{x}{1+x}.$ (2)

There is an additional uncertainty in this normalization factor due to the
differences in the concentration $c$ of the MW halo and our simulated hosts
that we do not account for. In the Appendix C we compare results obtained
using mass correction above with the correction based on $M_{\rm 200c}$ with
MW mass constraints from Deason et al. (2021) or $M(<100)$ from Shen et al.
(2022) and find that these different choices result in a similar median number
of satellites.

###  2.5 The GRUMPY Dwarf Galaxy Formation Model

Figure 3: Model prediction for the relation between stellar metallicity
$\rm\langle[Fe/H]\rangle$ and absolute V-band magnitude $M_{\rm V}$. This is
using $z_{\rm rei}=8.5$ and the galactic outflow parameters assumed in this
analysis. For comparison, we also show observational data of confirmed MW
dwarf galaxies taken from the compilation in Simon (2019) with updated
measurements of Antlia II and Crater II from Ji et al. (2021). Note that this
plot is for only demonstration purposes that our model broadly agrees with
observed dwarf galaxy stellar metallicities. As seen in Figure 7 of Kravtsov &
Manwadkar (2022), the form of this relation is sensitive to galactic outflow
parameters.

The GRUMPY (Galaxy formation with RegUlator Model in PYthon) is a regulator-
type model for dwarf galaxies introduced in Kravtsov & Manwadkar (2022). The
basic premise of this model is that the mass accretion history of each halo in
the simulation, $M_{200\rm c}(t)$, is used to follow evolution of the ISM gas
mass, stellar mass, and mass in metals in these two components using a system
of ordinary differential equations expressing mass conservation and mass
transfer between different components. Using realistic MAHs extracted from
high-resolution simulations allows to account for the intrinsic scatter of
observed dwarf galaxy properties arising from the scatter in their MAHs.

The model includes prescriptions for UV heating after reionization and
associated gas accretion suppression onto small mass haloes, galactic
outflows, model for gaseous disk and its size, molecular hydrogen mass, star
formation, etc. The model includes a prescription for the evolution of half-
mass radius of the stellar distribution. This is important for modelling the
surface brightnesses of dwarf galaxies and thus their detectability in a
particular observational survey. Integrating the galaxy model with the MAHs
extracted from simulations, we obtain the time evolution of the baryonic
components, including star formation histories and metallicity evolution.

The galaxy model parameters used in this paper are the fiducial parameters in
Kravtsov & Manwadkar (2022), except for the galactic outflow parameters
$\eta_{\rm norm}$ and $\eta_{\rm p}$ which we assume to be $\eta_{\rm
norm}=1.8$ and $\eta_{\rm p}=-0.45$. We modify these parameters555The fiducial
parameters assumed in Kravtsov & Manwadkar (2022) were $\eta_{\rm norm}=3.6$
and $\eta_{\rm p}=-0.35$ to better match the bright end of the MW satellite
luminosity function. Such modifications are reasonable as they are within the
current scatter/uncertainty of the mass loading factor (e.g., Mitchell et al.,
2020; Pandya et al., 2021).

Our galaxy model produces realistic mean stellar metallicities and
metallicity-$M_{\star}$ relations in good agreement with observations of MW
satellites, as shown in Figure 3. The normalization and slope of this relation
is sensitive to galactic outflow parameters (see Kravtsov & Manwadkar, 2022).
We thus choose the outflow model parameters to reproduce the observed relation
and a number of other observations of dwarf galaxies as demonstrated in our
previous paper. Note that our model galaxies simultaneously match the observed
stellar mass–metallicity relation and satellite luminosity functions (as we
show in Section 3.1) without needing preventative feedback due to preheating
as previous studies (e.g., Lu et al., 2017) have suggested (see Section 4 in
Kravtsov & Manwadkar, 2022, for discussion of why this is the case).

###  2.6 Forward modelling observational satellite samples

We construct mock dwarf galaxy catalog using the galaxy formation model run
with the mass accretion histories of all haloes that satisfy the criteria
described in Section 2.1. The star formation histories and stellar metallicity
evolution predicted by the model are used along with the Flexible Stellar
Population Synthesis (FSPS; Conroy et al., 2009; Conroy & Gunn,
2010)666https://github.com/cconroy20/fsps model to compute magnitudes of our
galaxies in different bands assuming a Chabrier (2003) initial mass function
and standard FSPS parameters.

Specifically, in this project we compute the absolute V-band magnitude,
$M_{\rm V}$, and azimuthally averaged effective half-light radius, $r_{1/2}$.
In practice, we convert from 3D stellar half-mass radius $r_{\rm 50,3D}$
predicted by the model to project half mass radius, assuming that stellar
distribution is approximately spherical: $r_{1/2}\approx r_{\rm 50,3D}/1.34$
(Wolf et al. 2010). The 1.34 conversion factor is for a spherical galaxy, but
is expected to depend on ellipticity (see, e.g., Section 4.2.2 and eq. 10 in
Behroozi et al., 2021): for a galaxy with $c/a=0.5$ the factor is
$\approx$1.54 instead of 1.33 – the $\approx$15% difference that we ignore in
our study.

We then forward model mock dwarf galaxy populations to compare with observed
satellite population in specific surveys. To take into account uncertainty in
the observer location, we choose 200 random observer locations at a physical
distance of 8 kpc from the host halo centre. As the Caterpillar simulations do
not model the central baryonic disk potential, we do not constrain the random
observer locations to any plane. For each observer location, we rotate the
entire reference frame so that the LMC analog matches with the true sky
coordinates of the LMC (RA = 80.89∘, Dec = -69.75∘). In the case of hosts with
no LMC analogs, we choose random pointings of survey footprint in the sky.
Thus, for each random observer location and each dwarf galaxy in our mock
catalog, we have heliocentric distance ($d_{\rm helio}$), absolute V-band
magnitude ($M_{\rm V}$), effective half-light radius ($r_{1/2}$) and sky
coordinates. Similar to Nadler et al. (2020), we only consider galaxies with
$r_{\rm 1/2}>10$ pc, which is the threshold used to define observed dwarf
galaxy sample to avoid confusion with star clusters (e.g., Drlica-Wagner et
al., 2020). We discuss implications of adopting this threshold and possibility
that a significant fraction of observed stellar systems with $r_{\rm 1/2}<10$
pc may be dwarf galaxies in Section 4.1.

We use the resulting galaxy catalogs and prescriptions for the detection
probability for the DES, PS1, SDSS, DELVE-WIDE, HSC-SSP and Vera C. Rubin
Observatory/LSST survey to compare with observations and make future
predictions. The prescriptions we use in our analysis are detailed in the
remainder of this section.

####  2.6.1 The SDSS Selection Function

To model the SDSS satellite population, we use the selection function
described in Koposov et al. (2009), which is expressed in terms of threshold
magnitude of detection as a function of heliocentric distance:

$M_{\rm V}<\frac{1.1-\log_{10}(d_{\rm helio}/\text{kpc})}{0.228}$ (3)

We assign zero detection probability for satellites with a magnitude fainter
than above threshold or with the $V$-band surface brightness, $\mu_{\rm
V}>30\rm\,mag\,arcsec^{-2}$ where

$\mu_{\rm V}=M_{\rm V}+21.5721+2.5\log_{10}(2\pi r^{2}_{1/2})$ (4)

for $r_{1/2}$ in pc and $\mu_{V}$ in mag/arcsec2.

The SDSS selection function is thus a step function with probability of $1.0$
if conditions are satisfied and zero otherwise. Furthermore, we apply the
selection thresholds only to galaxies within the SDSS DR9 imaging footprint of
the total area of 14555 deg2 (Ahn et al. 2012)

We use the approach of Jethwa et al. (2018) to jointly model the classical
satellites detected prior to SDSS and SDSS-detected systems. Namely, we assume
satellites with $M_{\rm V}\leq-11$ are detected over the entire sky, while
fainter satellites are modelled using the SDSS selection criterion described
above and are restricted to the SDSS footprint.

####  2.6.2 The DES & PS1 Selection Function

We use the DES and PS1 survey selection function described in Drlica-Wagner et
al. (2020) and implemented in a publicly available
code.777https://github.com/des-science/mw-sats The selection functions are
based on a machine learning method trained to compute satellite detection
probability as a function of absolute magnitude $M_{\rm V}$, heliocentric
distance $d_{\rm helio}$, azimuthally averaged half-light radius $r_{1/2}$,
and sky position (RA, Dec). The information about the local stellar density
and survey footprint is contained within the selection function code.

####  2.6.3 DELVE-WIDE Selection Function

The DELVE-WIDE survey888https://delve-survey.github.io/ (e.g., Drlica-Wagner
et al., 2021) is a survey that will cover the entire high-Galactic-latitude
range ($|b|>10^{\circ}$) in the southern hemisphere using the Dark Energy
Camera. DELVE-WIDE will reach a photometric depth that is similar to that of
the DES. We thus use the DES selection function discussed above to make
predictions for DELVE-WIDE dwarf galaxy searches.

Note that in the DES selection function, the detection probability calculation
is separate from footprint mask value calculation. So we simply modify the
footprint mask value calculation by using the DELVE-WIDE footprint instead of
the DES footprint. The total area of the DELVE-WIDE footprint is
$\sim$17000$\,\rm deg^{2}$, and includes both the $\sim$2200$\,\rm deg^{2}$
area around the Magellanic Clouds observed as part of the DELVE-MC survey and
the DES footprint (Drlica-Wagner et al., 2021).

Figure 4: Predicted satellite luminosity function for hosts with LMC analog
within $d_{\rm helio}<300$ kpc in the DES (left), PS1 (centre) and SDSS +
Classical (right) samples as a function of reionization redshift $z_{\rm
rei}=7.5,8.5,9.5,10.5$. The dark (light) orange shaded regions for $z_{\rm
rei}=8.5$ corresponds to the 68% (95%) of the model predictions due to host-
to-host scatter, bootstrap variations in the disruption model, etc. The
colored lines with different line styles denote the median luminosity
functions for different $z_{\rm rei}$ models. The solid black lines shows the
observed luminosity function for the respective survey. Note that the ‘SDSS +
Classical’ satellite luminosity function is combination of satellites observed
in SDSS footprint and all classical satellites ($M_{\rm V}\leq-11$) over
entire sky as discussed in Section 2.6.

####  2.6.4 HSC-SSP and LSST Y1 Selection Function

The selection functions for the ongoing HSC-
SSP999https://hsc.mtk.nao.ac.jp/ssp/ (e.g., Aihara et al., 2018a, b) and
upcoming LSST surveys (e.g., Ivezić et al., 2019) are not yet rigorously
defined and the detection efficiency in the presence of high-density of faint
background galaxies is not yet fully understood. We therefore adopt
approximate selection function using approximations for the detection
efficiency through simple surface brightness or magnitude cuts, as was done in
previous recent studies (e.g., Homma et al., 2018, 2019; Nadler et al., 2019;
Nadler et al., 2020).

Specifically, following Nadler et al. (2020), we assume that HSC-SSP detects
all satellites with $M_{\rm V}<0$, $r_{\rm 1/2}>10$ pc, and $\mu_{\rm V}<32$
mag/arcsec2 within a heliocentric distance of 300 kpc. LSST Year 1 (Y1) is
expected to have a detection efficiency similar to that of the HSC-SSP survey
due to similar limiting magnitude. So we assume the same selection function as
HSC-SSP for making predictions of dwarf galaxy observations in the LSST Y1
survey. Challenges of dwarf galaxy observations in such deep surveys are
discussed further in Section 4.3.

To forward model HSC-SSP observations, we use the HSC-SSP footprint definition
described in Table 5 and Figure 3 in Aihara et al. (2018a). At the moment,
only 676 deg2 of a total area of $\sim$1400$\,\rm deg^{2}$ in the HSC-SSP
footprint has been searched for satellite galaxies. To forward model the HSC-
SSP observations in this current (partially observed) area, we consider a
galaxy to be observed if – i) it passes the HSC-SSP detection threshold, ii)
it falls in the HSC-SSP footprint and iii) $\rm U(0,1)<676/1400=0.48$ where
$\rm U(0,1)$ denotes a random number drawn from the uniform distribution in
the interval $[0,1]$. In other words, we only consider a $676/1400=0.48$
fraction of galaxies that pass the HSC-SSP selection function. As the HSC-SSP
selection function described above is only an approximation, we do not include
HSC-SSP satellite observations in our main analysis. However, we discuss the
forward modelled predictions for HSC-SSP in Section 4.3.

In the Big Sky approach for the LSST101010https://www.lsst.org/ survey
described in Olsen et al. (2018), the survey footprint will cover the area the
entire equatorial declination range of $\delta<30^{\circ}$. However, as
central portions of the Galactic plane within this footprint will be removed,
we approximate this by considering only the high-Galactic-latitude range
($|b|>15^{\circ}$) within the LSST footprint. This corresponds to a total
footprint area of $\sim$23500$\,\rm deg^{2}$.

##  3 Model Results and Comparisons to Observations

Here we present results for the population of satellite galaxies resulting
from forward-modelling described in the previous section. When possible, we
compare model predictions with observed satellite population in the DES, PS1,
and SDSS surveys. The list of observed dwarf galaxies used in these
comparisons is taken from the compilation in Drlica-Wagner et al. (2020). Note
that our comparison sample of the MW satellites excludes the recently
discovered massive satellites Antlia II (Torrealba et al. 2019b) as it is not
detected in any of these surveys. Also, even though satellites Kim 2 and
Crater I/Laevens 1 have $r_{\rm 1/2}>10$ pc and have been detected with high
significance in the DES and PS1 surveys respectively, we do not include them
in our satellite samples because these systems are considered to be likely
star clusters (Kirby et al. 2015; Weisz et al. 2016). As a result, these two
objects have been listed among probable star clusters in Table 5. The
exclusion/inclusion of a single satellite in the DES and PS1 satellite samples
respectively does not change our conclusions.

We begin by comparing luminosity function and radial distribution of
satellites with respect to halo centre predicted by the forward modelling to
observations in Sections 3.1 and 3.3, respectively. Unless stated otherwise,
in all of the comparisons with observations we use the two hosts with LMC
analogs, as discussed in Section 2.3. We explore sensitivity of these
statistics to reionization redshift of the Milky Way volume, $z_{\rm rei}$
(Section 3.1.1), and to the presence of the LMC-like object (Section 3.1.2),
and derive the approximate range of $z_{\rm rei}$ consistent with observed
abundance of ultra-faint galaxies (Section 3.1.3). We use this range of
$z_{\rm rei}$ to derive the implied relations between galaxy stellar mass and
peak halo mass in our model (Section 3.2.1) and halo occupation distribution
for galaxies with $M_{\rm V}\leq 0$ (Section 3.2.2). Furthermore, we discuss
the temporal variation of satellite radial distributions, under the influence
of LMC infall for instance, in Section 3.3.3.

Given the importance of size (or, equivalently, surface brightnesses) for a
galaxy of a given luminosity for its detectability, we also consider
distribution of $r_{\rm 1/2}$ sizes of dwarf satellites in Section 3.4.

###  3.1 Luminosity Functions

In this section we present comparisons of predicted luminosity functions for
different observational surveys and their dependence on the reionization
epoch.

####  3.1.1 Luminosity functions and redshift of reionization

Figure 5: Predicted satellite luminosity function for $z_{\rm rei}=8.5$ within
$d_{\rm helio}<300$ kpc in the DES (left), PS1 (centre) and SDSS $+$ Classical
(right) samples for different LMC scenarios. The orange and red solid lines
shows the median luminosity function for hosts with ‘LMC analog’ hosts and
‘without LMC’ respectively. The orange contours carry the same meaning as in
Figure 4 and the red dashed-lines bracket the 68% confidence intervals of the
‘without LMC’ scenario. The teal colored dotted and dot-dashed lines
correspond to the two alternative LMC scenarios considered in Section 3.1.2.
Our forward modelling results prefer a LMC that is on its first infall and
near its pericentre (‘LMC analog’ scenario).

Figure 4 shows the predicted satellite luminosity function for the DES, PS1
and SDSS + Classical satellites in models with reionization redshifts of
$z_{\rm rei}=7.5,8.5,9.5,10.5$ and compares them with the observed LF in each
survey. The colored lines with different line styles show the median
luminosity functions for the respective reionization redshift $z_{\rm rei}$,
while the dark (light) shaded regions denote the 68% (95%) confidence
intervals.

As was discussed above, reionization is known to directly impact the amplitude
of the faint end of the luminosity function. This is expected because the UV
background after reionization heats the intergalactic medium and prevents gas
accretion onto haloes below a certain time dependent characteristic mass scale
(e.g., Gnedin, 2000; Okamoto et al., 2008). This suppresses buildup of stellar
mass in low mass haloes and consequently the number of faint satellites. This
can be seen in Figure 4, which shows that the amplitude of satellite LF at
$M_{\rm V}\gtrsim-6$ is increasingly suppressed for higher $z_{\rm rei}$. The
trend with $z_{\rm rei}$ reflects higher degree of suppression when haloes
fall below the critical mass for accretion earlier. haloes of mass larger than
critical mass ($M_{\rm peak}\gtrsim 10^{9}\,M_{\odot}$) are nearly insensitive
to UV heating, and satellite LF is thus insensitive to $z_{\rm rei}$ for
$M_{\rm V}\lesssim-6$. The bright end of the luminosity function, however, is
sensitive to galactic outflow modelling (Kravtsov & Manwadkar, 2022).

It is interesting to point out that the luminosity of $M_{\rm V}\gtrsim-8$
used to delineate UFDs and classical dwarf galaxies (e.g., Simon, 2019) is
very close to the luminosity below which galaxies are significantly affected
by UV heating and reionization redshift. This reiterates the importance of
UFDs as probes of reionization.

The faint end LF amplitude can thus be used to probe the redshift of
reionization of the Lagrangian region of the Milky Way. Results shown in
Figure 4 indicate that $z_{\rm rei}$ in the considered range are in general
agreement with observed LF. The $z_{\rm rei}\sim 8$ model matches the observed
LF most closely and we will consider this model as our fiducial model in this
study. We return to how $z_{\rm rei}$ for the MW environment could be
constrained in more detail in Section 3.1.3.

####  3.1.2 Effect of the Large Magellanic Cloud

A number of UFDs have been discovered around the Magellanic Clouds (e.g.,
Koposov et al., 2015; Koposov et al., 2018; Drlica-Wagner et al., 2015b, 2016;
Torrealba et al., 2018; Cerny et al., 2021a) consistent with the hypothesis
that the Magellanic Clouds arrived with their own population of dwarf galaxy
satellites (D’Onghia & Lake, 2008). Indeed, Nadler et al. (2020) recently
showed that the presence of the LMC has a significant impact on the MW
satellite population, especially in the vicinity of the LMC. It is therefore
important to account for the presence of LMC when forward modelling satellite
populations in surveys with footprint near this object.

To study the LMC’s impact on satellite population, we consider the satellite
LF and radial distribution (in Section 3.3.2) in four different LMC scenarios:
i) our fiducial case of hosts with a LMC analog as defined in Section 2.3 (LMC
analog), ii) hosts without any LMC mass object accreted (wo/LMC), iii) hosts
that accreted a LMC early on so it has undergone multiple pericentric passages
by $z=0$ (Early LMC) and iv) hosts that have an LMC mass object on its first
infall, but at distances $>100$ kpc (‘Far LMC’). For the “Early LMC” and “Far
LMC” scenarios we rotate the reference frame such that the LMC halo matches
with true sky coordinates of LMC. For the host halos without LMC, we assume a
random orientation of the host halo relative to a survey footprint.

Figure 5 shows the predicted satellite luminosity functions at $z_{\rm
rei}=8.5$ in the DES, PS1 and SDSS$+$Classical samples for the four different
LMC scenarios. There is no statistically significant impact of the LMC on the
PS1 and SDSS+Classical satellite luminosity function, except the trivial
effect that brightest satellites are lacking in the SDSS+Classical model in
hosts without LMC mass satellites. However, the figure shows that presence of
a LMC analog boosts the amplitude of the faint-end LF predicted for the DES
survey. This is not surprising as LMC lies close to the DES footprint (see,
e.g., Figure 1 in Drlica-Wagner et al. 2020). Indeed, if we reflect satellite
positions in declination in the LMC hosts, effectively moving LMC analog away
from the DES footprint, the amplitudes of the DES satellite LFs for hosts with
and without LMC (wo/LMC) become similar. This explicitly confirms that LMC has
significant impact on the DES satellite population due to its proximity to the
DES footprint.

Furthermore, we find that the ‘Far LMC’ and ‘Early LMC’ have little impact on
satellite LF and produce results similar to the host halos without LMC-sized
subhalo. This is because in the ‘Far LMC’ scenario the LMC associated
satellite population is still far from the host and hence we observe fewer
satellites due to observational survey selection limits. In the ’Early LMC’
scenario, the LMC has undergone multiple pericentre passages and so a
significant fraction of the associated satellite population has been disrupted
or dispersed of into the MW halo (e.g., Deason et al., 2015), reducing the
over-density of satellites around the LMC.

####  3.1.3 Implications for Reionization of the Local Environment

Figure 6: The number of satellites with $M_{\rm V}>-6$ and $d_{\rm helio}<300$
kpc predicted to be observed in total in the DES and PS1 surveys as a function
of reionization redshift $z_{\rm rei}$. Note that $z_{\rm rei}$ here is
reionization of the Local Volume, which can be different from the mean
reionization redshift of the Universe. The horizontal dot-dash line ($N_{\rm
sat}=25$) shows the number of satellites with $M_{\rm V}>-6$ observed in the
DES and PS1 surveys above a certain detection threshold as considered in
selection function. Even though the $z_{\rm rei}$ consistent with observed
satellite abundance is quite wide, it is in good agreement with the
reionization redshift range of $z_{\rm rei}=8.1^{+0.4}_{-0.4}$ (95% confidence
intervals) derived for the MW mass haloes by Zhu et al. (2019) using
simulations of reionization.

Given that the redshift of reionization directly impacts the number of UFD
satellites, the observed number of MW UFD satellites can be used to constrain
the reionization redshift of the Lagrangian volume of the Milky Way. To this
end, we use the total number of UFDs observed in the DES and PS1 surveys and
identify the range of $z_{\rm rei}$ for which predictions of the forward model
are consistent with the observed number.

Figure 6 shows the number of luminous satellites at $M_{\rm V}>-6$ we would
expect to observe in the DES and PS1 as a function of $z_{\rm rei}$. The
scatter about the median predicted $N_{\rm sat}-z_{\rm rei}$ relation is
dominated by the host-to-host scatter in the satellite (subhalo) counts, even
after accounting for uncertainty in the MW halo mass as described in Section
2.4.

Model predictions are consistent with the observed UFD number for reionization
redshifts of $z_{\rm rei}\sim 7-10$. The actual uncertainty in $z_{\rm rei}$
will be slightly larger because we do not account for uncertainties in
modelling of galactic outflows, which will impact the number of luminous
satellites. Even though our inferred $z_{\rm rei}$ range is quite wide, it is
interesting to see that our inferred $z_{\rm rei}$ range agrees with the
reionization redshift range of $z_{\rm rei}=8.1^{+0.4}_{-0.4}$ (95% confidence
intervals) derived for the MW mass haloes by Zhu et al. (2019) using the CROC
simulations of reionization. Our constraint is also in qualitative agreement
with results of the CoDa II simulations which show that outer regions of the
MW and M31 Lagrangian regions are reionized at $z\gtrsim 7$ (Ocvirk et al.,
2020, see their Fig. 11) and with the general expectation that the Lagrangian
regions around progenitors of the present-day collapsed structures reionize at
higher redshifts than the low-density regions and the universe as a whole
(e.g., Trac et al., 2022).

The sensitivity of the UFD abundance to reionization shows that this is one
avenue to probe reionization in the local volume when census of the UFDs will
improve with larger and deeper surveys like the DELVE-WIDE with the Dark
Energy Camera and LSST survey at the Vera C. Rubin Observatory.

###  3.2 Implications for the Galaxy-Halo Connection

####  3.2.1 Relations between peak halo mass, stellar mass and $V$-band
luminosity

Figure 7: Model prediction for the relation between peak halo mass $M_{\rm
peak}$ and stellar mass $M_{\star}$ for the Milky Way satellite galaxies. The
solid colored line denotes the median relation while the dark (light) shaded
region correspond to the 68% (95%) scatter around the median relation. The
purple dotted line shows the best fit relation derived from fitting the MW
satellite population by Nadler et al. (2020). The black dashed line shows the
approximation to the median relation given by Eqs. 7 in Appendix B. Note that
the significant down scatter in this relation at masses $M_{\rm peak}\lesssim
10^{9}\,M_{\odot}$ is due to some fraction of haloes not forming any
significant amount of stars.

Figure 7 shows relation between peak halo mass $M_{\rm peak}$ and stellar mass
$M_{\star}$. $M_{\rm peak}$ is defined as the maximum mass that a given halo
achieved during its evolution at some redshift $z_{\rm peak}$ and it is
measured within the radius enclosing density contrast of $200$ relative to the
critical density at $z_{\rm peak}$. The scatter around the median $M_{\rm
peak}-M_{\star}$ relation shown in the figure reflect both halo-to-halo
scatter and the range of $z_{\rm rei}$ consistent with observed abundance of
dwarf satellites (Figure 6). The figure shows that the $M_{\rm
peak}-M_{\star}$ relation arising in the model is not a single power law. This
is because its slope at $M_{\rm peak}\gtrsim 10^{10}\,M_{\odot}$ ($M_{\rm
V}\lesssim-11$) mainly reflects the scaling of the wind mass loading factor
with stellar mass, while at smaller masses the slope is determined by
suppression of gas accretion due to UV heating after reionization (see
Kravtsov & Manwadkar, 2022, for detailed discussion). The flattening of the
relation at the smallest $M_{\rm peak}$ values reflects the increasing
inefficiency of supernovae driven winds in blowing out gas from these galaxies
(see also Chen et al., 2022).

The dashed line in Fig. 7 shows an accurate analytical approximation to the
median $M_{\rm peak}-M_{\star}$ of the model, which is described in the
Appendix B. The figure also compares the relation in our model with the
relation derived by Nadler et al. (2020) using observed abundance of
satellites in the DES and PS1 surveys. In their analysis, a power law form of
the $M_{\rm peak}-M_{\star}$ was assumed with normalization that was fixed at
$M_{\rm peak}\gtrsim 10^{10}\,M_{\odot}$ to the normalization derived using
abundance matching with observed stellar mass function. Overall, the relation
in our model is close to the power law relation derived by Nadler et al.
(2020), but with some differences due to deviations from the power law
described above and in amplitude at the largest satellite masses (see Appendix
B for additional comparison). As we show in Figure 21, our relation is also
consistent with recent, state-of-the-art hydrodynamic simulations of dwarf
galaxies (e.g., Applebaum et al., 2021; Orkney et al., 2021; Santistevan et
al., 2022; Gutcke et al., 2022; Hopkins et al., 2022) at the halo mass range
probed in these simulations.

As our comparisons with observed satellite LF show, the non-power law $M_{\rm
peak}-M_{\star}$ relation in our model is consistent with observed
distribution of satellite luminosities. We note that the current data is not
sufficient to reliably exclude a single power-law $M_{\rm peak}-M_{\star}$
relation (Nadler et al., 2020), although as we show in Section 4.3 it is
easier to explain the abundance of satellite in the HSC-SSP survey with the
flattening of the $M_{\rm peak}-M_{\star}$ relation in our model.

####  3.2.2 Halo Occupation Distribution of $M_{V}<0$ Galaxies

Figure 8: Fraction of subhaloes that host galaxies with absolute V-band
magnitudes of $M_{\rm V}<0$ as a function of peak halo mass $M_{\rm peak}$.
The dark (light) shaded regions show the 68% (95%) confidence intervals on our
predictions for the halo occupation fraction for MW satellites. The blue
upward arrow shows the 95% confidence for the peak halo mass that hosts 50% of
detectable satellites from Nadler et al. (2020). For comparison sake, we also
show the median halo occupation fraction (black dotted line) if reionization
was not modeled. Figure 9: Comparisons of the model and observed distribution
of the heliocentric distance of MW satellites in the DES (left), PS1 (centre)
and SDSS + Classical (right) as a function of different reionization redshift
$z_{\rm rei}=7.5,8.5,9.5,10.5$. The distributions are normalized to the total
number of satellites within the heliocentric distance of 300 kpc. The solid
black line shows the observed satellite radial distribution in each of the
surveys, while colored lines show model predictions for different $z_{\rm
rei}$. The dark (light) orange shaded regions for $z_{\rm rei}=8.5$ correspond
to the 68% (95%) of the model predictions.

Suppression of gas accretion and star formation due to UV heating affects
luminosities of galaxies hosted in dark matter subhaloes and thus their
occupation distribution. In Kravtsov & Manwadkar (2022) we showed that the
fraction of haloes at a given $M_{\rm peak}$ that host a luminous satellite
($M_{\rm V}<0$) is sensitive to the redshift of reionization for $M_{\rm
peak}\lesssim 10^{9}\,\,M_{\odot}$. Thus, the range of $z_{\rm rei}$ favored
by the observed number of the DES and PS1 dwarf satellites implies a specific
constraint on the halo occupation fraction that we show in Figure 8. The
shaded regions around the median shown by the shaded region correspond to the
68% (95%) confidence intervals in the halo occupation fraction reflecting the
corresponding regions of the $z_{\rm rei}$ range. In addition, we show 95%
confidence upper bound111111Note that lower $\mathcal{M}_{50}$ corresponds to
higher halo occupation and vice-a-versa. on the halo mass corresponding to
$f(M_{V}<0)=0.5$ of $\mathcal{M}_{50}<8.5\times 10^{7}\,M_{\odot}$ derived by
Nadler et al. (2020). The model consistent with the observed number of DES and
PS1 satellites corresponds to $\mathcal{M}_{50}=3.6^{+2.6}_{-1.4}\times
10^{7}\,M_{\odot}$ (68% confidence intervals).

This constraint is qualitatively consistent with results of previous studies
(e.g., Macciò et al., 2010; Jethwa et al., 2018; Nadler et al., 2020), which
concluded that a significant fraction of UFDs must be hosted by subhaloes with
$M_{\rm peak}\lesssim 10^{8}\,M_{\odot}$. In particular, as can be seen in the
figure, the $f(M_{V}<0)$ of our model is consistent with the
$\mathcal{M}_{50}$ of Nadler et al. (2020). However, the range of halo masses
hosting $M_{V}<0$ galaxies allowed by the data that was quantified in the
latter study is quite broad (reflecting limitations of the current satellite
sample). The current model predicts the HOD form and $\mathcal{M}_{50}$ range
that are quite narrow. These reflects the specific assumptions made our model.

The halo occupation fraction implied by observations in the context of our
model, indicates that observed UFD satellites of the MW form in haloes as
small as $\sim 10^{7}\,\,M_{\odot}$, which is below the putative atomic
cooling limit (e.g., Haiman et al., 1997; Barkana & Loeb, 1999). We discuss
this further in Section 4.2.

###  3.3 Radial Distributions

As discussed in Section 1, radial distribution of satellites with respect to
host halo centre is expected to be sensitive to subhalo disruption modelling,
details of merger history, observational biases, and possibly even
reionization. In this section, we present a comparisons of predicted radial
distributions for different observational surveys and their dependence on
reionization epoch and presence of the LMC.

####  3.3.1 Radial distribution of satellites and redshift of reionization

Figure 9 shows the predicted satellite radial distributions for the DES, PS1
and SDSS+Classical satellites for different values of $z_{\rm rei}$ and radial
distributions of observed satellites in these surveys, both normalized to the
number of detected satellites within heliocentric distance of 300 kpc.
Overall, we see that the radial distributions are in good agreement with
observations and this agreement is achieved without modelling orphan galaxies
to account for subhaloes missing due to numerical effects.

For a more quantitative assessment, Table 2 presents the $p$-values of the
Kolmogorov-Smirnov (KS) test for the radial distributions predicted by the
model with different $z_{\rm rei}$ compared to observed radial distributions
in different surveys. The $p$-values for models with $z_{\rm rei}\sim 7-10$
are generally high in all survey comparisons for models that use host haloes
with LMC analogs, which indicates that the radial distribution of model
galaxies is in good agreement with observed distribution of the MW satellites.

As discussed in Section 3.1.1, earlier reionization (higher $z_{\rm rei}$)
results in an increased suppression in the number of luminous satellites. This
can be seen in Figure 9 where the amplitude of the radial distribution
decreases with increasing $z_{\rm rei}$. However, the overall shape or
concentration of the radial distribution does not change with $z_{\rm rei}$.
This is not surprising because the model assumes the same $z_{\rm rei}$ for
all subhaloes regardless of their spatial location or accretion redshift.
Although reionization, in general, is expected to proceed inhomogeneously
(e.g., Ocvirk et al., 2014) and satellites closer to the MW could be reionized
earlier, mixing of satellite orbits within MW will weaken any such trend
considerably. Moreover, the difference in $z_{\rm rei}$ across Lagrangian
region of the Milky Way is expected to be only $\Delta z_{\rm rei}\lesssim 1$
(e.g., see Figure 11 in Ocvirk et al., 2020) which would make any such trend
with $z_{\rm rei}$ small.

Table 2: Kolmogorov-Smirnov (KS) test results for model predicted radial distribution compared to observed radial distribution in the DES, PS1 and SDSS + Classical satellite samples. Model |  $p$-value
---|---
DES, LMC analog, $z_{\rm rei}=7.5$ | 0.65
DES, LMC analog, $z_{\rm rei}=8.5$ | 0.73
DES, LMC analog, $z_{\rm rei}=9.5$ | 0.90
DES, LMC analog, $z_{\rm rei}=10.5$ | 0.85
DES, without LMC, $z_{\rm rei}=8.5$ | 0.40
DES, far LMC, $z_{\rm rei}=8.5$ | 0.01
DES, early LMC, $z_{\rm rei}=8.5$ | 0.54
PS1, LMC analog, $z_{\rm rei}=7.5$ | 0.46
PS1, LMC analog, $z_{\rm rei}=8.5$ | 0.34
PS1, LMC analog, $z_{\rm rei}=9.5$ | 0.21
PS1, LMC analog, $z_{\rm rei}=10.5$ | 0.13
PS1, without LMC, $z_{\rm rei}=8.5$ | 0.03
PS1, far LMC, $z_{\rm rei}=8.5$ | 0.02
PS1, early LMC, $z_{\rm rei}=8.5$ | 0.08
SDSS + Classical, LMC analog, $z_{\rm rei}=7.5$ | 0.95
SDSS + Classical, LMC analog, $z_{\rm rei}=8.5$ | 0.99
SDSS + Classical, LMC analog, $z_{\rm rei}=9.5$ | 0.96
SDSS + Classical, LMC analog, $z_{\rm rei}=10.5$ | 0.93
SDSS + Classical, without LMC, $z_{\rm rei}=8.5$ | 0.49
SDSS + Classical, far LMC, $z_{\rm rei}=8.5$ | 0.07
SDSS + Classical, early LMC, $z_{\rm rei}=8.5$ | 0.92

####  3.3.2 Effect of the Large Magellanic Cloud

Figure 10: Predictions for the satellite radial distribution in DES (left),
PS1 (centre) and SDSS + Classical (right) satellite samples for $z_{\rm
rei}=8.5$ for different LMC scenarios. Similar to Figure 5, the orange and red
lines show the show the radial distributions for hosts with ‘LMC analog’ and
‘without LMC’ respectively. The teal colored dotted and dot-dashed lines
correspond to ‘Far LMC’ and ‘Early LMC’ scenario as discussed in Section
3.1.2. Similar to luminosity functions, our forward modelling results prefer a
LMC that is on its first infall and near its pericentre (‘LMC analog’
scenario).  Figure 11: The temporal variation in the concentration of the
predicted satellite radial distributions ($R_{\rm 50\%,survey}$) for $z_{\rm
rei}=8.5$ in the DES (left), PS1 (middle) and SDSS+Classical (right) samples.
Note that $\tau_{0}$ is the age of the Universe in Gyr and $R_{\rm
50\%,survey}$ is the heliocentric radius at which half of the detected
satellites within $d_{\rm helio}<300$ kpc are enclosed. The two colored lines
(red and orange) show the variation for the two hosts with LMC analogs. The
1$\sigma$ uncertainties on $R_{\rm 50\%,survey}$ are only shown for
Caterpillar 53, however, the uncertainties for Caterpillar 4 are of similar
scale. The solid grey line and shaded grey region show the median and 68%
confidence interval on the variation of concentration in hosts with no LMC
analogs. The vertical colored dot-dashed and dotted lines show the epoch when
the LMC analog is crossing $300$ kpc and at pericentre respectively. One
observes a significant decrease in $R_{\rm 50\%,survey}$ of DES satellites as
the LMC approaches pericentre.

As with the satellite LF, LMC has a significant impact on the predicted
satellite radial distribution for the DES, and a much smaller impact on the
predicted distributions for the PS1 and SDSS + Classical samples. In addition
to the amplitude of the radial distribution, the presence of the LMC affects
concentration of the radial distribution of the DES satellites, as can be seen
in the comparison of the normalized radial distributions in Figure 10.

The effect of the LMC on the radial distribution does depend on the location
and accretion redshift of the LMC analog. In the ‘Early LMC’ analogs (see
Section 2.3 for details on different LMC scenarios) already experienced
multiple pericentre passages. The satellites brought by the LMC analog are
likely already dispersed into the MW halo after $\gtrsim 2$ pericentre
passages (e.g., Deason et al., 2015). In the ‘far LMC’ analogs, the majority
of the LMC satellites are located at large distances and thus would not be
detectable (e.g., Zhang et al., 2019; Bose et al., 2020), which minimizes
their impact. Indeed the ‘Far LMC’ analogs have the least concentrated
satellite radial distribution among different LMC analogs. Not surprisingly,
the ‘LMC analog’ scenario in which the current distance to the LMC is matched
agrees best the observed DES satellite radial distribution.

The agreement of the radial distributions in hosts with different LMC analogs
can be quantitatively assessed using the KS $p$-values in Table 2. In the DES
survey, where the effect of LMC is expected to be most noticeable, the
observational data clearly prefers the ‘LMC analog’ scenario over the other
LMC scenarios we consider here. In other words, our forward modelling
framework along with the DES satellite observations provide strong evidence
for a LMC that is on its first infall into the MW and is very close to its
pericentre, in good agreement with conclusions of Nadler et al. (2020).
Interestingly, even though a difference in amplitudes of satellite luminosity
functions for PS1 and SDSS + Classical was not seen, the radial distributions
for these satellites are sligthly more concentrated in the ’LMC analog’
scenario, as can been seen from Figure 10 and $p$-values in Table 2.

In conclusion, predictions of our forward modelling of the MW satellite
population are in excellent agreement with the observed satellite radial
distributions in different surveys. Additionally, we show that the radial
distribution of the observed MW satellites clearly prefers the scenario in
which LMC is on its first infall, is close to the pericentre, and has
contributed a significant number of UFD satellites to the MW satellite
population.

.

####  3.3.3 Temporal variation of satellite radial distributions

Due to the relatively small number of observed satellites in a given survey
and their respective motions within the MW halo, the observed satellite radial
distributions could experience significant fluctuations over relatively short
time intervals of $\sim$1 Gyr. When comparing to observations we used model
predictions for a single time epoch $z=0$. However, in the case of the LMC-
analog hosts, we present results at the epoch ($z=0.057$) when LMC-analog is
on its first infall and is at a distance closest to 50 kpc.

We examine temporal variations of the radial distribution of satellites using
$R_{\rm 50\%,survey}$ – the heliocentric distance enclosing half of the
detected satellites within 300 kpc in a survey. To compute the $R_{\rm
50\%,survey}$ at each epoch, we follow the same forward modelling procedure as
described in Section 2.6 using the $M_{\rm V}$ and $r_{\rm 1/2}$ of the
subhaloes at $z=0$. We do this because we wish to focus on the satellite
motion in MW halo and so decouple the possible systematic of evolution in
satellite observable properties. Note that for Caterpillar 4 and Caterpillar
53 (the host haloes with LMC analogs), we shift our reference frame/survey
pointings such that the it matches with the true sky co-ordinates of the LMC.

Figure 12: (Left) The 2D distribution of all our model dwarf galaxies in
azimuthally averaged physical half light radius ($r_{1/2}$) and absolute
V-band magnitude ($M_{\rm V}$) space (that is size-luminosity plane) for
$z_{\rm rei}=8.5$. The grey lines show lines of constant surface brightness
$\mu_{\rm V}$. For comparison, we plot observed MW confirmed/candidate dwarf
galaxies (stars; from Drlica-Wagner et al. 2020), probable ultra-faint star
clusters (red squares; Table 5) and unconfirmed systems (red cirles) as DES
J0225+0304 from Luque et al. (2017) and DELVE 2 from Cerny et al. (2021a). We
also add the recently discovered MW dwarf galaxies Eridanus IV (Cerny et al.,
2021b) and Pegasus IV (Cerny et al., 2022). We also plot the data from other
Local Volume dwarf galaxies (e.g., Sharina et al., 2008; McConnachie, 2012;
Romanowsky et al., 2016; Carlin et al., 2016; Crnojević et al., 2019;
Martínez-Delgado et al., 2021; Mutlu-Pakdil et al., 2022; Carlin et al., 2021;
Sand et al., 2022). (Right) The 2D distribution of our model dwarf galaxies in
absolute V-band magnitude ($M_{\rm V}$) and surface brightness ($\mu_{\rm V}$)
space (that is surface brightness - luminosity plane). The grey lines show
lines of constant physical half-light radius $r_{1/2}$. Note that this plot
includes all confirmed/candidate MW dwarf galaxies and not just systems used
in our analysis (e.g., above plot includes Antlia II). We use the updated
measurements for Antlia II from Ji et al. (2021).

Figure 11 shows evolution of $R_{\rm 50\%,survey}$ for each survey over a
period of $\sim 2.5$ Gyr. The figure shows that as the LMC approaches its
pericentre, $R_{\rm 50\%,survey}$ predicted for the DES survey decreases (i.e.
concentration of the satellite radial distribution increases), reaching
minimum around the epoch of the LMC pericentre passage. However, no
significant variation around LMC pericentre is observed in the predictions for
the PS1 and SDSS + Classical samples.

Figure 11 also shows that $R_{\rm 50\%,survey}$ for DES increases around
$\tau-\tau_{0}\sim-1.3$ when the LMC is far from its pericentre resulting in
an extended satellite radial distribution. Over a very short period of $\sim
0.5$ Gyr, the $R_{\rm 50\%,survey}$ for DES satellite radial distribution is
doubled. Such behavior of satellite radial distributions based on position of
a massive subhalo (like LMC) within a host halo has been observed in other
studies (see, e.g., Figure 13 in D’Souza & Bell 2021).

###  3.4 Dwarf Galaxy Sizes

Size of a galaxy is one of its key properties. Among other things, size, for a
given luminosity, defines galaxy surface brightness and detectability, which
makes it important for our forward-modelling of UFD population. We therefore
present comparisons of the size distributions of our model galaxies to
observations in this section.

Detailed description of the size modelling in our framework is presented in
Kravtsov & Manwadkar (2022). Briefly, the model assumes that the gaseous disk
scale length $R_{\rm d}$ at any epoch $z$ is proportional to $R_{\rm
200c}(z)$: $R_{\rm d}=\chi_{\rm d}R_{\rm 200c}$. The extent of the
distribution of newly formed stars is computed from the radial surface density
of H2, which is used to compute radial star formation surface density profile.
Finally, the size of the overall stellar distribution is computed from the
total stellar mass profile formed during evolution assuming that there is no
significant radial re-distribution of stars during evolution.

Figure 13: Predicted half-light radius $r_{1/2}$ distribution for satellites
within $d_{\rm helio}<300$ kpc in DES (left), PS1 (centre) and SDSS +
Classical (right) samples as a function of reionization redshift $z_{\rm
rei}=7.5,8.5,9.5,10.5$. The distributions are normalized to the total number
of galaxies with sizes smaller than 5000 pc. The color scheme and shaded
regions carry the same meaning as in Figure 4. Note that we are only including
mock and observed galaxies with $r_{\rm 1/2}>10$ pc in the above forward
modelling.

The proportionality constant for each satellite is set by $\chi_{\rm
d}=0.06\times 10^{\mathcal{N}(0,0.25)}$, where $\mathcal{N}(0,0.25)$ is normal
probability distribution with zero mean and rms of $0.25$. Both the constant
of $0.06$ and scatter of $0.25$ dex were chosen to approximately match the
normalization of the size-stellar mass relation of dwarf galaxies in the Local
Volume (Kravtsov & Manwadkar, 2022).

Figure 12 shows the 2D distribution of our model dwarf galaxies in the
size–luminosity ($r_{\rm 1/2}-M_{\rm V}$) and surface brightness - luminosity
($\mu_{\rm V}-M_{\rm V}$) planes along with stellar systems (confirmed and
candidate dwarf galaxies, star clusters, and unconfirmed systems) detected
around Milky Way in recent surveys. For comparison, we also show dwarf
galaxies around other Local Group galaxies. The model prediction in Figure 12
is shown for $z_{\rm rei}=8.5$. However, we find that distribution is not very
sensitive to $z_{\rm rei}$.

The figure shows that sizes of model galaxies reproduces the wide range of
observed sizes of dwarf galaxies. Remarkably, the model also reproduces the
scaling of size and surface brightness with $M_{V}$ and subtle features in the
scaling of observed galaxies, although there was no tuning of the model to
match these features.

For example, there is a distinct break in the $r_{\rm 1/2}-M_{\rm V}$ relation
in observed dwarf galaxies at $M_{\rm V}\sim-6$. A similar break at $M_{\rm
V}\sim-6$ is also seen in the model galaxy distribution. This break occurs at
the luminosity below which dwarf galaxies become sensitive to reionization
(see Section 3.1.2). Such a break in the $r_{\rm 1/2}-M_{\rm V}$ relation is
expected due to reionization in the context of our size modelling because
dwarf galaxies with $M_{\rm V}\gtrsim-6$ have formed most of their stellar
population before reionization. Their size therefore reflects the small virial
radius of their parent haloes at these early epochs. This effect is likely
enhanced due to the fact that molecular hydrogren fraction is high only at the
central regions of these galaxies and this makes stellar distribution even
more compact.

Scatter of the model galaxy sizes is due to the scatter in $\chi_{\rm d}$
described above, and scatter of size and luminosity arising from the scatter
in the mass assembly history of their haloes in the presence of gas accretion
suppression due to UV heating. The overall scatter is $\approx 0.5-0.7$ dex,
consistent with observational constraints for MW satellites (Nadler et al.,
2020). The fact that the accretion rate around $M_{\rm c}(z)$ changes smoothly
contributes to the scatter. We find that in the scenario where accretion is
shut down immediately once halo mass drops below $M_{\rm c}(z)$ (as assumed,
e.g., by Bose et al., 2018), produces a significantly smaller scatter in the
$r_{\rm 1/2}-M_{\rm V}$ plane and smaller than the scatter exhibited by
observed galaxies in Figure 12.

Note, however, that our dwarf galaxy model does not include physical processes
that potentially can lead to increase of sizes due to feedback or tidal
heating, as is likely for galaxies like Antlia II (e.g., Ji et al., 2021).

Figure 14: The predicted luminosity function, radial distribution and half-
light radius $r_{\rm 1/2}$ distribution for clusters with $2\text{ pc}<r_{\rm
1/2}<10$ pc, $\mu_{\rm V}>24$ mag/arcsec2 and $M_{\rm V}<0$ that are detected
with high significance in DES and PS1. The cuts of $2\text{ pc}<r_{\rm 1/2}$
and $\mu_{\rm V}>24$ mag/arcsec2 are applied to only consider observed systems
that occupy the same region in the size-luminosity plane as our model dwarf
galaxies. The cut of $M_{\rm V}<0$ is applied to ensure that our model
predictions are not affected by numerical effects present in low-mass
subhaloes. The observed systems considered here are denoted in Table 5 with
the $\dagger$ symbol.

To make the comparison more quantitative, Figure 13 shows the projected half-
light radius $r_{\rm 1/2}$ distribution for satellites predicted by our model
for the DES, PS1, and SDSS+Classical samples for models with different $z_{\rm
rei}$ and compares the with observed $r_{\rm 1/2}$ distribution of galaxies in
these samples. As before, here we only include model galaxies with $r_{\rm
1/2}>10$ pc – the cut used to define the observational samples we are
comparing with. Figure 13 shows good agreement between our model predictions
and observed satellite size distributions. This agreement is non-trivial as
the size model normalization was calibrated to match the size-stellar mass
relation of nearby dwarf galaxies, not the distribution of sizes of the MW
satellites.

Interestingly, our model predicts that a sizable population of ultra-faint
dwarf galaxies with $M_{\rm V}\gtrsim-2$ and $5\text{ pc}<r_{\rm 1/2}<20$ pc.
A significant fraction of these are in the ultra-compact area of $r_{\rm
1/2}<10$ pc. This part of the $r_{\rm 1/2}-M_{\rm V}$ plane is occupied
predominantly by star clusters, but our model implies that some of the
currently unconfirmed stellar systems in this regime may be dwarf galaxies. We
further discuss potential implications of such population of small-size dwarfs
and its plausibility in Section 4.1.

##  4 Discussion

###  4.1 The smallest dwarf galaxies

Galaxies are usually defined as gravitationally bound stellar systems with
stellar and/or gas kinematics that cannot be explained by a combination of
baryons and Newton’s law of gravity (e.g., Willman & Strader, 2012). Such
definition distinguishes between star clusters and galaxies embedded in a dark
matter halo. This distinction is important because many of the faintest
stellar systems detected around MW lie in the region of the size-luminosity
plane ($L_{V}\lesssim 10^{5}\,L_{\odot}$, $r_{1/2}\lesssim 10\text{pc}$) that
is occupied both by star clusters and the UFD galaxies (see Figure 12).

Thus, to avoid confusion with star clusters, galaxy samples are often defined
using a lower limit of $r_{\rm 1/2}=10$ pc (e.g., Conn et al., 2018; Drlica-
Wagner et al., 2020). We applied such cut to our model galaxy sample in the
comparisons presented in the previous section. Nevertheless, our model
predicts that UFDs have sizes ranging continuously down to $r_{\rm 1/2}\sim 2$
pc (Figure 12). The implication is that a significant fraction of stellar
systems with $M_{\rm V}\gtrsim-2$ and $r_{\rm 1/2}<10$ pc can be dwarf
galaxies.

Figure 14 shows the predicted luminosity, distance and half-light radius
$r_{\rm 1/2}$ distribution for compact, faint stellar objects with $r_{\rm
1/2}<10$ pc and $M_{\rm V}<0$ detected with high significance in DES and PS1.
We compare our predictions with observations of objects with $2\text{
pc}<r_{\rm 1/2}<10$ pc, $\mu_{\rm V}>24$ mag/arcsec2 and $M_{\rm V}<0$. The
cuts of $2\text{ pc}<r_{\rm 1/2}<10$ pc and $\mu_{\rm V}>24$ mag/arcsec2 are
applied to only consider objects that occupy the same region in size-
luminosity plane as our model dwarf galaxies (refer to Figure 12). The cut of
$M_{\rm V}<0$ is applied to ensure that our model predictions are not affected
by numerical effects present in low-mass subhaloes. The observed sample used
in this comparison is denoted by $\dagger$ symbol in Table 5.

Figure 14 shows that the predicted abundance and distributions of distances
and sizes of such compact systems predicted in our model are in agreement with
observations. This agreement is non-trivial as our model predictions are
calibrated of dwarf galaxies at $r_{\rm 1/2}>10$ pc and only on normalization
of sizes at a given stellar mass. There is thus no a priori expectation of
agreement between observations and theoretical predictions for $r_{\rm
1/2}<10$ pc regime. Although the number of systems in this comparison is
small, the agreement indicates that such compact, faint stellar systems might
indeed be dwarf galaxies formed within small-mass halos.

We note that such compact faint stellar systems do form in recent high-
resolution of cosmological simulations of dwarf galaxies (e.g., see the dashed
line in Fig. 5 of Agertz et al., 2020), which reflects the small mass and
virial radius of their host halo at high $z$ where most of the stars form.
Observationally, measurements of velocity and metallicity dispersions of
compact low-luminosity stellar systems detected in the recent surveys would be
needed to confirm this possibility.

###  4.2 The Predicted MW Satellite Luminosity Function within 300 kpc

Figure 15: (Left) The predicted MW satellite luminosity function within
$d_{\rm helio}<300$ kpc for satellites with effective half light radii
$r_{1/2}>10$ pc. For comparison, we show predictions for the total number of
luminous MW satellites in other studies, namely Nadler et al. (2020), Jethwa
et al. (2018), Newton et al. (2018), Hargis et al. (2014) and Tollerud et al.
(2008). The solid black line shows the current observed MW satellite
luminosity function. (Right) The predicted surface brightness $\mu_{\rm V}$
distribution for MW satellites with $M_{\rm V}<0$ and $r_{1/2}>10$ pc. The
vertical dashed line indicates the predicted detection limit for LSST Y1
survey. The solid black line shows the current observed MW satellite surface
brightness distribution.

Figure 15 shows predictions of our model for the total MW satellite luminosity
function and surface brightness function within a heliocentric distance of 300
kpc and for satellites with effective half light radii $r_{1/2}>10$ pc. The
68% and 95% contours in Figure 15 are due to a combination of host-to-host
scatter and the inferred range of $z_{\rm rei}$ values from current data (see
Section 3.1.3). The model predicts $440^{+201}_{-147}$ (440${}^{+328}_{-242}$)
with 68% (95%) range satellites with $M_{\rm V}<0$, i.e. roughly an order of
magnitude larger than the currently known number. The right panel of Figure 15
shows that $>$99% of these satellites have surface brightness $\mu_{\rm
V}<32\,\rm mag\,arcsec^{-2}$ and thus should be detectable in the LSST survey
with the Vera C. Rubin Observatory (subject to its footprint constraints and
other selection effects) if background contamination (e.g., faint galaxies)
will not be overwhelming.

For comparison, the figure shows estimates of the expected satellite abundance
within this volume from several previous studies. Our estimate for the total
number of luminous satellites at the faint end $M_{V}>-3$ is a factor of $2-3$
higher then recent estimates of Newton et al. (2018) and Nadler et al. (2020),
but the difference is within 2 $\sigma$. This difference is mainly due to the
flattening of the $M_{\rm h,peak}-M_{V}$ relation in the UFD regime due to UV
heating (see Section 3.2.1, and is partly due to increased scatter of
luminosities at smaller masses, which increases the number of small-mass
subhalos that get scattered into observable range. We discuss and evaluate
this difference in the context of current estimates of the UFD satellite
abundance in the HSC-SSP survey in Section 4.3. Our estimates of the abundance
of dwarf satellites with $M_{\rm V}\lesssim-3$, on the other hand, are
consistent with previous studies.

The estimate of Jethwa et al. (2018) shown in Figure 15 are for their model
assuming 100% halo occupancy and fiducial power law SHMR model. The occupation
fraction and scatter assumed in their study are quite different from the
results of our model, but the overall UFD abundance they estimate is
consistent with our model results.

The prediction for the number of satellites with $M_{\rm V}<-2.7$ within 300
kpc from Tollerud et al. (2008) is the most discrepant from our model and
results of other studies. The reasons for this large discrepancy are twofold.
First, estimates from Tollerud et al. (2008) are based on the SDSS DR5
observations of 11 MW dwarf galaxies within a footprint of 8000 deg2, which
was rescaled to the entire sky. However, only 3 new satellites were discovered
within an additional 6500 deg2 of data in the SDSS DR9. Second, the SDSS
selection function used in Tollerud et al. (2008) is shallower (less
sensitive) compared to other SDSS selection functions used (e.g., Walsh et
al., 2009), as was also pointed out by Newton et al. (2018). This biases
satellite census prediction from Tollerud et al. (2008) to larger values.

In fact, our discussion on halo occupation fraction (Section 3.2.2) and
results from Jethwa et al. (2018) and Nadler et al. (2020) show that the
observed abundance of UFDs in the DES/PS1 surveys cannot be explained if star
formation is suppressed in haloes with $M_{\rm peak}\lesssim 3\times
10^{8}\,M_{\odot}$ (see also Graus et al., 2019). Recent numerical studies
show that haloes with masses as low as $M_{\rm peak}\lesssim 3\times
10^{7}\,M_{\odot}$ do form stars prior to reionization (e.g., Ricotti &
Gnedin, 2005; Bovill & Ricotti, 2009; Bland-Hawthorn et al., 2015, see also
Côté et al. 2018 for a compilation of recent galaxy formation results).
Suppression of star formation in the haloes with masses below “the atomic
cooling limit” corresponding to $T_{\rm vir}\sim 10^{4}$ K of $V_{\rm
circ}\lesssim 15$ km/s ($M_{\rm peak}\lesssim 3\times 10^{8}\,M_{\odot}$) is
thus strongly disfavored both on theoretical and observational grounds.

###  4.3 Prospects for current and future surveys

Figure 16: Our forecasts for satellite luminosity functions within $d_{\rm
helio}<300$ kpc in current HSC-SSP footprint (top left), total HSC-SSP
footprint (top right), DELVE-WIDE (bottom left) and LSST Y1 (bottom right).
Similar to before, these are satellite galaxies with effective half-light
radii $r_{\rm 1/2}>10$ pc and heliocentric distance $d_{\rm helio}<300$kpc.
The different colors denote different reionization histories and carry same
meaning as previous figures (Figure 4 for instance). For predictions for
current HSC-SSP observations (top left), we show comparison with currently
observed HSC-SSP satellite luminosity function (solid black line) and
agreement is observed in the total number of satellites observed. For total
HSC-SSP survey (top right), we show the satellite luminosity function of
satellites currently observed in HSC-SSP and other known satellites that fall
in its total footprint. For predictions for the ongoing DELVE-WIDE (bottom
left) survey and future LSST Y1 (bottom right) survey, we show the satellite
luminosity function of known satellites that fall within respective footprint
and can be detected by the respective survey.

Given that our model reproduces satellite luminosity functions in the DES, PS1
and SDSS surveys, we can use our model to forecast LFs that will be measured
in ongoing surveys, such as HSC-SSP (Homma et al., 2018, 2019) and DELVE-WIDE
(Drlica-Wagner et al., 2021), and the near-future LSST survey that will use
Vera C. Rubin Observatory (Ivezić et al., 2019). We present these forecasts
for different redshifts of reionization $z_{\rm rei}=7.5,8.5,9.5,10.5$ and for
our inferred range for $z_{\rm rei}$ (Section 3.1.3) in Figure 16. The black
solid lines in Figure 16 show the current observed luminosity function in the
corresponding survey footprint and accounting for the survey selection
function. This gives an indication of how many satellites are currently known
and how many more we expect to find in each of these surveys.

The top left panel in Figure 16 shows that our forecast for the observed
satellite luminosity function in the current HSC-SSP footprint of
$\sim$676$\,\rm deg^{2}$ is consistent to within 1$\sigma$ for the total
observed abundance of UFDs, but is somewhat lower than the abundance of
galaxies with $M_{V}\lesssim-2$ (still consistent at the 2$\sigma$ level).
Nadler et al. (2020) does a different version of this calculation where they
compute the total number of satellites observable in HSC-SSP, but not
observable in DES and/or PS1. While the observed number of satellites for such
a calculation is 4, Nadler et al. (2020) predicts a count of $1.75\pm 0.6$
(cf. also Homma et al., 2019). Performing the same calculation using our
model, we predict a count of $6^{+3}_{-2}$. Our higher prediction is
reflective of the higher occupation fraction $f(M_{V}<0)$ in haloes of $M_{\rm
peak}\lesssim 10^{8}\,M_{\odot}$ compared to Nadler et al. (2020). Assuming
Poisson statistics, the uncertainty in the observed satellite count is 2.
Thus, both model predictions are consistent with observations at the 1$\sigma$
level. Additional satellite observations in HSC-SSP will be helpful in
distinguishing between the different halo occupation fraction scenarios.

The top right panel in Figure 16 shows our forecast for the observed satellite
luminosity function in the entire HSC-SSP footprint of $\sim$1400$\,\rm
deg^{2}$. The solid black line is the LF for satellites detected in the survey
so far and other satellites that fall in the total HSC-SSP footprint, namely,
Pisces II, Leo V and DESJ0225+0304. The total number of satellites with
$M_{\rm V}<0$ and $r_{\rm 1/2}>10$ pc predicted by our model for the total
HSC-SSP survey is $13^{+7}_{-5}$ at the 68% level of confidence.

The bottom left panel in Figure 16 shows our forecast for observed satellite
luminosity function in the ongoing DELVE-WIDE survey and a comparison with the
current satellite LF in the DELVE-WIDE footprint denoted by solid black line.
Specifically, we Monte Carlo sample the detection probabilities using the
DELVE-WIDE selection function to get the median observed LF, which is plotted
as the black solid line. The total number of satellites with $M_{\rm V}<0$ and
$r_{\rm 1/2}>10$ pc expected in the final DELVE-WIDE survey according to our
model is $64^{+17}_{-13}$ with 68% confidence compared to 30 currently known
in the footprint with detection probability greater than 0.5. Most of the new
expected satellites are in the UFD regime ($M_{V}\gtrsim-4$).

The bottom right panel in Figure 16 shows our model forecast for the satellite
LF in the LSST Year 1 survey. As discussed in Section 2.6, we assume that LSST
Y1 selection function will be similar to that of the HSC-SSP, that is, it will
detect all satellites with $M_{\rm V}<0$, $r_{\rm 1/2}>10$ pc, $d_{\rm
helio}<300$ kpc and $\mu_{\rm V}<32\,\rm mag\,arcsec^{-2}$ The solid black
line is LF for all known satellites satisfying above detection criterion that
fall within LSST footprint. The total number of satellites with $M_{\rm V}<0$
and $r_{\rm 1/2}>10$ pc expected in the LSST Year 1 survey is
$273^{+119}_{-92}$ at the 68% confidence level. We thus predict that LSST Year
1 survey will bring our total observed census to $\approx 60\%$ of the MW
satellite population.

We caution that the above predictions for the HSC-SSP and LSST Year 1 surveys
are subject to systematic uncertainty due to uncertain selection function at
the extremely faint magnitude limit expected for these surveys. These
systematic uncertainties also make it very difficult to forecast the results
of the LSST survey based on 10 years of co-added data, which is expected to
have a 5-$\sigma$ point source magnitude of $r_{\rm lim}=27.5$. At this depth,
the number of unresolved background galaxies acting as contamination far
exceeds the number of stars (e.g., Willman, 2010), while their colors are
similar to those of the main sequence turn off stars in the local dwarf
galaxies. This can make detection of stellar over-densities tricky or
impossible compared to current surveys of shallower depth. Thus, the detection
efficiency of LSST 10 Year co-added survey will require detailed analyses of
the survey data and modelling of the corresponding affects. Thus, we do not
attempt to make predictions for LSST 10 Year co-added in this work.

Comparison of the LF predicted for the DELVE-WIDE and LSST surveys shows
difference in the LF shape at the faintest luminosities: the cumulative LF in
the former flattens, while LF in the LSST survey continues to increase as a
fairly steep power law down to $M_{V}=0$. This difference is due to different
surface brightness limits of these surveys.

Figure 17: Predicted MW satellite luminosity functions ($d_{\rm helio}<300$
kpc, $r_{\rm 1/2}>10$ pc) for different cuts on surface brightness $\mu_{\rm
V}$ (namely $\mu_{\rm V}<26,29,32\ \rm mag\,arcsec^{-2}$) and heliocentric
distance $d_{\rm helio}$ (namely $d_{\rm helio}<60,100,200$ kpc). For
comparison, we show the current observed satellite luminosity functions as
solid black lines for the same cuts on surface brightness or heliocentric
distance. The dashed line in each panel shows the model predicted median total
MW satellite luminosity without any cuts on surface brightness or heliocentric
distance. In other words, the dashed line is the same as the solid median line
from left panel in Figure 15.

Figure 17 show the satellite luminosity functions predicted by our model for
three different cuts in surface brightness or heliocentric distance. In all
cases, we apply a cut of $|b|>10^{\circ}$ to mimic the missing coverage around
the Galactic plane. The solid black lines in each panel show the current
observed satellite luminosity with those same respective cuts. The dashed teal
colored line shows the median predicted total MW satellite luminosity
function. As one might expect, the satellite luminosity function completeness
increases as a function of decreasing heliocentric distance. We appear to be
essentially complete at $d_{\rm helio}<60$ kpc except for $M_{\rm
V}\gtrsim-2$. As shown above, ongoing surveys like DELVE-WIDE and upcoming
surveys like LSST will help detect a significant fraction of these satellites.

The LF is flattened at faintest luminosities for the $\mu_{\rm V}<26\,\rm
mag\,arcsec^{-2}$ limit but steepens for the samples with fainter limits. One
might expect to be complete at the $\mu_{\rm V}<26\,\rm mag\,arcsec^{-2}$
limit after accounting for galactic plane incompleteness, however, our model
predicts a surplus of satellites relative to observations at $-2\gtrsim M_{\rm
V}\gtrsim-8$ with $\mu_{\rm V}<26\,\rm mag\,arcsec^{-2}$. For instance, the
DES and PS1 surveys are expected to be complete at $M_{\rm V}\gtrsim-4$ and
$\mu_{\rm V}<26$ at $d_{\rm helio}<300$ kpc. However, the difference could be
due to a lack of complete coverage at $|b|>10^{\circ}$ with depths equivalent
to DES/PS1. Uncertainties in size modelling, especially around the break in
the $r_{\rm 1/2}-M_{\rm V}$ could contribute to this as well.

###  4.4 Modelling Caveats

In the preceding sections we demonstrated that our model is able to correctly
and non-trivially reproduce a number of properties of the observed dwarf
galaxy population of the Milky Way and their distributions (see also
additional comparisons in Kravtsov & Manwadkar 2022). Nevertheless, the model
has a number of uncertainties or assumptions that might influence our
determination of MW satellite statistics and have to be discussed as caveats,
which is the focus of this section.

####  4.4.1 Galactic Outflow Modelling

Feedback-driven outflows are undoubtedly the key process shaping properties of
dwarf galaxies (e.g., Collins & Read, 2022). In the GRUMPY model the outflows
are parameterized using the stellar mass dependent mass loading factor
(Kravtsov & Manwadkar, 2022). In this work, we slightly modify the galactic
outflow parameters to be $\eta_{\rm norm}=1.8$ and $\eta_{\rm p}=-0.45$
compared to the values used in Kravtsov & Manwadkar (2022) of $\eta_{\rm
norm}=3.6$ and $\eta_{\rm p}=-0.35$. As discussed in Section 2.5, these
modifications are sensible and are done to better match the bright end of the
MW satellite luminosity function. The galactic outflow parameters affect the
bright end of the satellite luminosity function as they affect the high
stellar mass end of the stellar mass - halo mass relation as seen in Figure 3
of Kravtsov & Manwadkar (2022). Thus, different galactic outflow parameters
will result in a different satellite counts as well. Furthermore, Figure 7 in
Kravtsov & Manwadkar (2022) showed that agreement with observed stellar mass -
metallicity relations can be brought about with a range of galactic outflow
scaling relations or values of the wind metallicity $\eta_{\rm w}$ factor.
Thus, to account for effect of scatter/uncertainties in galactic outflows
parameters on our model predictions, we will first have to properly constrain
these outflow parameters using a likelihood based analysis on the bright end
of satellite luminosity function and stellar mass - metallicity relations
together. Thus, our uncertainties on predicted number of satellites and
$z_{\rm rei}$ are likely underestimated. We leave this investigation for
future work.

####  4.4.2 Uncertainties in Galaxy Size Modelling

Our dwarf galaxy size model is able to accurately reproduce the observed dwarf
galaxy statistics (for instance Figure 13). The basic assumption of our size
model is that the gaseous disk size is proportional to $R_{\rm 200c}$ with
some scatter. The linear relation between half-mass radius of stellar
distribution and $R_{\rm 200c}$ has been shown to exist for galaxies with
$r_{\rm 1/2}\gtrsim 100-200$ pc (e.g., Kravtsov, 2013), and a similar relation
for the disk scale length of gaseous disk and $R_{\rm 200c}$ for $\sim
L_{\star}$ galaxies. Such a relation between $R_{\rm 200c}$ and gaseous disk
size is also theoretically motivated (e.g., Fall & Efstathiou, 1980; Ryden &
Gunn, 1987; Mo et al., 1998) because the angular momentum acquired by haloes
during their collapse will influence the collapse of baryons onto the halo,
setting the galaxy size. Nevertheless, in our model this assumption is
extrapolated into the UFD regime with $r_{\rm 1/2}\lesssim 50$ pc.
Furthermore, the scatter assumed by our size model in $\chi_{\rm d}$ is 0.25
dex comparable to the scatter expected for the specific angular momentum,
while larger scatter might be expected in our size model due to other factors
such as tidal heating (e.g., Errani & Peñarrubia, 2020).

Another uncertainty stems from systematic uncertainties in the $\rm H_{2}$
modelling, as stars in our model form from the molecular gas and its extent
within the gas disk thus defines the extent of newly formed stellar
population. As shown in Appendix B in Kravtsov & Manwadkar (2022), the current
calibrations of the molecular gas models have significant uncertainty in the
low-metallicity regime.

Finally, in our model the size of the gas disk is assumed to be proportional
to $R_{\rm 200c}$ even when gas accretion is suppressed by UV heating. Thus,
the size of the gaseous disk continues to grow even when its mass does not
grow. This, however, does not result in increase in stellar half mass radius
because the extent of the molecular gas would actually be either constant due
to increasing metallicity or would shrink due to decreasing gas surface
density. We find that only a small subset of model dwarf galaxies undergoes
evolution in this regime which results in very faint surface brightness.

####  4.4.3 Merger Histories, Subhalo Disruption & MW’s peculiarities

As described in Section 2, for comparisons with observations we select host
haloes that have realistic LMC analogs and there are only two hosts of this
kind in the Caterpillar suite. This does not allow us to quantify the host-to-
host scatter and median or mean statistics of the satellite population. We
also cannot impose additional selection criteria using other MW properties,
such as presence of the SMC, the existence of a massive M31-like neighbor,
details of the MW merger history etc. To our current knowledge, the MW
underwent 3 signifcant merger events: i) Gaia-Enceladus-Sausage (GES), ii)
Sagittarius and iii) LMC. In fact, Bose et al. (2020) finds the incidence of
hosts with a GES-like and LMC-like mergers to be $\sim 3\%$. Reproducing such
properties of the Milky Way will require a much larger suite of the high-
resolution zoom-in simulations.

Additionally, the two hosts we select have rather similar merger histories, as
can be seen in the right panel of Figure 2. It is not clear whether this is a
coincidence of this is due to the LMC analog criterion and using two host
haloes with similar merger histories in our analysis could possibly bias our
results somewhat.

Figure 18: (Left 4 Panels): The cumulative subhalo radial distributions in
different $M_{\rm peak}$ ranges for the Caterpillar hosts at the different
resolution levels of LX12, L13 and LX14. The dashed black lines in each of
these correspond to the mean NFW profile for all the hosts. Their is
negligible difference in host concentration between different resolution
levels so we only show one line for mean NFW profile of host. The different
line styles correspond to different resolution suites. Note that the $5\times
10^{6}<M_{\rm peak}<10^{7}$ line for Caterpillar LX12 is not shown as it has
no subhaloes in that range. (Right Panel) The radial bias function between the
subhalo radial distribution and host NFW profile shown for different $M_{\rm
peak}$ ranges. The shading of the lines carries the same meaning as the
shading in the left 4 panels of the figure. The radial bias function (denoted
by ‘NFW Bias’) is defined as the ratio of $\rm dN/dx^{3}|_{\rm subhalo}$ and
$\rm dN/dx^{3}|_{\rm NFW}$ where $\rm dN/dx^{3}$ is the number of subhaloes
per unit shell volume as a function of $x=r/r_{\rm vir}$. The bias function
will be equal to one if it is unbiased with respective to the host NFW
profile. The above plots show that the effects of withering and artificial
disruption are prominent for small mass subhaloes as their radial distribution
gets increasingly shallower and biased away from the host NFW profile.

Details of mass assembly history of the host halo possibly translate into
subhalo disruption modelling uncertainties. As we use the Caterpillar
simulation suites, which do not model the central galactic disk potential, we
apply Nadler et al. (2018)’s subhalo disruption model to these hosts. The
subhalo disruption model from Nadler et al. (2018) was trained on only two
zoom-in simulations with very similar merger histories and so also has
systematic uncertainty. Furthermore, when using Nadler et al. (2018)’s subhalo
disruption model, we are implicitly assuming that all of the subhalo
disruption is due to MW’s central galactic disk. However, Jahn et al. (2019)
showed that LMC host halo depletes $\sim 30\%$ of its substructure within
$0.2r_{\rm 200m}$ prior to accretion onto MW-size host. This effect is small
compared to the disruption by the MW mass host and disk (Jahn et al. 2019) and
thus does not affect our results and conclusions, but such additional
disruption of LMC satellites by the LMC itself may need to be included in the
future models if they focus on the predictions in surveys around LMC.

Finally, recent studies of Webb & Bovy (2020) and Green et al. (2022) conclude
that effect of disk on subhalo disruption is much smaller than derived from
previous simulations due to insufficient resolution of the latter. This effect
is thus quite uncertain, which introduces uncertainty in the halo occupation
distribution of observable galaxies inferred within the context of our model.
An additional uncertainty is due to numerical effects affecting subhalo
disruption and their radial distribution and abundance discussed in the next
subsection.

####  4.4.4 Withering, Artificial Disruption & Orphan Galaxies

The predictions for the luminosity function and radial distribution of
satellite galaxies are based on the subhalo population down to masses of $\sim
10^{7}\,\,M_{\odot}$ in the Caterpillar simulations used in this study and can
potentially be subject to numerical effects. For example, Bose et al. (2020)
have recently argued that radial distribution of satellite galaxies cannot be
reproduced without accounting for numerical disruption of subhaloes by
modelling disrupted subhaloes with an “orphan galaxy” model. As shown in
Section 3.1 and Section 3.3, our model predicts the luminosity function and
radial distribution of satellites in reasonable agreement with observations.
The question is thus to what extent numerical disruption affects the subhalo
population on which these predictions are based.

In a recent study, Green et al. (2021) presented a comprehensive model for the
evolution of subhalo populations using the SatGen model of Jiang et al.
(2021). This model accounts for the effects of tidal mass loss and potential
loss of subhaloes when their bound mass becomes smaller than the resolution
limit of simulation and/or halo finder (“withering”) and premature disruption
of subhaloes due to numerical effects (“artificial disruption”). Green et al.
(2021) showed that withering has a large effect on the radial distribution of
subhaloes and the amplitude of their mass function, while artificial
disruption has a much smaller but still noticeable effect. They also showed
that when withering is corrected for, the radial distribution of subhaloes is
expected to be close to the radial distribution of dark matter in their host
halo.

Table 3: The Spearman correlation coefficient $\rho$ of $r/r_{\rm vir}$ with $M_{\rm acc}/M_{\rm host}$ and $\log_{10}(1+z_{\rm acc})$ for subhaloes with different $M_{\rm peak}$ cuts and in the different Caterpillar resolution suites. For comparison, in the rightmost two columns, we show these correlation coefficients from Bolshoi simulations for all subhaloes and their predictions from Green et al. (2021). ‘G21(W+D)’ denotes predicted value in the presence of W+D and ‘G21(All)’ denotes the same value but in the absence of W+D. Note that we do not quote a correlation coefficient in the [$5\times 10^{6},10^{7}$] bin for LX12 as no subhaloes are found there. The standard deviation $\sigma$ on all these correlation quantities are negligible ($\sigma<0.01$). The exception is that we have $\sigma\approx 0.02$ for the $[10^{9},)$ bin at all resolutions and $[10^{7},10^{8}]$ bin for LX12 due to low subhalo statistics. This standard deviation is computed using the bootstrap method. Correlation Quantities |  $M_{\rm peak}(\,M_{\odot})$ bins | LX12 | LX13 | LX14 |  |
---|---|---|---|---|---|---
| $[5\times 10^{6},10^{7}]$ | - | -0.07 | +0.00 | Bolshoi | -0.14
$r/r_{\rm vir}-M_{\rm acc}/M_{\rm host}$ | $[10^{7},10^{8}]$ | -0.11 | -0.11 | +0.00 | G21(W+D) | -0.15
| $[10^{8},10^{9}]$ | -0.09 | +0.01 | +0.06 | G21(All) | +0.01
| $[10^{9},)$ | -0.01 | +0.03 | +0.05 |  |
| $[5\times 10^{6},10^{7}]$ | - | -0.27 | -0.49 | Bolshoi | -0.33
$r/r_{\rm vir}-\log_{10}(1+z_{\rm acc})$ | $[10^{7},10^{8}]$ | -0.28 | -0.38 | -0.56 | G21(W+D) | -0.24
| $[10^{8},10^{9}]$ | -0.35 | -0.55 | -0.62 | G21(All) | -0.46
| $[10^{9},)$ | -0.53 | -0.60 | -0.62 |  |

We have examined population of satellites with $M_{\rm V}<0$, which
approximately corresponds to subhaloes with $M_{\rm peak}\gtrsim 4\times
10^{6}\,M_{\odot}$ according to Figure 8121212This $M_{\rm peak}$ limit
corresponds to the median halo mass at which $f(M_{\rm V}<0)<0.01$.. However,
most of the currently observed satellites (and hence our inferred statistics)
are at $M_{\rm V}\lesssim-1.5$ hosted by subhaloes of $M_{\rm peak}\gtrsim
2\times 10^{7}\,M_{\odot}$ (see Figure 7). We will thus examine the effects of
withering and artificial disruption on subhaloes of such masses. For purposes
of brevity, we will refer to ‘withering and artificial disruption’ together as
‘W+D’.

Firstly, note that among the subhaloes that accrete at the same epoch,
withering should have a larger effect on the radial distribution of $5\times
10^{6}\,M_{\odot}<M_{\rm peak}<10^{7}\,\,M_{\odot}$ subhaloes, which are close
to the mass resolution threshold, than on subhaloes with larger $M_{\rm
peak}$. Artificial disruption, on the other hand, is expected to affect all
subhaloes that lose $\gtrsim 95\%$ of their initial mass independent of that
mass (Green et al., 2021) and thus should not affect radial distribution of
subhaloes. We can thus probe potential existence of numerical effects by
examining radial distribution of subhaloes of different $M_{\rm peak}$.

Figure 18 shows the cumulative subhalo radial distributions for different
$M_{\rm peak}$ ranges in the Caterpillar LX12, LX13 and LX14 suites. For
comparison, we also show the NFW profile that matches the mean dark matter
density profile of the Caterpillar hosts. The figure shows that the radial
distribution becomes more concentrated with increasing $M_{\rm peak}$ in line
with general expectation of the effects of withering outlined above and in
agreement with the results of recent study by Grand et al. (2021, see their
Fig. 10). At the same time, the difference in the radial distributions of
subhaloes of different mass is very small for $M_{\rm
peak}>10^{8}\,\,M_{\odot}$. For $10^{7}<M_{\rm peak}/\,M_{\odot}<10^{7}$
subhaloes the difference is larger, but is still $\lesssim 1.6$ for $r\geq
0.1r_{\rm vir}$.

Furthermore, Green et al. (2021) show that effects of W+D (withering being the
dominant factor) are imprinted in correlations of subhalo properties, such as
$M_{\rm acc}/M_{\rm host}$ and $\log_{10}(1+z_{\rm acc})$ with their
normalized distance to the host centre, $r/r_{\rm vir}$. This is because W+D
effects are expected to be larger for subhaloes orbiting closer to the host
halo centre and stronger for haloes that spend more time within the host. In
particular, they showed that their SatGen model reproduces the Spearman
correlation coefficient for the correlations between $r/r_{\rm vir}$ and
$M_{\rm acc}/M_{\rm host}$ or $\log_{10}(1+z_{\rm acc})$ in the Bolshoi
simulation and show what these coefficients are expected to be if W+D were not
present.

In Table 3, we present the correlation coefficients of $r/r_{\rm vir}$ with
$M_{\rm acc}/M_{\rm host}$ and $\log_{10}(1+z_{\rm acc})$ for subhaloes with
different $M_{\rm peak}$ cuts and in different Caterpillar resolution suites.
We compare these numbers with the same correlation coefficients from SatGen
model and the Bolshoi simulation. Note that the G21(All) result applies when
numerical effects are completely absent, while G21(W+D) is for when they are
present at the level of the resolution of the Bolshoi simulation. The numbers
marked as ”Bolshoi” give the actual correlation coefficient measured for the
Bolshoi haloes by Green et al. (2021). Thus, correlation coefficients close to
G21(All) indicate that W+D effects are small, while those close to G21(W+D) or
Bolshoi indicate that these effects are significant.

We find that the $r/r_{\rm vir}-M_{\rm acc}/M_{\rm host}$ correlation for all
LX14 subhaloes with $M_{\rm peak}>10^{7}\,M_{\odot}$ is consistent with a weak
positive or zero correlation, similar to the prediction of the SatGen model
when there is no W+D (denoted by ‘G21(All)’). If it were dominant in
simulations, it would induce a negative correlation between halo mass ($M_{\rm
acc}/M_{\rm host}$) and radial distance ($r/r_{\rm vir}$), as it would remove
smaller mass haloes in the inner regions of the host, preferentially leaving
them in the outskirts. This is indeed seen in the SatGen prediction with W+D
and also for LX12 and LX13 simulations. We find similar results if we consider
$M_{\rm peak}/M_{\rm host}$ instead of $M_{\rm acc}/M_{\rm host}$.

For the $r/r_{\rm vir}-\log_{10}(1+z_{\rm acc})$ correlation, we find a
relatively strong negative correlation for subhaloes well above the mass
resolution limit, which is similar to prediction of SatGen model with no W+D.
When W+D is present, the negative correlation in $r/r_{\rm
vir}-\log_{10}(1+z_{\rm acc})$ will be weakened as these effects
preferentially remove subhaloes that were accreted earlier and that are
orbiting in the inner regions and thus would have lost more mass on average.
We indeed observe weakening of the negative correlation for the smallest mass
bins in LX 12 and LX 13 suites, where correlation coefficients become similar
to the predictions of the SatGen model when W+D is present or the value
measured in the Bolshoi simulation.

Another indication that effect of artificial disruption on the subhalo samples
used in our model is small is that there is no difference between subhalo mass
functions in the LX13 and LX14 suites in the well-resolved high-mass regime,
$\mu=M/M_{\rm host}>10^{-3}$. Given that artificial disruption does not depend
on the initial mass of subhalo, but affects subhaloes that lose $\gtrsim 95\%$
of their initial time, the time it takes for a halo to artificially disrupt
will increase with increasing halo mass. However, Figure 1 shows that this is
not the case even in the inner $r<50$ kpc region of the host halo.
Additionally, Green et al. (2021) found artificial disruption to affect the
subhalo mass function amplitude by at most $\sim$20%. So, even if artificial
disruption is present in lower mass subhaloes, it will bias our inferred
$z_{\rm rei}$ by at most $\sim$0.5 (can be inferred by raising contours in
Figure 6 by a factor of 1.2).

The main conclusion that can be drawn from the results discussed in this
section is that for a given mass bin, the numerical effects decrease as
resolution is increased. Likewise, in a given simulation resolution suite, the
numerical effects increase in smaller mass bins. Both Figure 18 and Table 3
indicate that the W+D effects in the smallest mass bin in the LX14 suite we
use in our modelling $[10^{7},10^{8}]$ are present, but their effect on the
subhalo population is relatively small. In the large mass bins these effects
are negligible.

These results imply that there is no need to include a separate model for
orphan galaxies in our analysis of LX14 simulation, as is needed when one uses
subhalo samples with masses close to the withering limit of the simulation
(e.g., Bose et al., 2020).

###  4.5 Comparison with previous models

Semi-analytic models of dwarf galaxies have been used in many studies over the
last three decades (e.g., Kauffmann et al., 1993; Somerville, 2002; Benson et
al., 2002a, b; Starkenburg et al., 2013; Lu et al., 2016, 2017; Bose et al.,
2016; Bose et al., 2018; Bose et al., 2020). Such models were used
productively to model dwarf satellites observed in the Local Group using mass
assembly tracks extracted from the modern generation of high-resolution zoom-
in $N$-body simulations (e.g., Starkenburg et al., 2013; Bose et al., 2018).

Our model uses a “regulator” like model that strip down semi-analytic
modelling down to the minimal number of equations (e.g., Finlator & Davé,
2008; Krumholz & Dekel, 2012; Lilly et al., 2013; Feldmann, 2013; Peng &
Maiolino, 2014). In particular, we exploit the fact that simulations showed
that accretion onto dwarf galaxies is predominantly in the cold-mode to avoid
explicit modelling of cold and hot gas phases within dwarf haloes. We also use
the fact that due to fairly steep stellar mass–halo mass in the dwarf regime
the mergers do not contribute significantly to stellar mass of dwarf galaxies
to avoid explicit modelling of mergers and simplify calculations.

This regulator type of models have been shown to be successful in matching
many key observational scaling relations and their evolution (see,e.g.,
Tacconi et al., 2020, for a review) and they were used for interpretation of
halo mass–stellar mass, stellar mass–metallicity correlations, and star
formation histories in the dwarf galaxy regime (Tassis et al., 2012;
Ledinauskas & Zubovas, 2018, 2020).

At the same time, our model employs refined modelling of several keys
processes informed by recent theoretical models and simulations results. For
example, we model effects of UV heating of gas during and after reionization
by modulating gas accretion onto haloes explicitly rather than using a hard
cutoff mass. We also explore different reionization redshifts, because the
Lagrangian volume of the Local Group is generally expected to reionize at a
different time than average volume in the Universe (e.g., Zhu et al., 2019).
We also use star formation modelling based on a model of molecular hydrogen,
which in turn is based on model for gas distribution and galaxy sizes
calibrated by observations. Our model also uses outflow prescription
calibrated using simulations (see Kravtsov & Manwadkar, 2022, for details).

A study by Chen et al. (2022) that presented the A-SLOTH model of dwarf
galaxies appeared after our paper was submitted. The authors used the model to
study $M_{\star}-M_{\rm h}$ relation and halo occupation distribution of
observed dwarf satellites of the Milky Way. Although some aspects of the
overall approach of A-SLOTH model is similar to GRUMPY, details of modelling
are rather different. For instance, a rather sophisticated model of both ISM
structure and star formation and of stellar feedback is used in the A-SLOTH
model with separate modelling of Population III and II stars and explicit
modelling of photo-heating by UV radiation of stars and momentum and energy
injection by supernovae. Sophistication of the modelling and the
correspondingly larger number of differential equations and parameters makes
A-SLOTH closer to the classic semi-analytic models. The approach of the GRUMPY
model, on the other hand, is to keep the number of equations to a minimum and
use simple, simulation calibrated parametrizations of key processes. As we
showed in this and previous paper, the model based on this approach can
simultaneously reproduce many key observed properties of dwarf galaxies.

##  5 Summary and Conclusions

We use a simple regulator-type model for modelling dwarf galaxies presented
Kravtsov & Manwadkar (2022) to study the luminosity function, radial
distribution, and size distribution of the dwarf galaxy satellite population
of the Milky Way. Although fairly simple, the model reproduces stellar masses,
gas mass, gas/stellar-phase metallicities, star formation rate and diversity
in star formation histories of dwarf galaxies. We forward model observations
of the model population accounting for observational biases in surveys through
their respective selection functions and current constraints on the MW halo
mass and presence of the Large Magellanic Cloud (LMC).

In particular, we compare model predictions with dwarf galaxy observations in
the DES, PS1, and SDSS surveys and make predictions for ongoing/future
surveys, HSC-SSP, DELVE-WIDE and Vera C. Rubin Observatory/LSST. We show that
our model is able to reproduce the observed MW satellite luminosity function,
radial distributions and the half-light radius $r_{\rm 1/2}$ distribution for
reasonable choices of model parameters, such as reionization redshift. We
present the halo mass-luminosity relation predicted and halo occupation
distribution for the MW satellite population.

Our main results and conclusions are:

  1. 1.

The amplitude of the ultra-faint end of the MW satellite luminosity function
is sensitive to the redshift of reionization $z_{\rm rei}$ at $M_{\rm
V}\gtrsim-6$ (in agreement with findings of Bose et al., 2018), while $M_{\rm
V}\lesssim-6$ is insensitive to $z_{\rm rei}$. Our model thus indicates that
the empirical boundary between “classical dwarf galaxies” and UFDs corresponds
to the physical boundary below which galaxies are strongly affected by
reionization, as was first suggested by Ricotti & Gnedin (2005) before the
UFDs were discovered.

  2. 2.

We constrain the connection between stellar mass $M_{\star}$ and peak halo
mass $M_{\rm peak}$ for MW satellites to be Equation 7 and show that this
relation has a clear break imprinted by reionization. This deviation from
power-law relation results in a high fraction of luminous galaxies in haloes
with $M_{\rm peak}<10^{8}\,\,M_{\odot}$.

  3. 3.

As a consequence, the model predicts a higher halo occupation fraction in low
mass haloes ($M_{\rm peak}<10^{8}\,M_{\odot}$) compared to recent studies. We
find that the the peak halo mass at which 50% of haloes host luminous
satellites is $M_{\rm peak}\sim 3.1\times 10^{7}\,\,M_{\odot}$. This is a
prediction of our model and the abundance of subhaloes in the Caterpillar
simulations used in this study and will need to be tested with future
observations. In particular, we find that such high halo occupancy in low mass
haloes in our model is helps to reproduce the number of dwarf galaxy
observations in the HSC-SSP survey.

  4. 4.

The high occupation fraction in $10^{7}-10^{8}\,\,M_{\odot}$ haloes implies
that suppression of star formation in the haloes with masses below “the atomic
cooling limit” of $V_{\rm circ}\lesssim 10$ km/s ($M_{\rm peak}\lesssim
3\times 10^{7}\,M_{\odot}$) is strongly disfavored.

  5. 5.

We demonstrate that our model based on the LX14 Caterpillar suite, which
accounts for the disruption of subhaloes due to disk using Nadler et al.
(2018) model, is able to match the observed radial distribution of the MW
satellites without modelling orphan galaxies. Under the assumption of
homogeneous reionization in the MW environment, $z_{\rm rei}$ only affects the
normalization of the satellite radial distribution but the concentration of
satellite radial distribution is insensitive to $z_{\rm rei}$.

  6. 6.

We show that our simple theoretical model can simultaneously reproduce
distribution of stellar half-mass radii of observed MW dwarf satellites and
scaling of the half-mass radius with stellar mass down to the faintest
detected galaxies. Comparison of predicted distribution of model galaxies to
observations suggests that some of the observed stellar systems with half mass
radii $<10$ pc may be low-luminosity dwarf galaxies, which can be tested with
the next generation of extremely large telescopes.

  7. 7.

We confirm findings of Nadler et al. (2020) that the presence of the Large
Magellanic Cloud (LMC) influences the observed satellite count in the DES
footprint due to its proximity to the LMC. Hosts with an LMC analog on its
first infall and near pericentre have $\sim$2 times larger number of ultra-
faint dwarf galaxies detected in DES footprint, compared to hosts with no LMC
analogs, an early accreted LMC or a LMC far from its first pericentric
approach. However, the LMC has a negligible impact on amplitude of satellite
luminosity function in PS1 and SDSS.

  8. 8.

The presence of the LMC affects the concentration of satellite radial
distributions with the effect being most significant for DES satellites. The
concentration of the satellite distance distribution increases significantly
when LMC analog passes through the pericenter of its orbit. We show that the
observed radial distribution concentration in DES strongly prefers a LMC on
first infall and near pericentre.

  9. 9.

Our fiducial model predicts a total of $440^{+201}_{-147}$ MW satellites (at
the 68% confidence level) with $r_{\rm 1/2}>10$ pc, $M_{\rm V}<0$, and $d_{\rm
helio}<300$. A significant fraction (>99%) of these satellites have surface
brightnesses $\mu_{\rm V}<32$ mag/arcsec2 and should be detectable in upcoming
surveys such as the LSST survey.

  10. 10.

We also present predictions for the number of dwarf satellites that will be
discovered in specific ongoing and future surveys: $13^{+7}_{-5}$ for the HSC-
SSP, $64^{+17}_{-13}$ for the DELVE-WIDE and $273^{+119}_{-92}$ for the LSST
Y1 (Big Sky approach). Note that these are numbers for the total count within
survey footprint and not for how many more dwarf galaxies will be discovered.

The success of this dwarf galaxy model in matching a wide variety of
properties of MW satellites including luminosity functions, radial
distributions and galaxy size distributions sets the stage for interpreting
the plethora of dwarf galaxy observations by LSST, DELVE-WIDE, HSC-SSP etc. in
upcoming years. This study shows the importance of the ultra-faint dwarf
galaxy population in not only understanding not only the host galaxy’s
reionization and merger history, but also constraining the very physical
processes that drive their formation.

Our code for producing all the figures in this paper is available at
https://github.com/kibokov/grumpy_mw_forward_model. The data required for
producing the figures is available at
https://data.mendeley.com/datasets/zmwh6wxyv3/1. The GRUMPY model pipeline is
available at https://github.com/kibokov/GRUMPY.

## Acknowledgements

We would like to thank the referee, Ethan Nadler, for detailed constructive
comments that improved presentation of this paper. We thank Alexander Ji and
the Caterpillar collaboration for providing halo tracks of the Caterpillar
simulations used in this study. We are grateful to Alex Drlica-Wagner, William
Cerny, Vasily Belokurov and Sidney Mau for useful discussions about
observations of dwarf galaxies and star clusters and to Frank van den Bosch
and Phil Mansfield for discussions about artificial halo disruption that
helped with interpretation of our results. We are also grateful to Alex
Drlica-Wagner for a thorough reading of the manuscript draft and detailed
comments that improved the paper. This work was supported by the National
Science Foundation grants AST-1714658 and AST-1911111 and NASA ATP grant
80NSSC20K0512.

Analyses presented in this paper were greatly aided by NumPy (Harris et al.
2020), matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), Jupyter
(Kluyver et al. 2016), pandas (Wes McKinney 2010), Astropy (Astropy
Collaboration et al. 2013), and Healpy (Zonca et al. 2019). This research has
made extensive use of Astrophysics Data Service (ADS) and arXiv preprint
repository. Finally, we are grateful to Grumpy Cat for helping us to carry on
through the past year.

## Data Availability

The data underlying this article will be shared on reasonable request to the
authors.

## References

  * Agertz et al. (2020) Agertz O., et al., 2020, MNRAS, 491, 1656
  * Ahn et al. (2012) Ahn C. P., et al., 2012, ApJS, 203, 21
  * Aihara et al. (2018a) Aihara H., et al., 2018a, PASJ, 70, S4
  * Aihara et al. (2018b) Aihara H., et al., 2018b, PASJ, 70, S8
  * Applebaum et al. (2021) Applebaum E., Brooks A. M., Christensen C. R., Munshi F., Quinn T. R., Shen S., Tremmel M., 2021, ApJ, 906, 96
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Balbinot et al. (2013) Balbinot E., et al., 2013, ApJ, 767, 101
  * Barkana & Loeb (1999) Barkana R., Loeb A., 1999, ApJ, 523, 54
  * Behroozi et al. (2013a) Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A., Primack J. R., 2013a, ApJ, 763, 18
  * Behroozi et al. (2013b) Behroozi P. S., Wechsler R. H., Conroy C., 2013b, ApJ, 770, 57
  * Behroozi et al. (2021) Behroozi P., Hearin A., Moster B. P., 2021, MNRAS,
  * Belokurov et al. (2007) Belokurov V., et al., 2007, ApJ, 654, 897
  * Belokurov et al. (2010) Belokurov V., et al., 2010, ApJ, 712, L103
  * Belokurov et al. (2014) Belokurov V., Irwin M. J., Koposov S. E., Evans N. W., Gonzalez-Solares E., Metcalfe N., Shanks T., 2014, MNRAS, 441, 2124
  * Benson et al. (2002a) Benson A. J., Lacey C. G., Baugh C. M., Cole S., Frenk C. S., 2002a, MNRAS, 333, 156
  * Benson et al. (2002b) Benson A. J., Frenk C. S., Lacey C. G., Baugh C. M., Cole S., 2002b, MNRAS, 333, 177
  * Besla et al. (2007) Besla G., Kallivayalil N., Hernquist L., Robertson B., Cox T. J., van der Marel R. P., Alcock C., 2007, ApJ, 668, 949
  * Bland-Hawthorn et al. (2015) Bland-Hawthorn J., Sutherland R., Webster D., 2015, ApJ, 807, 154
  * Bose et al. (2016) Bose S., Hellwing W. A., Frenk C. S., Jenkins A., Lovell M. R., Helly J. C., Li B., 2016, MNRAS, 455, 318
  * Bose et al. (2018) Bose S., Deason A. J., Frenk C. S., 2018, ApJ, 863, 123
  * Bose et al. (2020) Bose S., Deason A. J., Belokurov V., Frenk C. S., 2020, MNRAS, 495, 743
  * Bovill & Ricotti (2009) Bovill M. S., Ricotti M., 2009, ApJ, 693, 1859
  * Brown et al. (2012) Brown T. M., et al., 2012, ApJ, 753, L21
  * Brown et al. (2014) Brown T. M., et al., 2014, ApJ, 796, 91
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Bullock et al. (2000) Bullock J. S., Kravtsov A. V., Weinberg D. H., 2000, ApJ, 539, 517
  * Busha et al. (2010) Busha M. T., Alvarez M. A., Wechsler R. H., Abel T., Strigari L. E., 2010, ApJ, 710, 408
  * Cannon et al. (1977) Cannon R. D., Hawarden T. G., Tritton S. B., 1977, MNRAS, 180, 81P
  * Carlin et al. (2016) Carlin J. L., et al., 2016, ApJ, 828, L5
  * Carlin et al. (2021) Carlin J. L., et al., 2021, ApJ, 909, 211
  * Carlsten et al. (2020) Carlsten S. G., Greene J. E., Peter A. H. G., Greco J. P., Beaton R. L., 2020, ApJ, 902, 124
  * Cerny et al. (2021a) Cerny W., et al., 2021a, ApJ, 910, 18
  * Cerny et al. (2021b) Cerny W., et al., 2021b, ApJ, 920, L44
  * Cerny et al. (2022) Cerny W., et al., 2022, arXiv e-prints, p. arXiv:2203.11788
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chen et al. (2022) Chen L.-H., Magg M., Hartwig T., Glover S. C. O., Ji A. P., Klessen R. S., 2022, MNRAS, 513, 934
  * Collins & Read (2022) Collins M. L. M., Read J. I., 2022, Nature Astronomy, 6, 647
  * Conn et al. (2018) Conn B. C., Jerjen H., Kim D., Schirmer M., 2018, ApJ, 852, 68
  * Conroy & Gunn (2010) Conroy C., Gunn J. E., 2010, ApJ, 712, 833
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * Correa Magnus & Vasiliev (2022) Correa Magnus L., Vasiliev E., 2022, MNRAS, 511, 2610
  * Côté et al. (2018) Côté B., Silvia D. W., O’Shea B. W., Smith B., Wise J. H., 2018, ApJ, 859, 67
  * Crnojević et al. (2019) Crnojević D., et al., 2019, ApJ, 872, 80
  * D’Onghia & Lake (2008) D’Onghia E., Lake G., 2008, ApJ, 686, L61
  * D’Souza & Bell (2021) D’Souza R., Bell E. F., 2021, MNRAS, 504, 5270