Prospects for Near-Far Constraints on Reionization from Simulated LSST Observations of Milky Way Dwarf Galaxies and Stellar Streams

Introduction

Resolved stellar populations of the Milky Way (MW) and Local Group (LG) provide a unique window into understanding galaxy formation and evolution. Photometric and spectroscopic surveys continue to extend observations of dwarf galaxies to fainter luminosities and enable detailed reconstructions of their star-formation histories (SFHs), establishing them as effective near-field probes of the early Universe (E. Tolstoy et al. 2009; M. Boylan-Kolchin et al. 2015; D. R. Weisz et al. 2014a,c; D. R. Weisz & M. Boylan-Kolchin 2017; J. D. Simon 2019). The ancient stellar populations preserved in these

Corresponding author: Aritra Kundu

[email protected]

systems carry the fossil record of their high-redshift progenitors at $z \gtrsim6$, which are thought to have played a key role in reionizing the Universe. The forthcoming Legacy Survey of Space and Time (LSST) at the NSF-DOE Vera C. Rubin Observatory (Ž. Ivezić et al. 2019) is poised to substantially expand the census of near-field resolved stellar populations by detecting more faint satellites and low-surface-brightness stellar streams in the MW halo, and field galaxies outside its halo, offering an unprecedented opportunity to understand the early Universe through near-field cosmology. In particular, exploiting this window to probe the Epoch of Reionization (EoR) requires understanding both the connection between present-day systems and their high-$z$ progenitors and the limits of the upcoming surveys.

Reionization is a major phase transition in the history of the Universe when the ultraviolet (UV) radiation from the first stars, galaxies, and accreting black holes gradually photoionized the neutral hydrogen in the intergalactic medium (IGM) (R. Barkana & A. Loeb 2001; S. R. Furlanetto et al. 2006). Although the properties of the primary ionizing sources remain uncertain, the most likely scenario is that star-forming galaxies are the main drivers of reionization, with the process completing between $z \sim5.5$ and $z \sim6$ (P. R. Shapiro et al. 1994; P. Madau et al. 1999; C.-A. Faucher-Giguère et al. 2008; G. D. Becker & J. S. Bolton 2013; G. D. Becker et al. 2015; Y. Zhu et al. 2021; S. E. I. Bosman et al. 2022; Y. Zhu et al. 2023). The Hubble Space Telescope (HST) has directly observed reionization-era galaxies down to $M_{\mathrm{UV}} \sim-13$ in fields behind lensing clusters (S. L. Finkelstein et al. 2015; R. J. Bouwens et al. 2017; R. C. Livermore et al. 2017; H. Atek et al. 2018; R. J. Bouwens et al. 2022b), and the James Webb Space Telescope (JWST) is extending these measurements to $z \gtrsim10$ (S. L. Finkelstein et al. 2023; G. C. K. Leung et al. 2023; P. G. Pérez-González et al. 2023; S. L. Finkelstein et al. 2024; N. J. Adams et al. 2024; Y. Harikane et al. 2025) and slightly fainter magnitudes (e.g., $M_{\mathrm{UV}} \sim-12$) (I. Chemerynska et al. 2025; H. Atek et al. 2025).

Characterizing the population of the lowest-mass galaxies at high redshifts, and establishing its connection to the nature and abundance of present-day low-mass galaxies, is central to understanding galaxy formation in the past and present. A significant fraction of the ionizing photon budget at high redshift may come from galaxies beyond the reach of HST or JWST: some reionization models predict that faint yet numerous low-mass galaxies, potentially with $M_{\mathrm{UV}}$ as low as $\sim-8$, contribute the bulk of the ionizing photons (M. Kuhlen & C.-A. Faucher-Giguère 2012; J. H. Wise et al. 2014; B. E. Robertson et al. 2015; H. Atek et al. 2015; X. Ma et al. 2018a,b; S. L. Finkelstein et al. 2019; H. Atek et al. 2024; S. Mascia et al. 2024; C. Simmonds et al. 2024; A. Harshan et al. 2024; Z. Wu & A. Kravtsov 2024). These low-mass systems are not only sensitive probes of the quenching due to reionization (J. S. Bullock et al. 2000; T. M. Brown et al. 2014; H. Katz et al. 2020; Y. Wang et al. 2021; K. B. W. McQuinn et al. 2024a,b), but also of the impact of stellar feedback (J. H. Wise & T. Abel 2008; J. S. Bullock & M. Boylan-Kolchin 2017; T. Kimm et al. 2017; L. V. Sales et al. 2022) and the nature of dark matter (E. O. Nadler et al. 2021; J. Sipple et al. 2025). Furthermore, observations of faint galaxies behind lensing clusters carry significant systematic uncertainties from modeling the cluster mass distributions, intracluster light contamination, and other effects (H. Atek et al. 2018; R. J. Bouwens et al. 2022a). Even fainter sources are ever more difficult to measure directly, and many of them are expected to remain inaccessible even to JWST (M. Boylan-Kolchin et al. 2015, 2016; D. R. Weisz & M. Boylan-Kolchin 2017). Hence, indirect methods are required to constrain the faint end of the high-$z$ UVLF, which is important for quantifying the role of low-luminosity galaxies in reionization, and for understanding galaxy formation and evolution more generally.

The “near-far” approach is one such promising avenue that uses the present-day ($z = 0$) MW satellites and LG dwarf galaxies as fossil records of their $z \gtrsim6$ progenitors (M. Boylan-Kolchin et al. 2015, 2016; D. R. Weisz & M. Boylan-Kolchin 2017; D. R. Weisz et al. 2014a). The color-magnitude diagrams (CMDs) of nearby faint galaxies obtained from resolved HST/JWST photometry allow accurate reconstruction of their SFHs (D. R. Weisz et al. 2014a,b), from which the UV luminosities of their $z \gtrsim6$ progenitors can be inferred using stellar population synthesis models (D. R. Weisz et al. 2014c; D. R. Weisz & M. Boylan-Kolchin 2017). This indirect technique provides access to magnitudes as faint as $M_{\mathrm{UV}} \sim-10$, extending below the limits of direct high-$z$ observations.

To date, applications of the near-far technique to probe the faint end of the UVLF at $z \gtrsim6$ have been restricted to surviving MW satellites and intact LG dwarf galaxies (M. Boylan-Kolchin et al. 2016; D. R. Weisz & M. Boylan-Kolchin 2017; D. R. Weisz et al. 2014c; E. D. Skillman et al. 2017; A. Savino et al. 2023; K. B. W. McQuinn et al. 2024c). However, this ignores the significant population of low-mass galaxies that have been disrupted by galactic tides by the present day and are now in the form of stellar streams and phase-mixed galaxies (I. B. Santistevan et al. 2020; N. Panithanpaisal et al. 2021; D. Horta et al. 2024). In A. Kundu et al. (2026), we used the FIRE-2 simulations to show that these disrupted substructures account for up to $\sim50\%$ of the total stellar mass budget of the proto-MW/LG and its environment at $z = 6$–9. Their inclusion improves the normalization of the reconstructed high-$z$ UVLF by factors of $\sim2$–3 and reduces the halo-to-halo scatter in the recovered faint-end slope by $\sim20$–40%$. These results indicate that disrupted substructures are crucial for a comprehensive near-far analysis, but the feasibility of incorporating them depends on whether they can be detected in sufficient numbers in current and upcoming surveys, and whether they preserve fossil-record information comparable to that of intact dwarf galaxies.

Disrupted substructures such as stellar streams and phase-mixed galaxies are difficult to detect because of their low surface brightness ($\mu\gtrsim27$ mag arcsec$^{-2}$; K. V. Johnston et al. 2008). However, recent wide-area surveys such as the Gaia mission (Gaia Collaboration et al. 2016, 2018, 2021), the Dark Energy Survey (DES) (B. Flaugher et al. (2015); T. M. C. Abbott et al. [1, 2]), the Southern Stellar Stream Spectroscopic Survey ($S^{5}$; T. S. Li et al. 2019), the Large sky Area Multi-Object fiber Spectroscopic Telescope (LAMOST; G. Zhao et al. 2012), and the Apache Point Observatory Galactic Evolution Experiment (APOGEE; S. R. Majewski et al. 2017) have discovered a substantial population of streams and amassed a large number of high-quality observations of their stars. The total number of known stellar streams from dwarf galaxies and globular clusters in the MW is now approaching $\sim100$ (K. Malhan et al. 2022; C. Mateu 2023). These detections have advanced our understanding of the MW’s assembly history (e.g., A. Helmi et al. 1999; R. Ibata et al. 2001; H. J. Newberg et al. 2002; S. R. Majewski et al. 2003; C. J. Grillmair 2006; H. J. Newberg et al. 2009; E. Balbinot et al. 2016; N. Shipp et al. 2018; K. Malhan et al. 2018; V. Belokurov et al. 2018; T. S. Li et al. 2019; L. Necib et al. 2020; K. Malhan et al. 2021; A. Helmi et al. 2018; M. Haywood et al. 2018; G. C. Myeong et al. 2018a,b; N. F. Martin et al. 2022; J. W. Johnson et al. 2023). However, the full accretion history of the MW, encoded in its stellar streams, is expected to extend beyond the sensitivity of these surveys, to surface brightnesses fainter than 30 mag arcsec$^{-2}$ (N. Shipp et al. 2023). With LSST (Ž. Ivezić et al. 2019), we will be able to probe this regime, enabling the detection of substantially fainter substructures within the MW halo and obtaining their resolved CMDs to sufficient depth for reconstruction. Recent predictions suggest that LSST could detect $\sim89$ MW satellite galaxies within its footprint (K. Tsiane et al. 2025), up from the $\sim65$ currently known (e.g., A. B. Pace 2025, and references therein), and could significantly increase the number of known stellar streams (N. Shipp et al. 2023). Realizing the potential of these new data for the near-far technique requires accurate predictions for what LSST will actually recover.

In A. Kundu et al. (2026), we presented a preliminary detectability estimate of our simulated intact galaxies and stellar streams, finding that $\sim50\%$ of streams passed our simple detectability test, which did not distinguish between streams detectable through their tidal tails and those detectable only through their bound progenitors. In this paper, we extend that work by constructing mock LSST catalogs of both intact satellites and stellar streams from the same simulations, and by estimating their detectability through a more realistic approach.

We focus on intact galaxies and stellar streams in this paper and defer the detectability of fully phase-mixed galaxies to future work, since the detection methods differ substantially once a galaxy’s stars are fully incorporated into the stellar halo of the host. The mock catalogs generated in this work are released publicly alongside this paper, which serves as documentation of the public data. While our analysis is focused on probing the faint-end slope of the UVLF at $z \gtrsim6$ through the near-far technique, the catalogs have a broad range of applications for studies of hierarchical structure formation, tidal disruption, orbital dynamics, and the nature of dark matter. We describe the FIRE-2 simulations and the classification of present-day simulated galaxies in § 2; in § 3 we walk through the construction of our mock catalogs. In § 4 we present the detectability analysis of our mock streams and intact galaxies and highlight their present-day properties in § 5. We compare the mock and simulated SFHs and analyze them in § 6; in § 7 we present our reconstructed UVLFs and the faint-end slope obtained using the detectable substructures.

Simulations

FIRE-2 Simulations

We use cosmological zoom-in baryonic simulations from the Feedback in Realistic Environments (FIRE) project (A. Wetzel et al. 2023, 2025)¹ to construct our catalog of present-day intact and disrupted galaxies. The simulations are run with the GIZMO code (P. F. Hopkins 2015; P. F. Hopkins et al. 2018a), which uses a meshless finite-mass (MFM) Lagrangian method for accurate hydrodynamics. Baryonic processes such as star formation and stellar feedback are implemented using the FIRE-2 numerical prescriptions (P. F. Hopkins et al. 2018a), with realistic radiative heating and cooling processes in the multiphase interstellar medium (ISM) across temperatures from 10 K to $10^{10}$ K. Star particles form out of gas parcels satisfying a Jeans instability criterion, and inherit the mass and elemental abundances from their parent gas cells. Each star particle represents a single-age, single-abundance stellar population, whose evolution is based on stellar population synthesis models from STARBURST99 (C. Leitherer et al. 1999) assuming a Kroupa initial mass function (IMF; P. Kroupa 2001). FIRE-2 models stellar feedback processes including supernovae, radiation pressure, continuous mass loss, photoionization and photoelectric heating in a time-resolved manner.

The simulations adopt a spatially uniform UV/X-ray ionizing background that evolves with redshift (C.-A. Faucher-Giguère et al. 2009). This model has two important deficiencies for our purposes: it neglects the inhomogeneous nature of reionization, which may affect the degree of halo-to-halo scatter (J. Borrow et al. 2023), and assumes an earlier reionization redshift ($z \sim10$) than recent measurements suggest ($z \sim7$; Planck Collaboration et al. 2020). However, all simulations used in this work consistently employ this model; hence these assumptions will not introduce any simulation-to-simulation difference in our analysis. Additionally, these simulations include a spurious heating of neutral gas at $\lesssim1000$ K at $z \gtrsim10$ due to cosmic rays, which suppresses star formation in low-mass halos at high redshifts. After reionization begins ($z \lesssim9$), however, this has a negligible effect (P. J. Gandhi et al. 2024), so we restrict our analysis to $z \lesssim9$. For a more detailed description of the physics implemented in FIRE-2 simulations, see P. F. Hopkins et al. (2018a,b).

The simulations successfully generate isolated MW-like galaxies and MW/M31-like galaxy pairs, reproducing many observed properties of the MW (X. Ma et al. 2017; P. F. Hopkins et al. 2018a; R. E. Sanderson et al. 2020; S. Ansar et al. 2025) and its satellite galaxies (A. R. Wetzel et al. 2016; S. Garrison-Kimmel et al. 2019a,b; J. Samuel et al. 2020, 2021). The central massive MW-like galaxy that forms self-consistently in these simulations efficiently destroys satellite galaxies and low-mass subhalos that orbit within the inner $\sim30$ kpc through gravitational tidal forces (S. Garrison-Kimmel et al. 2017), yielding a radial distribution of surviving satellite galaxies consistent with observations (J. Samuel et al. 2020). The tidal disruption of well-resolved simulated satellite galaxies also leads to a realistic population of stellar streams in terms of number and mass (N. Panithanpaisal et al. 2021; D. Horta et al. 2023; N. Shipp et al. 2023; A. Kundu et al. 2026), although they appear to be distributed at larger radii than the known MW streams (N. Shipp et al. 2023). The discrepancy in the radial distribution of streams can be partially attributed to resolution effects that can change whether disrupted galaxies in the inner halo are classified as streams or phase-mixed, but not whether they disrupt or not. A higher-resolution halo in FIRE have shown the radial distribution of streams in that halo to be broadly consistent with that of the MW streams (Y. J. Meziani et al. 2026). Furthermore, streams with still-bound components are likely to survive longer in higher-resolution simulations than lower-resolution ones (A. H. Riley et al. 2025; Y. J. Meziani et al. 2026), but they still have tidal tails and are not classified as satellites in the simulations. We will discuss the implications of this finding for our results in § 5.1. Streams from globular clusters do not naturally form in the FIRE-2 simulations due to limited resolution.

In A. Kundu et al. (2026), we constructed a catalog of surviving intact galaxies, stellar streams, and phase-mixed galaxies from 13 MW-mass halos in the core FIRE-2 simulations: seven isolated MW-mass halos from the Latte suite (A. R. Wetzel et al. 2016), and three LG-like MW/M31 pairs from the ELVIS on FIRE suite (S. Garrison-Kimmel et al. 2019a), along with their surrounding low-mass galaxies (field dwarfs). All simulations follow a $\Lambda$CDM model based on the cosmological parameters from the Planck Collaboration et al. (2020) results. The star/gas particles have an initial mass resolution of $\sim3500$–$7100\ M_\odot$ and a spatial resolution of 1–4 pc, while the dark matter (DM) particles have mass resolution of 19000–39000 $M_\odot$ and a spatial resolution of $\sim40$ pc. The main host halos have total masses of $M_{200\mathrm{m}} \approx1$–$2 \times10^{12}\ M_\odot${^6} and stellar masses of $M_\star\approx10^{10}$–$10^{11}\ M_\odot$, spanning the measured values for the MW (R. E. Sanderson et al. 2020). For each simulation, cosmological zoom-in initial conditions are generated at $z = 99$ from the MUSIC code (O. Hahn & T. Abel 2011), with the zoomed regions selected within periodic cosmological boxes of sizes 70–172 Mpc. In each simulation, 600 snapshots are saved from $z = 99$ to $z = 0$, with a typical spacing of $\sim25$ Myr after $z \sim2$, enabling accurate tracking of disrupted galaxies.

Classification of Present-day Galaxies in the Simulations

Using the halo catalogs based on the Rockstar 6D phase-space halo finder (P. S. Behroozi et al. 2013a) and the merger trees constructed using CONSISTENT-TREES (P. S. Behroozi et al. 2013b), we track each subhalo (and its star particles) that falls within the virial radius of the main host galaxy, starting from a high redshift ($z_{\mathrm{MR_{3:1}}}$) to the present day. This is the redshift at which the stellar-mass ratio of the most massive halo to the second-most massive halo is 3:1, which is when the proto-host galaxy emerges (I. B. Santistevan et al. 2020; D. Horta et al. 2024). We refer the interested reader to L. Necib et al. (2019); J. Samuel et al. (2020); P. J. Gandhi et al. (2024); A. Kundu et al. (2026) for discussions about the halo finding criteria and star particle assignments at both low- and high-$z$. For more details about the algorithm used to track the galaxies, we refer to N. Panithanpaisal et al. (2021); A. Kundu et al. (2026); L. S. Mandacarú Guerra et al. (2026). We classify the present-day galaxies into four categories:

1. Surviving satellite galaxies: Intact low-mass galaxies within the virial radius (radii)³ of the isolated (paired) MW-like (MW/M31-like) host(s) at $z = 0$.

2. Stellar streams: Disrupted low-mass galaxies which are torn and stretched along their orbits by tidal forces from the host galaxy.

3. Phase-mixed galaxies: Disrupted low-mass galaxies which were once stellar streams, but are now merged with the host galaxy and occupy its equilibrium phase-space distribution.

4. Field galaxies: Intact galaxies that are outside the virial radius (radii) of the isolated (paired) MW-like (MW/M31-like) host(s) at $z = 0$.

We use three criteria based on the physical properties of streams to identify present-day stellar streams within $R_{200\mathrm{m}}$ of the main host galaxy in each simulation. First, the number of star particles in the stream candidate must be larger than 100, which is equivalent to a stellar mass of $\sim5 \times10^{5}\ M_{\odot}$. This is a practical limit on the required resolution to apply the second and third criteria, which are based on the application of Liouville’s theorem, but is also a plausible lower limit on surface brightness for most streams. As we will discuss in § 5.1, this selects streams with much lower stellar mass than required to create detectable tidal tails above our surface brightness cutoff (34 mag arcsec$^{-2}$).

Stream stars are unbound from their progenitor and spread out in position as their local velocity dispersion decreases. Thus we require that the maximum pairwise separation between any two star particles in the stream candidate must be greater than 120 kpc, while its local velocity dispersion must be below a stellar-mass dependent threshold (see Equation 2 of N. Panithanpaisal et al. (2021)). If both the positions and velocities of the stars are widely distributed, we consider them to be phase-mixed, while if both are narrowly distributed then we conclude that the galaxy is still intact. To summarize, we classify a galaxy as (a) a stream if it satisfies all three criteria; (b) a phase-mixed galaxy if it satisfies the first two criteria; and (c) a satellite galaxy if it satisfies just the first criterion. However, we manually reduce the stellar mass requirement to $10^{4.5}\ M_{\odot}$ for the case of intact galaxies, since they are more compact in comparison to streams, and therefore are easier to detect. Additionally, for all our classified substructures within the virial radius, we exclude star particles associated with the present-day host galaxies themselves. We also exclude star particles associated with the present-day host galaxies themselves. For more details, we refer the reader to N. Panithanpaisal et al. (2021); A. Kundu et al. (2026).

For the field galaxies, we select all the luminous halos with $M_{\star} \ge10^{4.5}\ M_{\odot}$ outside $R_{200\mathrm{m}}$ of the main host galaxy within a region uncontaminated by low-resolution DM particles. This differs from the selection used in A. Kundu et al. (2026), where we imposed a distance cut of 2 Mpc from the center of each isolated MW-like galaxy, or from the geometric center of mass of the LG-like pairs. These criteria ensure that our sample does not suffer from numerical issues due to very low-mass halos and safely excludes any contaminated halos, while extending our field galaxy sample to larger distances. However, this leads to the field galaxies being drawn from a common distribution in case of the paired hosts, the consequences of which we discuss in § 7. Table 1 shows the number of intact and tidally-disrupted galaxies in each of our 13 host galaxies in the simulations.

Simulation Name	$M_{200\mathrm{m}}$	$M_{\star,90}$	$R_{200\mathrm{m}}$	$z_{\mathrm{MR3:1}}$	Satellite Galaxies	Stellar Streams	Phase-mixed Galaxies	Field Galaxies
m12i$^{A}$	1.2	6.3	336	3.18	20 (5)	13	17	7
m12b$^{B}$	1.4	8.5	358	1.59	19 (3)	13	10	18
m12f$^{C}$	1.7	7.9	380	3.77	28 (4)	12	29	26
m12c$^{B}$	1.3	5.8	351	1.32	33 (9)	11	7	55
m12m$^{D}$	1.6	11.0	371	1.68	31 (10)	13	6	40
m12r$^{E}$	1.1	1.7	321	1.02	24 (5)	5	6	20
m12w$^{E}$	1.1	5.7	319	2.92	31 (12)	5	14	21
Romeo$^{B}$	1.3	6.6	341	5.70	23 (5)	20	37	81
Juliet $^{B}$	1.1	3.8	321	4.72	22 (10)	15	29
Romulus$^{F}$	2.1	9.1	406	1.70	38 (7)	14	25	58
Remus$^{F}$	1.2	4.6	339	0.98	20 (3)	15	6
Thelma$^{F}$	1.4	7.1	358	4.29	20 (8)	11	24	62
Louise$^{F}$	1.1	2.6	333	3.32	27 (7)	12	24

Creating Mock LSST Catalogs

We use py-ananke⁴ (A. Thob et al. (2024), a Python package based on the ananke framework (R. E. Sanderson et al. (2020)), to generate mock astrometric and photometric observations of synthetic stars in the satellites, streams, phase-mixed galaxies, and field galaxies identified in § 2.2, as observed by LSST from within the host galaxy. We create a total of 13 mock catalogs: 7 from isolated MW halos and 6 from the three LG-scale simulations, each of which contains two MW-like halos. Here we provide a working summary of the steps employed to run the pipeline in this work and refer the interested reader to R. E. Sanderson et al. (2020); A. Thob et al. (2024) for more details.

Preparing the simulation data

For each host galaxy, we establish a galactocentric coordinate system where the observer is placed on a solar viewpoint, specifically, a phase-space position for the local standard of rest (LSR). First, we determine the galactocentric coordinate frame, centered on and aligned with the stars in the galaxy such that the galaxy disk lies in the $X$-$Y$ plane and the disk angular momentum lies Note—Simulation Name: References listed below; the top seven are isolated MW-like and bottom six (three pairs) are paired MW/M31-like. $M_{200\mathrm{m}}$ and $R_{200\mathrm{m}}$: total mass and spherical radius within which the mean density is 200 times the matter density of the Universe. $M_{\star,90}$: stellar mass within a spherical radius that encloses 90% of the stellar mass within 20 kpc. $z_{\mathrm{MR3:1}}$: the redshift at which the stellar-mass ratio of the most massive halo to the second-most massive subhalo is 3:1, which is when the dominant host galaxy emerges. Numbers of surviving satellite galaxies ($M_{\star} \gtrsim10^{4.5}M_{\odot}$; in parentheses, the number above the stream and phase-mixed mass cut ($M_{\star} \gtrsim5 \times10^{5}M_{\odot}$)), stellar streams ($M_{\star} \gtrsim5 \times10^{5}M_{\odot}$), phase-mixed galaxies ($M_{\star} \gtrsim5 \times10^{5}M_{\odot}$), and field galaxies ($M_{\star} \gtrsim10^{4.5}M_{\odot}$) are determined as described in § 2.2.

References—(A) A. R. Wetzel et al. (2016), (B) S. Garrison-Kimmel et al. (2019a), (C) S. Garrison-Kimmel et al. (2017), (D) P. F. Hopkins et al. (2018a), (E) J. Samuel et al. (2020), (F) S. Garrison-Kimmel et al. (2019b).

along the $Z$ axis. As discussed in R. E. Sanderson et al. (2020), we use the young stars (age $<1$ Gyr) to define the disk angular momentum direction, since this selection tracks the plane of the thin gaseous disk but is more stable over time. The center position of the host galaxy is determined by applying the shrinking-spheres method to the star particles iteratively and the center-of-mass velocity is computed using the star particles within 15 kpc of the center of the host DM halo (identified by the halo catalog) (R. E. Sanderson et al. (2020)).

In this frame we assume that the Sun is located at $R_{\odot} = 8.2$ kpc (J. Bland-Hawthorn & O. Gerhard 2016). Then, we compute the moment of inertia tensor using young star particles (age $<1$ Gyr) inside the Solar radius. This establishes the $Z$ direction perpendicular to the plane of the galactic disk. As the norm for standard Galactocentric coordinates, we choose the orientation of the total angular momentum of the galactic disk to point in the $-Z$ direction, so that the host galaxy rotates clockwise.

If the LSR position is on the $-X$ axis, the rotation of the galaxy carries it to the $+Y$ direction, consistent with the standard conventions used for Galactic coordinate systems.

After establishing a galactocentric coordinate frame, we define a solar viewpoint and LSR in it. We initially start with a Solar position at $(-8.2, 0, 0)$ kpc and select all the star particles within 200 pc of this position. This allows us to define a sufficiently “local” LSR, with enough particles to avoid sampling noise. We then compute the median position and velocity using all the star particles within this “local Solar neighborhood” and use this to define our LSR position and velocity. This phase-space coordinate serves as our 6D observer position for LSST.

Generating the density field

py-ananke uses the local density of the selected star particles in phase space to determine the volume over which to distribute the mock stars sampled from that particle. We apply EnLink-2.6 (S. Sharma & K. V. Johnston 2009) to the positions, velocities, and masses of the selected star particles to compute the phase-space density at the location of each simulated star particle⁵. The six-dimensional densities are estimated using the adaptive, entropy-based density estimation scheme described in S. Sharma & K. V. Johnston (2009) with 16 nearest neighbors, and smoothed using bootstrapping. The estimated densities are converted into equivalent single-dimensional, “cubic” densities in the position and velocity subspaces and then inverted to return a local estimate for the bandwidth of a spherical Epanechnikov kernel distribution in position and velocity, such that sampling this distribution reconstructs the input density distribution. The resulting kernel lengths are written to a file which we use as inputs to py-ananke. Because the density estimate uses 16 nearest neighbors, we restrict the sample to substructures with more than 20 star particles. This only removes some of the intact low-mass galaxies with stellar masses between $10^{4.5}\ M_\odot$ and $10^5\ M_\odot$ and not the streams, phase-mixed galaxies and higher-mass intact galaxies. We classify these low-mass galaxies as undetectable, but keep them in the simulated population.

Creating a mock LSST catalog using py-ananke

The galactocentric positions and velocities and the list of kernel lengths for the selected star particles are passed to py-ananke along with their masses, ages, and iron abundances [Fe/H], which are used to determine the number of stars spawned from the single-age, single-metallicity stellar population represented by each particle. We also pass additional information from the parent star particle on to each synthetic star: chemical abundances for He, C, N, O, Ne, Mg, Si, S and Ca, $\alpha$ abundances (mean abundance using O, Mg, Si, Ca), and the distance of the star particle from the host-galaxy center at formation (also known as the birth radius).

py-ananke samples stellar masses from a Kroupa IMF (P. Kroupa 2001) and calculates their absolute magnitudes in the LSST filters using PARSEC isochrone models (A. Bressan et al. 2012) matched to the metallicity and age of each star particle. Specifically, we use a grid of isochrones for non-rotating single stars generated using PARSEC version 1.2S (Y. Chen et al. 2014; J. Tang et al. 2014; Y. Chen et al. 2015), incorporating COLIBRI models of asymptotic giant branch evolution (P. Marigo et al. 2013; P. Rosenfield et al. 2016; G. Pastorelli et al. 2019, 2020) with $\eta_{\mathrm{inTPC}} = 10$ (P. Marigo et al. 2017) and $\eta_{\mathrm{Reimers}} = 0.2$. Circumstellar dust is modeled using scaling relations from P. Marigo et al. (2008) for the recommended multi-component dust compositions: 60% silicate/40% AlOx for M stars and 85% AMC/15% SiC for C stars (M. A. T. Groenewegen 2006). The isochrone grid spans 71 age values evenly spaced from $\log_{10}(\mathrm{age}/\mathrm{yr}) = 6.6$ (3.98 Myr) to $\log_{10}(\mathrm{age}/\mathrm{yr}) = 10.1$ (12.6 Gyr) and 34 metallicities from $Z = 0.0001$ to $Z = 0.03$ (approximately equivalent to $-2.3 < [\mathrm{M}/\mathrm{H}] < 0.2$). The age and metallicity of each star particle are mapped to the closest isochrone in the grid for sampling; in other words, we do not interpolate between isochrones. For the few star particles that fall outside of the range spanned by the model isochrones, we use the nearest edge of the grid.

We also define our survey parameters for py-ananke, based on projections for LSST. We set the observer parameter to the location (in kpc) of the observer in galactocentric coordinates as calculated in § 3.1, which roughly corresponds to $(-8.2, 0, 0)$ kpc. We set fsample=1 for a full-sized survey. We choose rshell=4000 kpc to sample stars from star particles in the tidal arms of the streams that extend beyond the host halo, as well as the field galaxies. We calculate magnitudes of the sampled stars in the six bands of the LSST filter system (photo_sys = ‘LSST’) and specify the limiting magnitude of the survey in one LSST filter with the cmd_magnames and cmd_box parameters. More details about these parameters are available at — https://py-ananke.readthedocs.io/en/main/

Since LSST has different limiting magnitudes in its different filters, which stars are selected after sampling from the IMF depends on the filter and magnitude limit we choose. We choose to accept stars into the survey that are detected in at least one of the six bands; since the red filters are more efficient, we set the limiting magnitude in the $z$ band to 29 for the preliminary survey, before making the final cut. This is somewhat fainter than the published 10-year co-added apparent magnitude limit for LSST in this band (25.6; F. B. Bianco et al. 2022) in order to allow for the spread in distances of the stars spawned by each star particle and the stochastic rounding used during the sampling of the IMF to generate an integer number of stars from each particle. We achieve significant speedup by estimating the lower limit in mass that corresponds to the faint limit of the survey at the location of each star particle and limiting the IMF sampling to this range, but for galaxies with very few simulated star particles this simplification can combine with the distance spread and stochastic sampling to erroneously produce zero stars above the strict limiting magnitude of the survey. Choosing an apparent magnitude which is a bit fainter than our survey limit ensures that we get some synthetic stars for these very low-mass galaxies. Before analyzing our mock catalogs, we then re-select stars with apparent magnitudes brighter than the LSST faint limit in at least one band. Hence, to generate a mock catalog using py-ananke it is advisable to use the most sensitive filter for defining the magnitude limit, and to choose a preliminary apparent magnitude limit a bit fainter than the desired final faint limit in that filter.

Running py-ananke using the above methods generates mock LSST catalogs for the substructures identified in each of the 13 simulated MW analogs in Table 1, assuming perfect observation of their stellar populations without any foreground stars. For a few very low-mass galaxies, we get zero synthetic stars based on our survey parameters; we assume these to be undetectable. The data contained in the mock catalogs is described in detail in Appendix § 13. In brief, we generate sky positions, parallaxes, proper motions, and radial velocities; absolute and apparent magnitudes in the six LSST bands; stellar parameters and chemical abundances; and indices mapping each star to its parent star particle in the simulation.

To standardize file size and for greater accessibility, the resulting catalogs are saved in a series of files corresponding to a multiscale Healpix tessellation of the sky in Galactic coordinates. The tessellation follows the Multi-Order Coverage (MOC) mapping convention defined by P. Fernique et al. (2022) where a Healpix NUNIQ indexing scheme allows for each Healpix cell to be uniquely identified across all resolution levels. The determination of the tessellation is done iteratively through the resolution levels using the distribution of star particles from the simulation and their initial masses that were used as input of py-ananke. Starting at level 0 with the full set of particles, the process checks all Healpix cells of that level (12 of them at level 0) and compares for each of them the total mass from the particles they contain to a threshold mass of $10^{8}\,M_{\odot}$: if the total is below the threshold, the cell’s NUNIQ id is saved and the corresponding particles are removed from the set of particles to consider in the next iteration. Once all cells of the level have been considered, if particles still remain in the set, the process iterates the same procedure in the next resolution level with the remaining particles, and continues to do so in the deeper levels until no particles are left. The result is a set of NUNIQ ids corresponding to the full multiscale Healpix tessellation that exhaustively covers the entire sky. While this process does not allow for the consideration of particles kernels overlapping over multiple cells where their center are sometimes not located — considering this would require a more involved algorithm — it provides a good first approximation of a multiscale tessellation that matches the expected distribution of synthetic stars for accessibility purposes.

The mock catalogs generated in this work contain perfect information for all the synthetic stars that would be observed by LSST, but do not account for either galactic extinction and reddening or LSST’s photometric errors. We do not expect extinction and reddening to have a large effect on our results since most of the streams are well above or below the galactic disk. In the future, we plan to use LSST’s StreamObs⁶ software to model LSST uncertainties. Although our mock catalogs adopt the 10-year co-added magnitude cuts for completeness, we use the 1-year magnitude cuts for making predictions in this work since the proper motions of the stars will likely play a significant role in assigning them to streams.

Data availability

All our LSST synthetic observations of satellite galaxies, stellar streams, phase-mixed galaxies, and field galaxies, are available via Globus.⁷ Anyone can access the data after logging in with a free Globus account. We include the substructures from all the 7 isolated MW-like hosts (simulation names listed as m12i, m12b, m12f, m12c, m12m, m12r, and m12w), and the 3 MW/LG-like hosts (simulation names listed as RJ, RR, TL). Within each host directory we have the synthetic observation of a classified substructure, with its own subdirectory in the format: <simname>_<cflag>_<ID>_Rubin_lsr_axes. In case of the isolated suite, simname corresponds to the simulation name, e.g., m12i. In case of the paired suite it also corresponds to the simulation name but with an additional flag of the host, e.g., RJ_host1 or RJ_host2, identifying the respective host. cflag can be either cat or hal, where cat includes the stellar streams, phase-mixed galaxies, and satellite galaxies with $M_{\star} \gtrsim5 \times10^{5}\,M_{\odot}$, and hal includes the field galaxies and the lower-mass satellites ($10^{5}\,M_{\odot} \lesssim M_{\star} \lesssim5 \times10^{5}\,M_{\odot}$). We include the ID of each identified substructure in the <simname>_object_ids.pkl in each simulation directory. The key names stream, sat, and phase_mixed belong to the category cat, while other_dwarf_indices and field_gal_indices belong to the hal category and corresponds to the halo indices from the halo catalogs. The key names very_low_mass_other_dwarf_indices and very_low_mass_field_dwarf_indices includes the low-mass halos with $M_{\star} \gtrsim10^{4.5}\,M_{\odot}$ but have less than 20 star particles. These halos do not have any synthetic observation due to our selection criterion of star particles. Each synthetic observation directory includes: the hdf5 files survey.<simname>_<cflag>_<ID>_Rubin_ls r_axes.N.h5, where $N$ is the Healpix partition ID; the survey parameter file survey_params; and the density file in ebf format.

Detectability of Streams and Intact Galaxies

We follow an approach similar to that of N. Shipp et al. (2023) to estimate the detectability of the mock streams and satellites, with a few important differences. N. Shipp et al. (2023) estimated the detectability of the streams and satellites for DES using the respective substructure classifications from N. Panithanpaisal et al. (2021). In this work, we use the classified substructures from A. Kundu et al. (2026), which contain some additional objects which are massive ($M_{\star} \gtrsim10^{9}\,M_{\odot}$) and morphologically complex, and also the low-mass satellites ($10^{4.5}\,M_{\odot} \lesssim M_{\star} \lesssim10^{6}\,M_{\odot}$) — both of which are excluded from N. Panithanpaisal et al. (2021). Additionally, we also estimate the detectability of field galaxies based on the criteria discussed in B. Mutlu-Pakdil et al. (2021) for LSST-based searches of low-mass galaxies beyond the MW and LG. We do not estimate the detectability of the phase-mixed galaxies and defer that to future work.

Stream Detectability

As described in N. Shipp et al. (2023), we determine the detectability of the streams following a method motivated by matched-filter stream searches in the DES data (N. Shipp et al. 2018). The detectability of the streams can be estimated using the average surface brightness from the observed stars within the survey magnitude limits. Streams that have low luminosity and are more diffuse have fainter surface brightness, making it more difficult to disentangle them from background and foreground stellar populations. This method also takes into account that more distant streams will have fewer detectable stars, and hence will have a fainter observed surface brightness. Instead of using a converted $V$-band magnitude to compute the average surface brightness of the streams as N. Shipp et al. (2023) did for DES, we directly use the LSST $g$-band, which is close to the $V$-band. Since it is difficult to estimate the surface brightness limit for resolved stellar populations in LSST before the availability of real survey data, we keep the surface brightness limit as $\mu_g \leq34\,\mathrm{mag\,arcsec}^{-2}$, the same as the surface brightness limit in $V$-band used in N. Shipp et al. (2023). We restrict our results to the magnitude limits for 1-year exposure and hence obtain conservative estimates with respect to the survey depth. Data from additional years and improved detection methods will reach deeper magnitudes and fainter surface brightnesses. Users can explore this avenue, as our mock catalogs go deeper in magnitudes than the 1-year limits.

In order to determine the average surface brightness of the tidal tails of a mock FIRE stream, we start by masking all the stars within $4r_{1/2}$ of the bound progenitor (if any), using the $r_{1/2}$ parameter derived in § 4.2. Then we follow the exact approach of N. Shipp et al. (2023) and convert to a coordinate system in which the stream lies approximately along the equatorial plane, where $\phi_1$ is the coordinate along the stream and $\phi_2$ is the coordinate perpendicular to the stream. The stars belonging to a long, convoluted stream may not share the exact same orbital plane. Hence, we estimate the best-fit RA and DEC of the orbital pole $(\alpha,\delta)$ for each stream by minimizing the quantity $D(\alpha,\delta) = \sum_i \left| \frac{\pi}{2} - D_i(\alpha,\delta) \right|$, where $D_i(\alpha,\delta)$ is the angular distance from the pole to each stream member star.

To capture the variation in surface brightness along the stream, we divide it into segments of equal length (15 deg) and compute the detectability of each one individually. Then, to compute the $\phi_2$ position of each endpoint, we take stars that are within 3 deg of $\phi_1^{\mathrm{end}}$. If this region contains fewer than 50 stars, we exclude this endpoint as it is likely empty, and automatically flag its neighboring segments as undetectable. For sufficiently populated regions, we identify the corresponding $\phi_2^{\mathrm{end}}$ to be the peak in the $\phi_2$ distribution of the selected stars. After determining the detectability of each of these segments individually, we classify a stream as detectable provided it has at least one detectable segment (after masking the bound progenitor). Any system with more than one detectable segment is considered to be a single detectable stream — consistent with how disconnected but physically associated segments of known Milky Way dwarf galaxy streams are now treated as one single stream. Examples include ATLAS-Aliqa Uma ((S. E. Koposov et al. 2014; N. Shipp et al. 2018; T. S. Li et al. 2021)), Orphan-Chenab ((C. J. Grillmair 2006; V. Belokurov et al. 2007; N. Shipp et al. 2018; S. E. Koposov et al. 2019)), and Sagittarius ((R. A. Ibata et al. 1994; S. R. Majewski et al. 2003; S. Hasselquist et al. 2019; B. D. Johnson et al. 2020)).

Within each segment, we apply another coordinate transformation such that the great circle joining the two endpoints aligns with the equatorial plane. Then, we measure the segment width in $\phi_2$ by fitting the peak interval — defined as the smallest $\pm1\sigma$ interval enclosing the peak of the distribution. For our analysis, we restrict to only the inner $\pm1\sigma$ range to avoid bias in the detectability estimate from the diffuse components at large $\phi_2$.

Unlike N. Shipp et al. (2023), who converted DES $g$- and $r$-band magnitudes to the $V$ band, we estimate the total luminosity directly in the LSST $g$-band (absolute magnitude $M_g$). Using the total luminosity and area of each segment, we calculate the average surface brightness ($\mu_g$). We use a detectability limit of $\mu_g \le34$ mag arcsec$^{-2}$. In Figure 1, we show some examples of mock streams projected along the stream coordinates. The top left panel shows a stream with both a detectable bound progenitor and detectable tidal tails and the top right panel shows a stream with no bound remnant but with most of the tidal tails being detectable. The bottom left panel shows a stream which is entirely detectable, and the bottom right panel shows a stream which only has a detectable bound part without any tidal tails — meaning it will be detectable as a satellite.

We emphasize that our detectability results may change depending on the actual surface brightness limits LSST will achieve. Additionally, we do not consider background galaxies to estimate detectability. This requires star-galaxy separation, particularly at faint magnitudes, where unresolved galaxies can contaminate the foreground stellar populations. Hence, our estimates are optimistic by neglecting background contamination, and users of our catalog are encouraged to explore to what extent the results can change when unresolved background galaxies are included.

Intact Galaxy Detectability

We also estimate the detectability of intact galaxies and stream progenitors in our catalog at $z = 0$. This category includes the satellites identified by our pipeline ($M_\star\gtrsim5 \times10^5\ M_\odot$; § 2.2), the additional low-mass satellites ($10^{4.5}\ M_\odot\lesssim M_\star\lesssim5 \times10^5\ M_\odot$) and the field galaxies identified from the Rockstar halo catalogs. Some of the streams at $z = 0$ still have bound parts, which are identified as bound halos in the halo catalogs. We call them “stream progenitors” and also estimate their detectability. If the stream has no detectable tidal tail using the method described in § 4.1, but has a detectable bound part, we classify it as a satellite in our subsequent analysis.

Using the mock LSST observations of these objects, we follow a similar approach to K. Bechtol et al. (2015); A. Drlica-Wagner et al. (2015, 2020); N. Shipp et al. (2023) and fit an ellipsoidal Plummer (1911) profile to measure the half-light radius ($r_{1/2}$). However, unlike N. Shipp et al. (2023), we do not use the code ugali⁸ to compute the parameters. We instead fit an ellipsoidal Plummer profile to the spatial distribution of detected stars using Nelder-Mead optimization (F. Gao & L. Han 2012) to get $r_{1/2}$ in the $V$-band, for which we convert the $g$- and $i$-band LSST magnitudes to visual magnitudes using the relations in M. N. Porter et al. (2025). Furthermore, we calculate the total luminosity ($M_V$) by taking the median distance of the mock member stars as the distance to the intact galaxy. Unlike N. Shipp et al. (2023), we do not correct for the contribution to the luminosity from unobserved stars (due to LSST magnitude limits) using the isochrones in ugali. We will show in § 5.2 that despite these approximations and simplifications, our results still agree well with the observed properties of MW satellites above our simulation resolution limits. However, we encourage the users of our mock catalogs to try more sophisticated methods for calculating these properties.

The detectability of satellites can be estimated using an analytical approximation as a function of absolute magnitude, half-light radius, and distance. This was presented by A. Drlica-Wagner et al. (2020) for DES and then updated for LSST in K. Tsiane et al. (2025). We use the derived parameters above to estimate the detectability of satellites and stream progenitors within the virial radius of the host galaxy. Combining these results with those for detectable streams (§ 4.1), we can divide our simulated streams into streams with detectable tidal tails, detectable progenitors only (streams mistakenly classified as satellites), or no detectable component.

For the field galaxies we follow a different approach. Using our estimates of the absolute magnitude and half-light radius, we compare with the size-luminosity limits presented in B. Mutlu-Pakdil et al. (2021) for dwarf galaxies out to $\sim3.5\ \mathrm{Mpc}$ using LSST. We use their 50% detectability regions in $r_{1/2}$–$M_V$ in the two closest distance bins ($< 1.5\ \mathrm{Mpc}$, $< 3.5\ \mathrm{Mpc}$) and flag the galaxies within these regions as detectable. We highlight that their detectability estimates are based on 10 years of LSST data, while our properties are calculated using 1-year magnitude cuts. This means our detectability estimates of field galaxies are conservative and will only improve if using the 10-year magnitudes.

Present-day Properties of Streams and Intact Galaxies

Detectability Estimates of Streams and Intact Galaxies

Using our detectability estimates for the stream populations and intact galaxies, we compute the number of detectable systems. In Figure 2, we show the number of detectable streams in each of the 13 host galaxies. The bottom segment of each bar shows the number of streams with detectable tidal tails, the middle segment represents the number of streams with only detectable progenitors (i.e., would be detected as satellites), and the top segment shows the number of undetectable streams (i.e., neither the tails nor the bound progenitor is detectable). The total height of each bar shows the total number of streams with $M_\star\gtrsim5 \times10^5\ M_\odot$ identified in each host galaxy. The median number of streams with detectable tidal tails is 3, which is consistent with the predictions of N. Shipp et al. (2023). Across the 13 simulations, we find an average of $\approx29\%$ of the streams to be detectable with tidal tails. We find that many of the streams in this stellar mass range and at this luminosity would remain undetectable even with LSST (within our chosen surface brightness limit). However, we would also mistake many of the streams for intact satellites, since their tidal tails are too diffuse to detect. We find on average $\approx38\%$ of the streams across the 13 hosts to be detectable as intact satellites. This can have important implications when reconstructing the high-$z$ UVLFs using the near-far technique. If the bound component contains only a small fraction of the total stellar mass of the stream, the inferred UV magnitude of its progenitor at high-$z$ can be biased unless the tidal tails are properly accounted for. We briefly explore the impact on the SFHs in § 6, but defer a detailed analysis to future work. We expect the percentages of streams with detectable tidal tails and streams detectable as satellites to change when using a higher-resolution simulation (within the same stellar-mass limit) as some of the systems classified into one category in the lower-resolution simulation can be re-classified into a different category in the higher-resolution simulation. For example, Y. J. Meziani et al. (2026) find two more stellar streams above the same mass resolution in the higher-resolution m12i simulation than the lower-resolution simulation used here, which would change our detectability estimates somewhat. These two streams are likely labeled as phase-mixed objects by our pipeline in the lower-resolution simulations we study here. On the other hand, we expect the resolution to have minimal effect on gravitationally bound objects (satellite galaxies and the parts of stream progenitors that remain bound). Thus, while the number of detected streams we predict may be a lower limit, the overall contribution of partially or fully disrupted galaxies remains significant, even if a few of them are classified differently at higher resolution.

We next consider the intact galaxies and the bound components of streams. Across the 13 hosts, an average of $\approx59\%$ of the satellites with $M_\star\gtrsim10^{4.5}\,M_\odot$ are detectable. The undetectable satellites are predominantly low-mass systems close to our resolution limit, including those with fewer than 20 star particles for which we do not generate mock stars (§ 3.2). Increasing the stellar mass cutoff to $M_\star\gtrsim10^5\,M_\odot$ improves the detectability fraction to $\approx96\%$, indicating that the undetectable satellites are almost entirely confined to $10^{4.5}\,M_\odot\lesssim M_\star\lesssim10^5\,M_\odot$ and are therefore limited primarily by simulation resolution rather than by LSST sensitivity. For the streams that have a surviving bound progenitor at $z = 0$, $\approx89\%$ of the progenitors with $M_\star\gtrsim10^{4.5}\,M_\odot$ are detectable. The fraction increases to $\approx95\%$ when we consider the progenitors with $M_\star\gtrsim10^5\,M_\odot$.

For the field galaxies, we average over the 7 isolated halos and the 3 LG-like pairs. In case of the paired suite, we consider the field galaxies to be outside the virial radii of the two main hosts, but use the properties as estimated from the more massive host to avoid overcounting. This is similar to the observed local field dwarfs — observing the field galaxies of both MW and M31 from MW. We find that $\approx42\%$ of the field galaxies with $M_\star\gtrsim10^{4.5}\,M_\odot$ out to $\sim3.5\,\mathrm{Mpc}$ are detectable, decreasing from $\approx47\%$ within 1.5 Mpc to $\approx28\%$ between 1.5 and 3.5 Mpc. For $M_\star\gtrsim10^5\,M_\odot$, these fractions increase to $\approx64\%$, $\approx71\%$, and $\approx46\%$, respectively. These fractions are summarized in Table 2.

Substructures	$M_\star\gtrsim10^{4.5} M_\odot$	$M_\star\gtrsim10^5 M_\odot$
Satellites	59%	96%
Bound stream progenitors	89%	95%
Field galaxies ($< 1.5$ Mpc)	47%	71%
Field galaxies (1.5–3.5 Mpc)	28%	46%
Field galaxies ($< 3.5$ Mpc)	42%	64%
Streams ($M_\star\gtrsim5 \times10^5 M_\odot$)
With detectable tidal tails	29%
Detectable only as satellites	38%

In Figure 3, we show the cumulative number of detectable galaxies as a function of distance from the observer and compare them against that of the MW satellites (left panel) and local field ($< 3\,\mathrm{Mpc}$) galaxies (right panel) taken from the Local Volume Database (LVDB)⁹ ((A. B. Pace 2025)). The lines show the median number of satellites (orange), stream progenitors (blue), satellites + stream progenitors (purple), and field galaxies (green), with the shaded regions showing the 68% halo-to-halo scatter across the 13 host galaxies (10 for field galaxies). The solid lines and the darker shaded regions represent the detectable population, while the dashed lines and the lighter shaded regions represent the simulated population. We also show the observed distribution of MW satellites and the local field dwarfs (black) with $M_\star\ge10^{4.5} M_\odot$ taken from LVDB. In the left panel, the observed MW satellite distribution is consistently above the median distribution of the detectable population and is within the 68% halo-to-halo scatter only beyond $\sim220$ kpc.

However, it agrees well with the distribution of the combined population of detectable satellites and stream progenitors (streams detectable as satellites) within the halo-to-halo scatter beyond $\sim150$ kpc, indicating that some of the observed MW satellites could be bound remnants of disrupted systems that would still appear as satellites even with LSST. The small discrepancy within $\sim150$ kpc is mainly due to our sample selection for generating the mock catalog. As discussed in § 3.2, we generate the mock catalogs for substructures with more than 20 star particles. Thus, the simulated low-mass intact galaxies and the bound stream progenitors with stellar masses within $10^{4.5}\,M_{\odot} \lesssim M_{\star} \lesssim10^{5}\,M_{\odot}$, but with fewer than 20 particles, are undetectable. This biases the distribution from the detectable population low, while the simulated distribution agrees well with the observed MW satellite distribution within 150 kpc. We further confirm this by selecting the low-mass intact galaxies, stream progenitors, and the MW satellites from LVDB with $M_{\star} \geq10^{5}\,M_{\odot}$. Figure A1 shows that the observed distribution agrees well within the 68% halo-to-halo scatter of the combined population of detectable satellites and stream progenitors over the entire distance range. In this case, the simulated distribution is very similar to the detectable distribution across all three populations.

We find these results to be consistent with simulated radial distributions in J. Samuel et al. (2020), but also note some discrepancies (see their Figure 3). We note that our median curve is below the observed MW satellite distribution within 150 kpc, while they have a better agreement. This is because our observed MW satellite population from LVDB with $M_{\star} \geq10^{5}\,M_{\odot}$ is not the same as theirs, since LVDB includes more recently discovered satellites which is missing in their analysis. Furthermore, their median radial distributions are combined across the last 1.3 Gyr in the simulations, while we only show them at the present day. Additionally, the underprediction from the detectable/simulated population within 150 kpc could be physical. The presence of the LMC in the MW makes its radial distribution more concentrated within 150 kpc at the present day than the typical simulated profiles (E. Patel et al. (2024)). Since not all of the FIRE host galaxies have an LMC-like counterpart, the median curve underpredicts the observed distribution. A few host galaxies which have an LMC-like satellite agree well with the observed distribution (see the thin lines in Figure A1).

In the right panel of Figure 3, the radial distribution of detectable field galaxies underpredicts the observed local field distribution beyond $\sim800$ kpc. The simulated distribution agrees well, indicating most of our field galaxies at these distances would go undetected. Our zoom-in simulations, particularly the isolated suite, do not capture the entire local field population within the resolution limit. Since zoom-in simulations are run at a higher-resolution after selecting a small region from a low-resolution cosmological volume, the field galaxies beyond a certain distance are unresolved. Additionally,

because the zoomed-in volume covers only a fraction of the cosmological volume, it fails to sample the entire field population. Although not shown here, the median from just the detectable field galaxies in the paired suite agrees well out to $\sim1.5$ Mpc, beyond which they too fail to capture the observed distribution. As discussed previously, we use the properties of the field galaxies as estimated from the more massive host in the paired suite to avoid overcounting.

We show the present-day cumulative stellar-mass functions (SMFs) of the streams and satellites with estimated detectability constraints in Figure 4. In the top row, we show the SMFs constructed using only the detectable substructures and in the bottom row, we show the total SMFs obtained from the simulations. The thin lines show the SMFs from each halo, while the thick lines and the shaded regions show their median and 68% halo-to-halo scatter. The blue and orange lines are from the isolated and paired suites, respectively. We compare the SMFs of the simulated and detectable streams in the left column. We find that most of the low-mass streams ($M_{\star} \lesssim10^{7.5}\,M_{\odot}$) are undetectable. In the bottom row of the middle column, we show the SMFs constructed using only the substructures identified as satellites in the simulation. But in the top row, the detectable substructures also include the streams without detectable tidal tails (only the bound progenitor is detectable). We find that the simulated and detectable SMFs are similar, which means the undetectable intact satellites are largely compensated for by the streams detectable as satellites. In the right column, we show the combined SMFs of the streams and satellites. We find that the simulated and detectable SMFs of the combined population are similar to each other. While constructing the SMFs for the detectable substructures, we use their total stellar masses from the simulation by assuming perfect reconstruction of the SFH from only the stars in the detectable regions. Hence, our estimates of the present-day SMFs show the ideal scenario where we are able to extract the total present-day stellar mass of each detectable substructure. In future work, we will explore the feasibility of this assumption by modeling the CMDs to reconstruct the SFHs.

Size-Luminosity Relation

We compare our estimated luminosity and size for the intact galaxies and bound stream progenitors with the observed dwarf galaxies in the MW and LG. We use the observed V-band magnitude ($M_V$) and half-light radius of the MW satellites and local field ($< 3$ Mpc) galaxies from LVDB. Figure 5 shows the $r_{1/2}$–$M_V$ relation for our mock satellites (orange), stream progenitors (blue), and observed satellites of the MW (green). The filled markers show the detectable galaxies, while the open ones show the undetectable galaxies according to our satellite detectability scheme discussed in § 4.2. For illustrative purposes, we also show the 50% detectability curves at distances of 50 kpc, 100 kpc, and 400 kpc, using the analytic relation from K. Tsiane et al. (2025) for LSST.

We see good agreement between the mock and observed size-luminosity relation for brighter ($M_V \lesssim-6$) and larger ($\log_{10}(r_{1/2}/\mathrm{pc}) \gtrsim3$) galaxies. This suggests our ellipsoidal Plummer profile fits are a good approximation in these regions. The agreement with the observed size-luminosity indicate similar conclusion as in C. Klein et al. (2024), showing that the half-light radii from synthetic observations are better estimates for the observed sizes than particle-based half-mass radii or half-light radii assuming a constant mass-to-light ratio. Hence, even though FIRE-2 galaxies appear “puffier” than observed galaxies when using the particle-based size, our synthetic observations alleviate that discrepancy within the resolution limit. The tidally transformed stream progenitors tend to be slightly larger than other satellites at the same $M_V$.

While there exist some fainter and smaller observed galaxies, our simulation resolutions are not high enough to probe this regime. Furthermore, our mock galaxies start to deviate from the observed trend around $M_V \gtrsim-6$ as they approach the simulation resolution limits ($M_\star\approx10^5\,M_\odot$ or $\approx20$ star particles). These galaxies are also the ones that comprise most of the undetectable fraction.

In Figure 6, we also compare the size-luminosity distribution of our mock field galaxies to observed dwarf galaxies in the local field taken from the LVDB. The local field dwarfs in the LVDB reside outside both the MW and M31, and within a distance of $\sim3$ Mpc. Similarly, we consider the galaxies outside the virial radius of the main (most massive) host galaxy in the isolated suite as field galaxies. As discussed in § 5.1, we use the properties as estimated from the more massive host in case of the paired suite.

The blue and orange points show the field galaxies within a distance of $D \leq1.5$ Mpc and $1.5 < D \leq3.5$ Mpc, respectively, from the observer. The blue and orange shaded regions show the 50% detectability regions for distance from the observer within 1.5 Mpc and 3.5 Mpc, respectively, taken from B. Mutlu-Pakdil et al. (2021), for LSST. The filled and open circles show the detectable and undetectable galaxies, respectively according to the 50% detectability region. As in B. Mutlu-Pakdil et al. (2021), we also classify any field galaxy above the brightest limit of the regions as detectable. Here too, we see the impact of resolution; our mock field galaxies start to deviate from observations around $(\log_{10}(r_{1/2}/\mathrm{pc}), M_{\mathrm{V}}) \approx(3, -6)$. However, almost all of these galaxies are below the 50% detectability regions. Given that our systems tend to be slightly lower in surface brightness than the observed field galaxies, we can consider our analysis as a lower limit on the number of detectable field galaxies.

Color-Magnitude Diagrams

In order to estimate the SFHs and UV luminosities at earlier redshifts, observers use the resolved color-

magnitude diagram (CMD) of a galaxy. We show example CMDs in the LSST $g$ and $r$ bands for some of our mock substructures in Figure 7. As discussed in § 3.3, these CMDs do not include LSST photometric uncertainties, extinction, or reddening, so the isochrones used to sample the stellar populations are still apparent in our results. Including photometric errors will blur these features. The top row shows the CMD constructed using the synthetic stars only in the detectable regions (see § 4), while the bottom row includes all the synthetic stars brighter than the selected LSST apparent magnitude cutoffs ($g < 25$ and $r < 24.7$). The left column is from a stream with detectable tidal tails, the middle column is from a stream progenitor with no detectable tails, and the right column is from an intact satellite. For the stream detectable with tidal tails, we include the stars within the region whose average surface brightness is $\mu_g \leq34$ mag arcsec$^{-2}$. For the stream progenitor (without detectable tidal tails) and satellites, we included all the stars within $4r_{1/2}$ as detectable. As discussed previously, this is also the region where we mask any bound progenitor of a stream. We find that the detectable parts of these substructures generally contain enough stars to map out the full CMD within LSST magnitude limits, but for some of the streams detectable as satellites we might lose the information from the fainter stars if they reside in the tidal tails. This is important when testing the accuracy of the reconstructed SFHs via CMD fitting. However, we show these CMDs only to illustrate the usefulness of mock catalogs and do not perform any CMD fitting to reconstruct the SFHs, as it is beyond the scope of this paper. In future data releases, we will include photometric errors and fit the CMDs using codes like MATCH (A. E. Dolphin 2002) to get more accurate “observed” SFHs.

Star-Formation Histories

SFHs from resolved photometric observations of the stellar populations in nearby galaxies are the key to inferring the UVLFs at $z \gtrsim6$ through the near-far technique. In order to correctly reconstruct the UVLFs, it is important to accurately reconstruct the SFHs as well. However, due to stellar evolution combined with complex merger histories, it is difficult to get the complete SFH. Furthermore, detectability of the galaxy plays a crucial role in how accurately we can reconstruct the complete SFH. Hence, we compute the SFHs from our mock observations using only detectable stars and compare them against the “true” SFH from the simulations.

We note that we do not use a CMD fitting approach like observers to reconstruct the SFHs, and hence our detectable stars encode the age and mass of their parent star particle with 100% precision. Hence, the results shown here are for best scenarios when we have accurate and precise information of the age and mass of stars.

In Figure 8, we show as an example the SFHs of the same stream, stream progenitor, and intact satellite shown in Figure 7. The blue curves show the SFHs computed directly from the entire set of simulation star particles, the orange curves include SFHs from all the mock stars within LSST magnitude limits, the green curves show the SFHs from all the stars within the detectable regions as discussed in § 5.3, and the red curves show the SFHs from all the detectable star particles which we identify by uniquely matching the detectable mock stars to their parent star particles. In the top row, we show the differential star-formation rate (SFR) as a function of age (or redshift). We renormalize the mock “observed” SFRs (orange and green) to a total mass comparable to the simulated SFH (blue), since the total mass of the mock stars, whether at present day ($m_{\mathrm{act}}$; see Appendix § 13) or at formation ($m_{\mathrm{mass}}$), is well below the actual total stellar mass of the substructure (faint orange and green curves in the top row of Figure 8). This is mainly because a substantial portion of the stellar mass falls below the magnitude limits imposed on our survey. Since we do not perform any CMD modeling, we cannot extrapolate the total stellar mass or luminosity of the simulated substructure from the mock stars alone. Instead, we calculate the total mass of the substructure using the formation mass of the star particles, which is what we use to calculate the SFH. Then, we calculate the total formation (or ZAMS) mass of the mock stars, provided by the smass quantity, in both the magnitude-limited and detection-limited cases. Mock SFRs scaled by the ratio of these two quantities are shown as the brighter orange and green curves, and are in good agreement with the “simulation” SFRs (blue). While this rescaling sets the total mass, it does not change the shape of the SFH, which depends on the fraction of each star particle’s formation mass that ends

up in detectable mock stars. Since this fraction varies strongly with the age of the star particle (see below), a single, age-independent factor cannot correct for it, and the mock SFRs overestimate later star formation and underestimate early star formation. However, this discrepancy is not prohibitively large for our calculation of the UV luminosities in most cases, as we will show in § 7. To replicate what CMD modeling using detectable stars would recover ideally, we also show the SFRs computed from only the detectable star particles (red). We find excellent agreement with the “simulation” SFRs from the entire system, particularly in case of the intact satellite, under the assumption that we get precise age and formation mass of the stars.

In the bottom row, which shows the cumulative SFH normalized to unity as a function of age, the relative shape of the “true” simulated SFH (blue) is largely recovered by the mock “observed” SFHs (orange and green). The remaining differences between the shape of the mock and simulated SFHs of the systems shown in this example are not primarily due to the missing stellar populations in the undetectable regions, as the SFHs derived from only the detectable star particles (red) closely match that of the entire system (blue). This shows that the detectable stars sample the full range of stellar ages. Instead, the differences arise mainly from our methodology to calculate the mock SFHs. The number of detectable mock stars per unit formation mass is much higher for younger star particles than for older ones. This is because in older populations, only evolved stars in the red-giant branch or near the main-sequence turnoff are detectable, whereas younger populations have more numerous, massive luminous stars. Summing the formation masses (smass) of the detectable stars in each age bin therefore underestimates earlier star formation relative to later star formation. CMD fitting of the detectable stars, which infers the total mass formed at each age while accounting for this age-dependent selection, should therefore recover the SFH of the entire system to a good extent.

Across the full set of mock catalogs, the trends depend on the classified substructure. For streams detectable with tidal tails, the mock SFHs agree well in shape with the simulated SFHs in all but a few cases. For the streams detectable as satellites, the mock SFHs are also generally representative of the entire stream. This suggests that the bound progenitors in these systems span the entire stellar population of the stream. This could be possible if the disrupted system has gone through multiple pericentric passages leading to mixed stellar populations between the bound component and the tidal tails. Therefore, such systems would provide unbiased estimates of the high-$z$ UV magnitudes. The intact satellites show the largest discrepancies, with several cases in which the shape of the mock SFH differs substantially from the simulated one (see Figure B2 for examples from another set of stream, stream progenitor, and satellite). However, when we reconstruct the SFHs from the detectable star particles (red), the satellites show excellent agreement with the simulation SFHs (blue). As Figure B2 shows, the largest discrepancies occur in systems whose star formation continues to the present day, where younger populations with massive detectable stars have much larger contributions to the SFH, relative to the older populations. This shows that our SFH reconstruction method using just the smass quantity cannot reproduce the “simulation” SFHs very well in these cases. In future work, we will derive SFHs through full CMD fitting to confirm this.

Ultraviolet Luminosity Functions

The main goal of the near-far technique is to probe the faint end of the UVLF during reionization. In A. Kundu et al. (2026), we explored how accurately we can recover the UVLFs of the proto-MW/LG at $z \gtrsim6$ by adding the progenitors of present-day disrupted galaxies along with the intact galaxies in the all-progenitors and fossil-record scenarios. The all-progenitors scenario uses the full merger history information to account for the fact that a present-day substructure can have multiple high-$z$ progenitors, while the fossil-record scenario combines the total stellar mass from all the progenitors of each present-day substructure into one. In this work, we only focus on the fossil-record scenario, since this is the main assumption of the near-far technique and can also recover faint-end slopes depending on the selection criteria of present-day galaxies (see A. Kundu et al. 2026 for details). Furthermore, we only add the stellar streams to our present-day sample of intact galaxies and ignore the phase-mixed galaxies, as we are not estimating their detectability in this work.

To calculate the rest-frame UV magnitude ($M_{\mathrm{UV}}$) at 1500 Å for the progenitor of each simulated present-day substructure, we track its constituent star particles to $z \gtrsim6$ and select the ones that already formed by that redshift. Then, using their masses and ages at that redshift, along with their metallicities, we compute the total $M_{\mathrm{UV}}$ of the progenitor, following the approach in G. Sun et al. (2023). In case of the mock observations, we follow a slightly different approach. We select all the “mock” or “detectable” stars that formed by our redshift of interest, take their formation mass (smass) and multiply by the factor that we calculated in Section 6. The rest of the approach is the same as the one we used in case of the “simulated” star particles. Once we have $M_{\mathrm{UV}}$ for the progenitor of each present-day substructure, we compute the differential UVLFs following the convex-hull volume calculation method described in A. Kundu et al. (2026).

In Figure 9, we show the reconstructed differential UVLFs at $z = 7$ from the fossil record using simulated star particles and present-day detectable mock stars. In the left panel, we consider only satellites as intact galaxies, while in the right panel we also add the field galaxies to the list of intact galaxies. We separate the two in order to check how well we can reconstruct the high-$z$ UVLFs using only the MW satellites. Furthermore, in the paired suite, the field galaxies are common to both the hosts. This can lead to a bias in the halo-to-halo scatter as these galaxies are drawn from the same distribution. The lines show the median across the 13 hosts, while the shaded regions (only shown for the detectable cases) show the 68% halo-to-halo scatter. The dashed blue and red curves show the reconstructed UVLFs from the progenitors of present-day simulated intact galaxies (satellites on the left panel; satellites and field galaxies on the right panel) with and without the streams, respectively. The solid red curve shows the reconstructed UVLF using only the detectable intact galaxies. The solid purple curve includes the streams detectable as satellites, while the solid blue curve also adds detectable streams. We compare all the reconstructed UVLFs with that of the proto-MW/LG as shown by the black dashed curve and the UVLF from the FIREbox$^{\mathrm{HR}}$ (olive) large-volume simulation (R. Feldmann et al. 2025), which we use as the benchmark for the universal UVLF during EoR. The olive shaded region shows the 68% confidence interval (CI) obtained via bootstrapping, and not a halo-to-halo scatter. The difference between the proto-MW/LG and FIREbox$^{\mathrm{HR}}$ UVLF is due to the methodological difference when calculating the distributions. We refer the reader to A. Kundu et al. (2026) for detailed discussion. Adding the streams detectable as satellites improves the normalization of the UVLF of detectable intact galaxies by a factor of $\sim1.5$. This shows that the “streams as satellites” population can contribute significantly to the reconstructed UVLF at high redshifts. Furthermore, when comparing simula-

tions with observations, these objects can be mistakenly compared with the satellite population of the simulations. The near-far reconstruction will need to account for the missing stellar populations in the tidal tails for cases where the bound components are biased tracers of the entire population. However, in § 6, we find the mock SFH for most of these systems agrees well with the simulated SFH for the entire stream. It is mostly the actual satellite galaxies that have larger differences between the mock and simulated SFHs. We emphasize that this could be due to our methodology for computing the SFHs and needs further investigation. Adding the streams with detectable tidal tails, we improve upon the normalization at $M_{\mathrm{UV}} \lesssim-10$, which is consistent with what we found in A. Kundu et al. (2026). The more massive streams are easily detectable with tidal tails and hence amplify the brighter end of the reconstructed UVLF at high redshifts. The difference in the normalization of the black and blue curves is due to the progenitors of the present-day phase-mixed galaxies and the main host galaxy, which we have not considered here. The shape of all the curves (both simulated and detectable) is qualitatively similar to that of the proto-MW/LG and the universal UVLF down to $M_{\mathrm{UV}} \approx-10$. Furthermore, adding the stream population improves the normalization by at least a factor of $\sim2$. Since LSST will detect more of these systems, they should be included while reconstructing the reionization-era UVLFs using the near-far technique.

For a quantitative comparison, we calculate the slope at the faint end of the UVLFs assuming a Schechter function (P. Schechter 1976) and performing a Bayesian modeling using the PyMC package (S. Abril-Pla et al. 2023), as described in A. Kundu et al. (2026). Figure 10 shows the faint-end slopes at $z = 7$ for both the simulated and detected cases. The points show the central estimate of the slope across all the 13 halos and the error bars indicate the bootstrapped halo-to-halo scatter (see A. Kundu et al. 2026 for details). The olive point and the shaded region show the best-fit value and uncertainty (not halo-to-halo scatter) in the faint-end slope from the FIREbox$^{\mathrm{HR}}$ UVLF. The black points show the faint-end slope of the proto-MW UVLF, while the red points show the slope from only the intact galaxy progenitors (left panel includes only satellites; right panel also includes the field galaxies). The blue points show the slope when streams are also added to the sample of present-day intact galaxies. The dashed error bars with open markers indicate the slopes estimated from the progenitors of the simulated substructures, while the solid error bars with filled markers denote the slopes calculated from the progenitors of only detectable substructures. The purple point includes detectable satellites and streams without tidal tails (streams detectable as satellites). While this is closer to the slopes measured observationally in the near-far technique, we can see that the accuracy is very similar to using only satellites. However, we see a decrease in the scatter which leads to a better constraint on the proto-MW slope. Adding the detectable streams yields more accurate estimates of both the proto-MW and universal faint-end slopes. Additionally, we reduce the halo-to-halo scatter by $\sim33\%$. As shown in the right panel, adding the detectable field galaxies to the sample of present-day intact galaxies provides similar slopes. However, our estimates are also biased by the common distribution of field galaxies in the paired hosts. Adding detectable streams from LSST can thus substantially improve both the accuracy of the recovered faint-end slope of the proto-MW UVLF during EoR and its robustness against halo-to-halo scatter.

Summary

In this work, we generate and publicly release 13 mock LSST catalogs of resolved stars in both intact and disrupted galaxies around MW-mass hosts in the FIRE-2 cosmological zoom-in simulations (7 isolated MW analogs and 6 hosts from 3 LG-like pairs) using py-ananke. The catalogs provide astrometry, photometry in the six LSST bands, stellar parameters, and chemical abundances for each synthetic star, along with indices linking them to their parent star particles (Appendix § 13). Using these catalogs, we estimate the detectability of satellites, stellar streams, and field galaxies, identified in A. Kundu et al. (2026), and test whether including detectable stellar streams in addition to the intact galaxies improves near-far reconstructions at the faint end of the UVLF during the EoR. Our results are summarized below.

1. Estimating the detectability of the tidal tails of the streams using a surface brightness limit of $\mu_g \le34$ mag arcsec$^{-2}$ and 1-year LSST magnitude limits, we find an average of $\approx29\%$ of our simulated streams with $M_\star\gtrsim5 \times10^{5}\,M_\odot$ are detectable with tidal tails (a median of three per host) across the 13 halos. Another $\approx38\%$ of the streams are detectable only through their bound progenitors and would therefore be classified as satellites. The remaining $\approx33\%$ are undetectable, and most streams with $M_\star\lesssim10^{7.5}\,M_\odot$ fall in this category.

2. We estimate the detectability of the intact galaxies (satellites and fields) and bound remnants of stellar streams (stream progenitors) using analytic detectability relations in the size–luminosity ($r_{1/2}-M_V$) plane. We find that an average of $\approx96\%$ of satellites, $\approx95\%$ of bound stream progenitors, and $\approx64\%$ of field galaxies (out to $\sim3.5$ Mpc) with $M_\star\gtrsim10^{5}\,M_\odot$ are detectable, with the numbers decreasing to $\approx59\%$, $\approx89\%$, and $\approx42\%$, respectively, for $M_\star\gtrsim10^{4.5}\,M_\odot$ owing to resolution limits. Our mock galaxies follow the observed size–luminosity relation of MW satellites and local field dwarfs for $M_V \lesssim-6$.

3. The cumulative radial distribution of detectable satellites alone underpredicts the observed MW satellite distribution, which lies within the 68% halo-to-halo scatter only beyond $\sim220$ kpc. Including streams detectable as satellites improves the radial distribution and agrees well with observations beyond $\sim150$ kpc. When we limit our sample to substructures with $M_\star\gtrsim10^{5}\,M_\odot$, the combined population agrees well with the observed distribution over the entire distance range within the halo-to-halo scatter. Our results suggest that some known MW satellites may be the bound remnants of disrupted systems. The remaining deficit within $\sim150$ kpc might be driven by the absence of an LMC analog in most of our MW-mass hosts. Detectable field galaxies reproduce the observed local field distribution out to $\sim800$ kpc, beyond which the zoom-in volumes become incomplete.

4. SFHs computed from only the detectable mock stars broadly recover the simulated SFHs for most streams and streams detectable as satellites, indicating that their bound remnants typically sample the full stellar population. However, we emphasize that our reconstructions are based on 100% pre-cise age and mass of stars, which is not possible in the observer-approach of CMD fitting. Thus, our results are optimistic reconstructions of the SFHs. Even under this assumption, discrepancies exist, the largest of which are typically found in satellites with star formation continuing to the present day. These discrepancies are mainly due to our simplified SFH reconstruction method using the formation mass and age of the mock stars, readily available in the catalogs. They are not from missing stellar populations, as SFHs from only the detectable star particles closely match the full simulated SFHs in most cases. Hence, reconstructions based on fitting the CMDs of LSST-detectable stars can recover the SFHs to a good extent, which we will test in future work.

5. The qualitative shapes of reconstructed UVLFs at $z=7$ using only detectable substructures are in good agreement with those from the simulated populations and with the proto-MW UVLF. We also find that adding the detectable streams improves the normalization of the reconstructed UVLFs — adding the streams detectable as satellites improves the normalization by a factor of $\sim1.5$, and adding all detectable streams (including streams with detectable tidal tails) by a factor of $\gtrsim2$ around $M_{\mathrm{UV}}\sim-10$.

6. From our estimation of the faint-end slopes of the UVLFs at $z=7$, we find that adding the streams with detectable tidal tails along with the intact galaxies and stream progenitors reconstructs a more accurate and precise faint-end slope of the proto-MW UVLF. Using only the intact galaxies we recover a faint-end slope steeper by $\approx0.2$ than the median proto-MW slope across the 13 halos. Adding the detectable streams reduces this offset by $\approx0.1$ and reduces the halo-to-halo scatter by $\approx33\%$.

Our results suggest that LSST-detectable substructures can feasibly be used in the near-far technique to infer the high-$z$ UVLF during the EoR, and stellar streams can improve the constraints on the reconstructions, even with 1-year magnitude limits. This opens up the possibility that proper motions for stream stars derived from LSST observations ((Ž. Ivezić et al. 2019; D. C. H. Gomes et al. 2025}) can be used to improve stream membership selections for the near-far technique, which can lead to a substantially cleaner selection. Because a substantial fraction of streams will be observed only as satellites, near-far analyses should allow for the possibility that some observed satellites are the bound remnants of larger disrupted systems, although in our simulations these remnants are generally representative of the full stellar population.

Our analysis includes some caveats. We assume perfect photometry without extinction and neglect contamination from unresolved background galaxies, which would require modeling star-galaxy separation; we expect this will have the effect of making the surface brightness cutoff slightly brighter. Our results might change depending on the actual surface brightness limit LSST will achieve. In future work, we will add photometric uncertainties to the catalogs and reconstruct SFHs through an observer-oriented approach of fitting the CMDs. We have also not estimated the detectability of the phase-mixed objects as it is beyond the scope of this work, but they are available in our mock catalogs. Adding the phase-mixed galaxies would further improve our reconstruction of the faint-end slope as we showed in A. Kundu et al. (2026). Beyond the near-far reconstructions of the high-$z$ UVLFs, these publicly available catalogs can support studies of tidal disruption, halo tracer populations such as RR Lyrae and M giants (Alvarado et al., in prep.), and the chemical abundances of intact and disrupted galaxies.

Acknowledgments

We would like to thank Robert Feldmann for providing us the FIREbox$^{\mathrm{HR}}$ UVLF. We would also like to thank Peter Ferguson for providing us the LSST resources and helpful discussions. We would also like to thank Yasmine Meziani and Pratik Gandhi for insightful comments which improved the quality of the manuscript. We acknowledge use of the marmalade cluster at the University of Pennsylvania, on which the analysis of this work was done. AK and RES acknowledge support from NSF grants AST-2007232 and AST-2307787. RES further acknowledges support from NASA awards 22-ROMAN22-0055, 22-ROMAN22-0011, and 22-ROMAN22-0013, Simons Foundation grant 1018462, and Sloan Foundation award FG-2023-20669.

Software: astropy ((Astropy Collaboration et al. 2013, 2018, 2022)), matplotlib ((J. D. Hunter 2007)), IPython ((F. Perez & B. E. Granger 2007)), NumPy ((S. Van Der Walt et al. 2011; C. R. Harris et al. 2020)), SciPy ((P. Virtanen et al. 2020)), GizmoAnalysis ((A. Wetzel & S. Garrison-Kimmel 2020a)), HaloAnalysis ((A. Wetzel & S. Garrison-Kimmel 2020b; A. R. Wetzel et al. 2016)), PyMC ((S. Abril-Pla et al. 2023)), healpy ((A. Zonca et al. 2019)), mhealpy ((I. Martinez-Castellanos et al. 2022)), scikit-learn ((F. Pedregosa et al. 2011)), py-ananke ((A. Thob et al. 2024)), NASA’s Astrophysics Data System (ADS), arXiv

Appendix

Radial distribution of stream progenitors and intact galaxies with higher stellar mass cutoff

Figure A1 shows the radial distribution of the stream progenitors and the intact satellites with $M_\star\geq10^5\,M_\odot$. Similar to Figure 3, we show the median distribution of the detectable populations across all the 13 host galaxies (solid lines). The shaded regions show the 68% halo-to-halo scatter. The thin lines show the detectable populations from each of the host galaxies. The dashed curves show the radial distribution of the simulated substructures, which includes the undetectable populations. The black curve shows the distribution of observed satellites in the MW with $M_\star\geq10^5\,M_\odot$, taken from LVDB. Similar to Figure 3, the observed MW satellite radial distribution agrees better with the distribution of combined population of detectable satellites and streams detectable as satellites, and is within the halo-to-halo scatter over the entire distance range. Additionally, we find that the distribution from the detectable systems is in excellent agreement with the simulated one for this stellar mass cutoff.

Examples of large differences between mock and simulated star-formation histories

In this section, we show a set of example SFHs where the discrepancy between the mock and simulated SFHs is large. Similar to Figure 8 we show the SFHs for a stream, a stream progenitor, and a satellite. The top row shows the SFRs as a function of age (or redshift) and the bottom row shows the SFHs normalized to unity. The “observed” SFHs from all the mock stars (orange) as well as the detectable mock stars (green) do not agree well with the “simulation” SFHs (blue). This is due to our SFH reconstruction method as discussed in § 6. The detectable mock stars span the entire age distribution of the simulated star particles which leads the SFH from detectable star particles to have excellent agreement with the simulated one. This suggests that CMD fitting from detectable LSST stars can recover the SFHs to a good extent.

Data in the mock catalogs

In this section, we describe in detail each of the fields in the mock catalogs. For each of the 13 host galaxies, we generate synthetic stars for every intact and disrupted galaxy. A given mock observation contains multiple files partitioned according to the Healpix tessellation scheme described in § 3.3. For a complete mock observation of a particular galaxy, the user needs to read all the partitioned files. A summary of the data products including their units is available in Table C1.

Data	Description	Units
Indices
partitionid	Index of the Healpix partition	$\cdots$
parentid	Array index of the parent star particle in the $z = 0$ simulation snapshot	$\cdots$
partid	Flag for a star at the exact phase-space coordinates of its parent particle (0 or 1)	$\cdots$
Phase-space Coordinates
ra, dec	Astrometric Right Ascension and Declination on the celestial sphere	deg
px, py, pz	Cartesian positions relative to the observer (LSR) position	kpc
vx, vy, vz	Cartesian velocities relative to the observer (LSR) velocity	km s$^{-1}$
rad	Radial distance to the observer	kpc
vr	Radial velocity to the observer	km s$^{-1}$
dmod	Distance modulus	$\cdots$
dform	Formation distance from the host galaxy center	kpc
pi	Parallax (see Equation (1))	mas
glat, glon	Galactic latitude and longitude on the celestial sphere	deg
mura, mudec	Astrometric proper motions along RA and DEC	mas yr$^{-1}$
mul, mub	Astrometric proper motions along Galactic longitude and latitude	mas yr$^{-1}$
Photometry
lsst_<filtername>	Intrinsic magnitude in LSST band <filtername> ($ugrizy$) in AB system	mag
lsst_<filtername>	Apparent magnitude in LSST band <filtername> ($ugrizy$) in AB system	mag
Stellar Parameters
mact	Present-day stellar mass after stellar evolution	$M_\odot$
smass	Zero-age main sequence (ZAMS) mass	$M_\odot$
mtip	Mass at the tip of the giant branch for the star’s age and metallicity	$M_\odot$
teff	Effective temperature	K
lum	Stellar luminosity	$L_\odot$
grav	$\log_{10}$ surface gravity	cm s$^{-2}$
age	$\log_{10}(\mathrm{age}/\mathrm{yr})$ from the parent star particle	$\cdots$
Chemical Abundances (derived from parent star particle)
<element>	[X/H] abundance of element X (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) relative to solar	dex
alpha	Mean of [O/H], [Mg/H], [Si/H], and [Ca/H], relative to solar	dex

Indices

Our catalog contains a number of different integers for identifying the synthetic stars and their associated generating star particle in the simulations.

• parentid - array index of the star particle in the simulation from which a given synthetic star was generated. Synthetic stars with the same parentid originated from the same parent star particle in the simulation. This index can be used to get the full properties of the star particle in the simulation snapshot at $z = 0$. We also note that the star particles in FIRE have a particle ID which is not unique across the particles and we do not use them in our pipeline. Hence, parentid just refers to the array index of the star particle in the $z = 0$ simulation snapshot.

• partid - indicates whether a synthetic star is located at the exact phase-space coordinates of its generating star particle. The first synthetic star to be generated is always assigned the coordinates of its parent star particle and has partid = 0, while the rest are assigned partid = 1.

• paritionid - indicates the partition index of the Healpix tessellation scheme. Synthetic stars belonging to the same Healpix partition have the same paritionid.

Phase-space Coordinates

We provide different types of phase-space coordinates associated with each synthetic star. We emphasize that we do not apply observational errors to these values and hence they are the true coordinates. Future data releases may contain the error convolved versions of them.

• ra, dec - the astrometric Right Ascension (RA) and Declination (DEC) of the synthetic star on the celestial sphere.

• px, py, pz - the Cartesian positions of the synthetic star relative to the observer (or LSR in this work) position.

• vx, vy, vz - the Cartesian velocities of the synthetic star relative to the observer (or LSR in this work) velocity.

• rad - the radial distance of the synthetic star to the observer.

• vr - the radial velocity of the synthetic star to the observer.

• dmod - the distance modulus of the synthetic star.

• dform - the formation distance of the synthetic star. FIRE simulations track the distance of the star particles from the center of the host galaxy at their time of formation. We provide this as an input to the pipeline and hence synthetic stars originating from the same parent star particle will have the same formation distance.

• pi - the parallax of the synthetic star. It is related to the distance modulus via

• glat, glon - the Galactic latitude and longitude of the synthetic star on the celestial sphere.

• mura, mudec - the astrometric proper motions in the direction of RA and DEC, respectively.

• mul, mub - the astrometric proper motions in the direction of Galactic longitude and latitude respectively.

Photometry

We have two types of photometric magnitudes for the LSST bands in our mock catalogs:

• lsst_<filtername>_Intrinsic - absolute AB magnitude of the synthetic star, where <filtername> is any one of the LSST ugrizy bands.

• lsst_<filtername> - apparent AB magnitude of the synthetic star. For now, this does not include galactic extinction or photometric errors, which future data releases will contain. Hence, the apparent and absolute magnitudes differ only by the distance modulus.

Stellar Parameters

We provide the values of some stellar parameters for the synthetic stars in our catalog, obtained from the isochrone models described in § 3.3. They do not have any observational errors associated with them.

• mact - present-day stellar mass of the synthetic star after stellar evolution.

• smass - mass of the synthetic star when it was on the zero-age main sequence (ZAMS).

• mtip - mass of the synthetic star at the tip of the giant branch for its given age and metallicity. All stars originating from the same parent star particle have the same mtip value.

• teff - effective temperature of the synthetic star in K.

• lum - stellar luminosity of the synthetic star in $L_{\odot}$.

• grav - $\log_{10}$ surface gravity of the synthetic star in $\mathrm{cm}\,\mathrm{s}^{-2}$.

• age - $\log_{10}(\mathrm{age}/\mathrm{yr})$, obtained directly from the age of the parent star particle. All stars originating from the same parent star particle have the same age.

Chemical Abundances

FIRE simulations track the abundances of 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) through the stellar yield tables and subgrid turbulent metal diffusion models (see P. F. Hopkins et al. (2018a) for details). We forward these abundances directly to each synthetic star originating from its parent star particle. Thus, all synthetic stars spawning from the same parent star particle share the same abundances. All abundances are given relative to hydrogen and normalized to solar values ((M. Asplund et al. 2009)). We also provide the alpha abundances, defined as the mean abundance using O, Mg, Si, and Ca, relative to solar. The list of chemical abundances is reported in Table C1.

Note— Abundances are relative to hydrogen and normalized to the solar values of M. Asplund et al. (2009). Elements He, C, N, O, Ne, Mg, Si, S, Ca, & Fe are respectively identified in the data files as helium, carbon, nitrogen, oxygen, neon, magnesium, silicon, sulphur, calcium, & feh. Extinction, reddening, and photometric errors are not included in the apparent magnitudes.