Introduction

The satellites of the Milky Way and its neighbors in the Local Group, thanks to their proximity, have often served as our basis of understanding satellite and dwarf galaxy formation and evolution. In the past decades there has been an explosion in our understanding of satellites around our own Milky Way (Mateo 1998;Koposov et al. 2008;Simon 2019;Drlica-Wagner et al. 2020, and references within) and Andromeda (Ibata et al. 2014;Martin et al. 2016;McConnachie et al. 2018, and references within). Further, in the age of ultrafaint galaxy detection, the low surface brightness end of the Milky Way's satellite distribution continues to grow (e.g., Drlica-Wagner et al. 2015;Koposov et al. 2015;Simon 2019). As we continue to discover fainter objects nearby, the question of the Milky Way's uniqueness becomes an important one. Applying what we learn locally to the Universe at large would not be appropriate if the Local Group could be considered "atypical."

To test for any potential discrepancy, surveys such as the "Satellites Around Galactic Analogs" (SAGA; Geha et al. 2017;Mao et al. 2021) and "Exploration of Local VolumE Satellites" (ELVES; Carlsten et al. 2020Carlsten et al. , 2021aCarlsten et al. , 2021bCarlsten et al. , 2022) ) study the satellite distributions of galaxies similar to our own, placing the Milky Way in a broader, cosmological context. The SAGA survey is an ongoing effort to compile spectroscopically complete satellite luminosity functions of 100 Milky Way analogs with distances between 20 and 40 Mpc, providing vastly improved statistics for the bright end of these satellite distributions (down to M R = -12.3). In complement to the SAGA survey's probing of distant Milky Way-like systems, the ELVES survey seeks to map the satellite distributions of the hosts within the Local Volume fully (<12 Mpc) down to M V = -9.

Working in tandem, SAGA and ELVES will provide a better understanding of both what a "typical" Milky Way-like halo will look like and what influences an environment like the Local Volume can impart. The SAGA survey has found that the luminosity function of the Milky Way is consistent with their observations of other systems, but that the host-to-host scatter in the number of satellites is large (Mao et al. 2021). SAGA also finds that the total number of satellites in a system correlates with the host's K-band luminosity. Similar to SAGA, the ELVES survey finds that satellite abundance correlates with host mass and that the Milky Way is typical for its mass. However, Carlsten et al. (2021b) find that the observed luminosity functions of local hosts are typically "flatter" than predicted by the cosmological model; the stellar-to-halo mass relation tends to underpredict bright satellites and overpredict faint ones, a result found also by Geha et al. (2017). These results highlight the power of a larger sample of galaxies and their satellites to provide context for understanding satellite dwarf galaxies.

One of the most interesting discrepancies to be highlighted so far is that the quenched fraction of Local Group (Milky Way and Andromeda) satellites is not in agreement with SAGA's results: the SAGA sample exhibits lower quenched fractions than those found in the Local Group. On the other hand, the ELVES survey finds higher quenched fractions among the Local Volume than in the SAGA sample, though still not as high as the Local Group. Although Mao et al. (2021) carefully attempt to quantify incompleteness in the SAGA survey, it remains an open question whether SAGA may be missing faint, red, or low surface brightness satellites which would be predominantly quenched (Carlsten et al. 2022;Font et al. 2022), or whether the Local Group is a true outlier in terms of quenched satellite fraction.

In general, various simulations of Milky Way-like galaxies tend to find good agreement in their resulting quenched satellite fractions, lying somewhere between the Local Group and Local Volume fractions (Akins et al. 2021;Engler et al. 2021;Karunakaran et al. 2021;Samuel et al. 2022). These simulations generally find that galaxies that infall into a host Milky Way with stellar masses above M * ∼ 10 8 M e are better able to retain their gas and continue star forming for extended periods. On the other hand, galaxies with stellar masses below 10 8 M e instead tend to experience ram pressure stripping that strips gas and quenches their star formation (SF), and the quenching timescales can often be quite short (<2 Gyr; Wetzel et al. 2015;Simpson et al. 2018;Simons et al. 2020;Akins et al. 2021). These results lead to high predicted quenched satellite fractions as luminosity decreases.

On the theoretical front, many analyses use zoom-in simulations of a handful of Milky Way analogs (e.g., Akins et al. 2021;Samuel et al. 2022), though Font et al. (2022) use the ARTEMIS suite of 24 cosmological Milky Way-mass zooms to interpret the ELVES and SAGA results. Font et al. (2022) find that applying a surface brightness limit to the ARTEMIS satellites can bring the quenched fractions and radial distributions into line with the SAGA results, suggesting that SAGA is missing faint surface brightness galaxies. Fainter surface brightnesses correlate with more quenching at a fixed luminosity in ARTEMIS, and thus bias the SAGA results if true. On the other hand, Engler et al. (2023) found that a surface brightness cut could not bring the TNG50 satellite quenched fractions fully into agreement with SAGA, though it did bring the simulation and observational results more into line. Engler et al. (2023) were able to use TNG50, a 50 Mpc on-a-side uniform cosmological volume, to study a larger sample of Milky Way analogs and look for statistical trends. In this work, we use ROMULUS25, a 25 Mpc on-a-side uniform cosmological volume with comparable resolution to TNG50, to study similar trends. We particularly focus on the questions of how host mass and large-scale environment impact both satellite counts and quenched fractions for our simulated Milky Way analogs.

The paper is outlined as follows. We begin in Section 2 by describing the ROMULUS25 simulation. In Section 3 we outline our various methods for identifying Milky Way analogs, as well as their satellites. In Section 4 we present our primary results, focusing on the general size of the satellite populations and their quenched fractions. We then discuss and summarize our results in Sections 5 and 6, respectively.

Simulation

For this work, we use the ROMULUS25 simulation (Tremmel et al. 2017). ROMULUS25 was run with CHANGA (Menon et al. 2015), which includes standard physics modules previously used in GASOLINE (Wadsley et al. 2004(Wadsley et al. , 2008(Wadsley et al. , 2017) ) such as a cosmic UV background (Haardt & Madau 2012) including selfshielding (Pontzen et al. 2008), SF, "blastwave" supernova (SN) feedback (Stinson et al. 2006), and low-temperature metal cooling (Bromm et al. 2001). CHANGA implements an updated smoothed particle hydrodynamics (SPH) routine that uses a geometric mean density in the SPH force expression, allowing for the accurate simulation of shearing flows with Kelvin-Helmholtz instabilities (Wadsley et al. 2017). Finally, a timedependent artificial viscosity and an on-the-fly time-step adjustment (Saitoh & Makino 2009) system allow for more realistic treatment of weak and strong shocks (Wadsley et al. 2017).

ROMULUS25 assumes a Λ cold dark matter (ΛCDM) model with cosmological parameter values following results from Planck (Ω 0 = 0.3086, Λ = 0.6914, h = 0.6777, and σ 8 = 0.8288; Planck Collaboration et al. 2016). The simulation has a Plummer equivalent force softening of 250 pc (a spline softening of 350 pc is used, which converges to a Newtonian force at 700 pc). Unlike many similar cosmological runs, the dark matter particles were oversampled relative to gas particles, such that the simulation was run with initially 3.375 times more dark matter particles than gas. This increased dark matter resolution allows for the ability to track the dynamics of supermassive black holes within galaxies (Tremmel et al. 2015). The result is a dark matter particle mass of 3.39 × 10 5 M e and gas particle mass of 2.12 × 10 5 M e . These relatively low dark matter particle masses decrease numerical effects resulting from two-body relaxation and energy equipartition, which occur when particles have significantly different masses, both of which can affect the structure of simulated galaxies (e.g.,Ludlow et al. 2019). ROMULUS25 has been shown to reproduce important galaxy and supermassive black hole scaling relations (Tremmel et al. 2017;Ricarte et al.2019;Sharma et al. 2022aSharma et al. , 2022b)).

Star Formation and Gas Cooling

Gas cooling at low temperatures is regulated by metal abundance as in Guedes et al. (2011), as well as SPHs that include both thermal and metal diffusion as described in Shen et al. (2010) and Governato et al. (2015; thermal and metal diffusion coefficients set to 0.3, see Tremmel et al. 2017Tremmel et al. , 2019 for an in-depth discussion). SF and associated feedback from SNe are crucial processes that require subgrid models in cosmological simulations like ROMULUS25. Following Stinson et al. (2006), SF is regulated with parameters that encode SF efficiency in dense gas, couple SN energy to the interstellar medium (ISM), and specify the physical conditions required for SF. These parameters were calibrated using several dozen zoom-in simulations of dwarf to Milky Way-mass galaxies (Tremmel et al. 2017) and are as follows. a gas particle that has a dynamical time t dyn ( ) () = --D p m m e 1 . 1 c t t gas star SF dyn 2. The fraction of SN energy coupled to the ISM, =  0.75 SN . 3. Minimum density, n å = 0.2 cm -3 , and maximum temperature, T å = 10 4 K, thresholds beyond which cold gas is allowed to form stars.

Star particles form with a mass of 6 × 10 4 M e , or 30% the initial gas particle mass. ROMULUS25 assumes a Kroupa initial mass function (Kroupa 2001) with associated metal yields and SN rates. Feedback from SNe uses the "blastwave" implementation (Stinson et al. 2006), with thermal energy injection and a cooling shutoff period approximating the "blastwave" phase of SN ejecta when cooling is inefficient.

Halo Identification

Amiga Halo Finder (AHF; Knollmann & Knebe 2009) was applied to ROMULUS25 to identify dark matter halos, subhalos, and the baryonic content within. AHF uses a spherical top-hat collapse technique (Bryan & Norman 1998) to calculate each halo's virial radius (R vir ) and mass (M vir ). Halos are considered resolved if their virial mass is at least 3 × 10 9 M e at z = 0. This corresponds to a dark matter particle count of ∼10 4 , and a stellar mass of at least 10 7 M e (star particle count of ∼150). Following Munshi et al. (2013), stellar masses were scaled by a factor of 0.5 as a photometric correction for a better comparison to values inferred from typical observational techniques, and all magnitudes use the Vega zero-point.

Analog and Satellite Identification

There is no concrete definition of what constitutes a Milky Way analog; observational surveys like SAGA and ELVES make sample cuts using K-band magnitudes as proxies for stellar mass, while simulations have access to more exact values for halo properties such as stellar mass and virial radius. In this work, we select samples of Milky Way analogs according to three different criteria sets in order to test if the selection criteria can influence the resultant satellite distribution. Our samples are defined as follows.

1. A general M vir restriction: any halo where 10 11.5 < M vir /M e < 10 12.5 . 2. A general M * restriction: any galaxy where 10 10 < M * /M e < 10 11 . This corresponds to the host stellar mass range outlined in Section 2.1.2 of SAGA II (Mao et al. 2021). A stellar mass of 10 10 M e also corresponds to the lower limit of the ELVES survey (Carlsten et al. 2022).

3. An M k +environmental restriction: any galaxy where -24.6 < M K < -23. Additionally, no neighbor within 300 kpc can have M K,neighbor < M K,MW -1.6. This corresponds to the K-band magnitude cut and environmental restrictions from SAGA II (Mao et al. 2021).

We also explore two different ways to identify a satellite galaxy. First, we consider galaxies within the host's virial radius down to a stellar mass of 10 7 M e , the resolution limit for ROMULUS25. This corresponds to amagnitude limit of M R ≈ -12.6. We note that the magnitude limit for the SAGA survey is M R = -12.3 (though they have four satellites below this limit; see Mao et al. 2021), so our samples do not probe the lowest-mass regions of the SAGA or ELVES sample spaces. In addition, we perform a selection where satellites are identified by being within 300 kpc of a Milky Way analog, rather than the analog's virial radius, as a more direct comparison to the SAGA and ELVES surveys. We note, however, that these surveys use 2D projected distances while in this work we use true 3D distances. In the event that a satellite is hosted by multiple analogs, it is ascribed to the most massive host. Any satellites that fall into the criteria of a Milky Way analog are not included in the satellite distribution. As a final step, any analogs that host a "satellite" more massive than themselves are removed from consideration. This cut is partially responsible for the slight variation in the number of Milky Way analogs under the same criteria when switching between the R vir and 300 kpc satellite identifications. Our sample of Milky Way analogs and satellites are summarized for each criteria set in Table 1.

Figure 1 shows the three Milky Way analog samples that we focus on in this work: M vir and M * with R vir and M K with 300 kpc. While the samples largely overlap, we find that none of them are simple subsets of the others. As they approach the boundaries of the selection cuts, the samples diverge from one another. For example, the stellar mass sample probes virial masses below the virial mass cut, and vice versa. This is the result of natural scatter within the stellar-halo mass relation, which was shown in Tremmel et al. (2017) to match observations (Moster et al. 2013;Kravtsov et al. 2018). Within the overlapping regions of the criteria, there are galaxies considered Milky Way analogs in some samples but not others. This occurs as a result of the environmental criteria in the SAGA sample, which could remove analogs that are still within the K-band magnitude limits.

In Figure 2, we compare the normalized distributions of hosts and satellites from our largest sample, M * with simulated R vir (in order to encapsulate the full magnitude range of our samples), to data from SAGA II and ELVES (Mao et al. 2021;Carlsten et al. 2022). We note that the ELVES satellites are weighted according to their likelihood estimates (P sat in Table 9 of Carlsten et al. 2022), so each satellite adds its likelihood as a count rather than 1. In panel (a), we see that our hosts' span in K-magnitude space matches well with the ELVES sample, while the SAGA II sample (by definition) resides in -24.6 < M K < -23. The peaks of the host distributions are in good agreement as well, though we note our peak is at a slightly dimmer magnitude than the observational data. In panel (b), we see that our satellite distribution is in very good agreement with the SAGA II data, though we have an interesting lack of satellites at M K ≈ -17. The ELVES data probe much dimmer satellites (due to the difference in observational limits), but when only considering satellites brighter than M K = -12, the ELVES sample is still more concentrated at low-mass satellites when compared to SAGA II and ROMULUS25. This is consistent with ELVES finding steeper luminosity functions (fewer high-mass satellites and more low-mass ones) in their sample when compared to SAGA, and might also contribute to the different quenched fractions found by the two surveys (see Section 4.2 for discussion).

Results

Figure 3 shows the V-band satellite luminosity function for our sample of Milky Way analogs alongside data from the Milky Way and several Milky Way-like systems. The outer gray region outlines the space occupied by our M * within R vir sample, while the black line and inner dark gray region indicate the mean and standard deviation. The Milky Way and M31 data are taken from Geha et al. (2017). The NGC 4258 and NGC 4631 data were taken from Carlsten et al. (2021b), the M94 data from Smercina et al. (2018), and the M101 data from Bennet et al. (2019). We find that our sample of Milky Way analogs is in good agreement with these observations. We note that the space occupied by our sample remains largely unaffected when changing the Milky Way analog criteria.

Host Effects on Satellite Accumulation

In order to study how the physical properties of our Milky Way analogs affect their satellite populations, we separated our sample according to their mass and environment. Figure 4 shows the average number of satellites hosted by the Milky Way analogs where the analogs are binned according to their stellar mass and minimum distance to a Milky Way-sized or larger halo, hereafter D MW+ . In calculating D MW+ , we consider the closest galaxy outside the system of the analog (i.e., not a satellite) that exceeds the minimum criteria of Milky Way analog under the given criteria. The text in each bin details N (the number of analogs in that bin) and σ (the standard deviation of the number of satellites hosted by analogs in that bin). A plot is shown for both our M vir with simulated R vir (left) and SAGA II comparison (right) samples. In all of our samples, the number of hosted satellites appears to increase with host mass, and potentially with decreasing D MW+ . However, this latter trend cannot be verified by eye as the box size of ROMULUS25 yields a lack of data in the upper regions of this plot (i.e., highly isolated hosts), so the apparent trend is not statistically significant.

While these macroscopic trends are present across all of our simulated samples, there are some notable differences in the distributions. We see that while the M vir definition includes analogs at a lower stellar mass, the number of analogs below M * = 10 10 M e is much larger in the SAGA II sample. Additionally, in the higher-mass bins there is fluctuation in both the number of analogs and hosted satellites due to the changing of the satellite selection radius from R vir to 300 kpc.

In an effort to quantify the "by-eye" trends seen in Figure 4, we looked at the specific frequency of the number of satellites hosted by our Milky Way analogs, S N , normalized to their mass and environment. We use the following specific frequency equations adapted from Harris & van den Bergh (1981):

Here, N sat is the number of satellites hosted by the Milky Way analog, D is D MW+ in units of Mpc, and M is log(M * /M e ). The normalization values of 1.5 Mpc and 10.3 were chosen to be roughly the averages of the M * with simulated R vir sample.

Figure 5 shows the specific frequencies normalized to mass and environment for our M vir and M * with simulated R vir sample, as well as our SAGA II comparison sample. In looking at the trend with mass, the S N values consistently increase with the stellar masses of the Milky Way analogs. These results, which are present in all of our Milky Way analog samples, indicate that stellar mass exerts a large influence on satellite accumulation. The SAGA and ELVES surveys both observe this trend of satellite abundance increasing with host mass, though the trends they find are slightly weaker than ours (see Section 5.1 for discussion). Further, a study of seven nearby Milky Way-like systems with the Hyper Suprime-Cam on the Subaru telescope observes this trend as well (Nashimoto et al. 2022). The trend of satellite abundance with host mass was also found by Font et al. (2021) using the ARTEMIS suite of zoomin simulations (Font et al. 2020), and by Engler et al. (2021) using the TNG50 simulation.

In looking at the trend with environment, Figure 5(b), we see some interesting behavior. The S N values increase somewhat linearly until D MW+ ≈ 3.5 Mpc, where future points go either to zero or extreme outliers. This would suggest that N sat increases as hosts become more isolated, but we note that a majority of our hosts (∼60%-70%) have D MW+ < 2 Mpc, so beyond this distance our samples get increasingly small, resulting in the large error bars and stochasticity of the higher-D MW+ points. Thus, we see no definitive trend of satellite accumulation with environment, though one might become present with a larger sample of more isolated hosts. However, we do not believe that we can fully rule out an environmental impact on satellite accumulation through a measurement of specific frequency. Figure 6 shows the average number of satellites hosted by our analogs when binned by their D MW+ measurement. We see that if we split our analogs into subsamples of "pairs" and "isolated" based on having D MW+ < 1 Mpc or >1 Mpc, the average numbers of hosted satellites, plotted as red diamonds, are notably different between the subsamples. This difference, which is present in all of our samples except M K +Env. with simulated R vir , is driven primarily by the low number of satellites hosted by analogs in the ∼3-5 Mpc bin. Although this is where our sample tapers off in D MW+ space, we can see that the number of analogs in this bin is certainly nonnegligible, and may be hinting at a strong environmental impact on satellite accumulation in more extreme isolation. The stochasticity of our trends in Figure 5(b) compared to the more direct information from Figure 6 leads us to believe that normalizing specific frequency to such an extremely variable parameter (in this case, largescale environment) does not yield a reliable measurement.

Host Effects on Satellite Quenching

In addition to studying the number of satellites hosted by our analogs, we also analyzed the quenched fraction of the satellites. When studying the quenched fraction ( f q ), we only consider satellites with a stellar mass of at least 10 8 M e , as ROMULUS25 may be subject to numerical overquenching below this mass (Wright et al. 2021). A galaxy is considered quenched if its instantaneous specific star formation rate (sSFR) is below 10 -11 yr -1 . The instantaneous sSFR is calculated using the expected SFR from gas particles meeting the temperature and density thresholds for SF given in Section 2.1. In Figure 7, we show our quenched fractions as a function of host K-band magnitude for the M K +Env. with 300 kpc satellite selection (our SAGA II comparison sample), and compare our results to data from the SAGA and ELVES surveys (Mao et al. 2021;Carlsten et al. 2022). For a direct comparison, we only consider SAGA and ELVES satellites with stellar masses above 10 8 M e . We note, however, that the SAGA and ELVES surveys' methods of determining quenching are different than ours: SAGA considers a satellite quenched if it lacks strong Hα emission (equivalent width (EW) of Hα < 2 Å) and ELVES considers a satellite quenched if it exhibits an early-type morphology, i.e., not exhibiting clear star-forming structures such as blue clumps or dust lanes (this correlates with color as well; see Carlsten et al. 2021a for an indepth discussion). Our sSFR quenched definition was shown (see Sharma et al. 2022a) to yield a good match to galaxies identified observationally as quenched using EW[Hα] < 2 Å and * > + D M 4000 0.6 0.1 log n 10

(as in Geha et al. 2012). While all three samples show the quenched fractions increasing with host brightness, our simulated sample exhibits slightly larger quenched fractions than the observational surveys, with the exception of the lowest-mass bin where the difference becomes significant (see Section 5.2 for discussion). The SAGA and ELVES data are in very good agreement up to the brightest magnitude bin, where the sample sizes are only one host for ELVES (M31) and two hosts for SAGA (NGC 5792 and NGC 7541). This agreement within the highmass satellite subset is interesting, as the SAGA and ELVES quenched fractions are quite different when considering their full samples. Carlsten et al. (2022) find that the quenched fractions of the Local Volume are significantly higher than the SAGA sample (their Figures 11 and 12), particularly in the low-mass satellite regime. In Figure 2(a), we see that the ELVES survey contains a much larger number of faint satellites when compared to SAGA, but also that ELVES hosts (along with those of ROMULUS25) probe fainter magnitudes.

In studying the ARTEMIS simulations, Font et al. (2022) found that the SAGA detection methods may be preferentially selecting star-forming or recently quenched satellites near their completeness limit, missing a notable population of quenched dwarfs. This detection bias could explain the difference between SAGA and ELVES low-mass satellites, both the abundances and quenched fractions. Following Font et al. (2022), in Figure 2(b) we apply an additional cut to our SAGA II comparison sample by requiring satellites to have μ eff,r < 25 mag arcsecond -2 . As in the ARTEMIS simulations, we find that this cut lowers the resultant quenched fractions, and brings our results (particularly the middle bins) into excellent agreement with SAGA and ELVES.

To quantify the trend of quenched fraction with mass seen in Figure 7, and to search for a trend with environment, we again used the specific frequency in Equation (2) with N sat replaced by f q . Figure 8(a) shows our quenched fraction specific frequencies for the M vir and M * with simulated R vir samples, as well as our SAGA II comparison sample. We find that, as with the number of hosted satellites, we see a trend of S N with host mass, indicating that larger hosts are expected to yield higher quenched fractions. However, we note that this trend is not as strong as the one seen in Figure 5(a). Interestingly, when applying the satellite surface brightness criteria in Figure 8(b), we see that our trend of quenched fraction with host mass is effectively erased. As the high-mass end of Figure 8(b) is strongly affected by this surface brightness cut, it seems that the preferentially quenched satellites below this threshold are more common in higher-mass hosts, which is consistent with Figure 8(a) implying a larger number of quenched galaxies in this regime. Our results agree with Engler et al. (2023) who, using the TNG50 run from the IllustrisTNG simulations, found that massive hosts exhibit systematically larger satellite quenched fractions. Further, Engler et al. (2023) found no difference in quenching between isolated and paired analogs when considering satellites within 300 kpc of their host (see Section 5.3 for discussion).

To look for a trend with environment, we examined the quenched fractions of our systems plotted against D MW+ . Figure 9 bins our analogs in D MW+ space, and shows the quenched fraction of all satellites hosted by analogs in each bin. Again, pairs are identified as having D MW+ < 1 Mpc. The figure shows data for our M * and M K +Env. samples (both R vir and 300 kpc), where the color of each point represents the number of satellites in each bin. We find that the average quenched fraction is higher among pairs than isolated analogs, though the magnitude of this difference is not ubiquitous across our samples (see Section 5.3 for discussion). In fact, while not shown in the plot, the M vir criteria exhibit no notable difference in the average quenched fractions of paired and isolated analogs. We also find that in the switch from 300 kpc to R vir when identifying satellites, hosts typically have either the same or lower quenched fractions and a higher satellite count. This indicates that restricting satellites to within 300 kpc for this host range is more likely to exclude satellites, and that the satellites beyond 300 kpc are predominantly star forming; though this is still only when considering satellites with M * > 10 8 M e . In applying the surface brightness cut in subplots (c) and (d) we find that, while the resultant averaged quenched fractions are lower, the difference between the isolated group and pairs notably increases in all samples.

Similar to Figure 6, this difference between our "pairs" and "isolated" subsamples is largely driven by the most isolated bin, where the quenched fractions are extremely low. This, again, could be alluding to a strong environmental effect on satellite quenching in the highly isolated regime that we do not quite capture in ROMULUS25. This result is at odds with the aforementioned study of TNG50 by Engler et al. (2023), as well as the study of FIRE-2 by Samuel et al. (2022), both of whom find no difference in quenching between isolated and paired analogs. However, the simulations studied in Samuel et al. (2022) are zoom-in simulations, and thus would not capture the large-scale isolation in which we begin to see differences between paired and isolated hosts.

While this work does not consider signs of conformity within satellites, we did perform some cursory analysis on how the quenched fractions and number of hosted satellites relate to the instantaneous SFRs of our analogs. We find that the quenched exhibit no notable trend with host SFR. In looking at the average number of hosted satellites, though, we do see a correlation of more populated systems having a higher host SFR. However, we believe this is just a reflection of the number satellites increasing with host mass, as scaling to sSFR largely removes this correlation.

Discussion

In Section 3, we discussed the various methods by which we identified Milky Way analogs and satellites. While shifting between these definitions has no effect on our conclusions, there are subtle impacts worth noting.

Satellites within R vir versus 300 kpc

In Figures 4 and 5 we showed host stellar mass to be a driving factor in satellite accumulation, but this trend is less prominent when using our SAGA II comparison sample. This appears to be the result of using 300 kpc to identify satellites, not the selection on K-band magnitude, as our M K +Env. with R vir sample actually exhibits the strongest trend. In fact, identifying satellites via a 300 kpc selection rather than R vir reduces the strength of the mass trend in all criteria (though the trend is still prominent). The weakening of the trends is the result of analogs in the high-mass regime (where the trends manifest), which have virial radii larger than 300 kpc and exclude satellites in this shift to 300 kpc. This shift is in agreement with the ARTEMIS simulations (Font et al. 2021), in which satellite abundance trends strongly with host mass, but the trend is weakened when SAGA observation selection criteria are applied (M r,sat < -12, μ eff,r < 25 mag arcsecond -2 , and within 300 kpc of the host).

When considering quenched fractions, our choice of satellite selection radius also seems to have a noticeable effect on our M K +Env. sample. In Figure 9(b), the switch from R vir to 300 kpc typically raises the quenched fraction while lowering the number of satellites (with the notable exception of the least isolated bin). Thus, within the context of satellites with M * > 10 8 M e , it seems applying a satellite cut of 300 kpc to the K-band magnitude analog selection primarily removes starforming satellites from massive hosts, and biasing the global quenched fraction high.

Since the 300 kpc selection results in a more centrally located satellite population, it is likely that these satellites had an earlier infall time and underwent more ram pressure stripping when compared to satellites near or beyond 300 kpc from the host. This effect is present in Figure 8 as well, wherein the SAGA II comparison sample exhibits the strongest trend of quenched fraction with host mass. These results are consistent with those from the TNG50 simulation (Engler et al. 2023), another large-volume, uniform-resolution simulation with comparable resolution to ROMULUS25.

Quenched Fraction Discrepancy

The shift from R vir to 300 kpc, however, does not explain why our quenched fractions are higher than those of SAGA and ELVES (Figure 7). Our SAGA II comparison sample uses 300 kpc as a selection radius, and our results indicate that if SAGA and ELVES had access to the virial radii of their hosts, their quenched fractions would be lower. Donnari et al. (2021) find that the adopted definition of quenching and using 2D projected distances can both affect the resultant quenched fractions. Notably, the quenched fractions of ROMULUS25 are in better agreement with the observations when satellites with μ eff,r < 25 mag arcsecond -2 are removed, in agreement with Font et al. (2022). The exception is the faintest bin, where a key factor may be the resolution of ROMULUS25. The lower resolution of the volume is unable to resolve a multiphase ISM, i.e., there is no extremely dense gas (Tremmel et al. 2019, 2020, andreferences within). Thus, all of the gas is "puffy" and overly susceptible to ram pressure stripping and quenching. Dickey et al. (2021) found that large-scale cosmological simulations overquench isolated galaxies below M * = 10 9 M e when compared to the Sloan Digital Sky Survey. The authors attribute this to overefficient feedback, which is typically tuned to recreate quenched fractions found in the Local Volume. In looking at ROMULUS25, Sharma et al. (2022a), also found that isolated dwarfs exhibit a higher quiescent fraction when compared to observations, but that this can be entirely attributed to the presence of massive black holes and their feedback. Although we are not studying isolated dwarfs in this work, it is still likely that these feedback properties are influencing our results. We note, however, that there are only six satellites in our SAGA II comparison sample with black holes and M * > 10 8 M e so this does not notably affect our results.

Isolated versus Paired Hosts

Figure 9 suggests that pairs exhibit higher quenched fractions than more isolated analogs, but there are some caveats preventing us from making a more robust statement about the environment's effect on the quenched fraction. First, we are only considering satellites with M * > 10 8 M e . Within our SAGA II sample, this is only ∼56% of our total satellite population and they are hosted by ∼72% of our Milky Way analogs with a nonzero satellite count (or ∼49% of all Milky Way analogs), so a large section of our population is being removed. Second, our simulation box size prevents us from having a large sample of highly isolated analogs; only ∼14% of our SAGA II sample analogs have D MW > 3 Mpc. Finally, by ignoring low-mass satellites, we are looking at the quenched fractions of several systems with few satellites (only one or two satellites). Around 43% of the high-mass satellite-hosting analogs in our SAGA II sample contain only one satellite above our resolution limit, so their quenched fractions can only occupy the extremes of 0 and 1, and in Figures 8 and 9 these systems are being considered equally alongside systems with as many as eight high-mass satellites (though the binning in Figure 9 should alleviate this issue). These combined effects yield a sample that is lacking low-mass satellites (and thus the analogs' full satellite distributions) as well as highly isolated hosts, making it difficult for us to extrapolate our results to the Universe at large.

Recently, Engler et al. (2023;TNG50) and Samuel et al. (2022;FIRE-2) found no difference in the satellite quenched fractions of paired and isolated hosts in their simulations. Further, Garrison-Kimmel et al. (2019;FIRE-2) found that the satellites of isolated Milky Way-mass galaxies have nearly identical SF histories to satellites of Milky Way analogs in Local Group-like pairs. However, these results were only considering satellites within 300 kpc of the host. In looking further out to 300-1000 kpc, Engler et al. (2023) find that paired, Local Group-like hosts exhibit significantly larger quenched fractions than their isolated counterparts.

Conclusions

Using the ROMULUS25 simulation, we have created various samples of Milky Way analogs along with their satellite distributions. We explored the role of host mass and environment on satellite numbers and quenched fractions. Our results can be summarized as follows.

1. When testing various criteria for defining a Milky Way analog, from more theoretically motivated (M vir ) to more observationally motivated (M å and SAGA-like), we find that the resultant samples do not fully overlap. Within the overlapping regions, galaxies may also be defined as analogs in one sample but not in another due to environmental criteria (see Table 1 and Figure 1). 2. The number of satellites hosted by a Milky Way analog increases predominantly with host stellar mass, while environment may have a significant impact in high isolation (D MW+ > 3 Mpc; see Figures 456). 3. The quenched fraction (for satellites with M * > 10 8 M e )

of our analogs increases with host mass (see Figures 7(a) and 8(b)), but applying a surface brightness cut to satellites can erase this trend (see Figure 8(d)). 4. Being in a pair may yield higher satellite quenched fractions, but it is hard to draw statistically robust results given the small volume of ROMULUS25 and the fact that we can only study satellites down to M * = 10 8 M e to avoid numerical overquenching (see Figure 9).

We find that the distributions of both the Milky Way and M31 are well explained by our sample, with M31 being at the highly populated edge of our sample space. This is in agreement with the SAGA and ELVES surveys, where ELVES found the Local Volume to be slightly more populated and exhibiting a steeper luminosity function when compared to the full SAGA sample. Additionally, we are in agreement with ELVES in finding that the stellar mass of a Milky Way analog seems to be the dominant factor in both the number of hosted satellites and the number of quenched satellites. Interestingly, in our study of quenching, we find that the SAGA and ELVES results are in good agreement for satellites with M * > 10 8 M e , suggesting that their discrepancy in quenched fraction comes from lower-mass satellites, which we are unable to probe here due to numerical effects that artificially quench simulated galaxies. However, our results support the notion put forward in Font et al. (2022) that SAGA is missing a large population of low surface brightness satellites near its detection limit that are preferentially quenched.