Gas and Star Formation in Satellites of Milky Way Analogs

We have imaged the entirety of eight (plus one partial) Milky Way (MW)–like satellite systems, a total of 42 (45) satellites, from the Satellites Around Galactic Analogs II catalog in both Hα and H i with the Canada–France–Hawaii Telescope and the Jansky Very Large Array. In these eight systems we have identified four cases where a satellite appears to be currently undergoing ram pressure stripping (RPS) as its H i gas collides with the circumgalactic medium (CGM) of its host. We also see a clear suppression of gas fraction (M HI/M *) with decreasing (projected) satellite–host separation—to our knowledge, the first time this has been observed in a sample of MW-like systems. Comparisons to the Auriga, A Project Of Simulating The Local Environment, and TNG50 cosmological zoom-in simulations show consistent global behavior, but they systematically underpredict gas fractions across all satellites by roughly 0.5 dex. Using a simplistic RPS model, we estimate the average peak CGM density that satellites in these systems have encountered to be logρcgm/gcm−3≈−27.3 . Furthermore, we see tentative evidence that these satellites are following a specific star formation rate to gas fraction relation that is distinct from field galaxies. Finally, we detect one new gas-rich satellite in the UGC 903 system with an optical size and surface brightness meeting the standard criteria to be considered an ultra-diffuse galaxy.


INTRODUCTION
In recent years two major efforts, the Satellites Around Galactic Analogs (SAGA, Geha et al. 2017;Mao et al. 2021) and the Exploration of Local VolumE Satel-lites (ELVES, Carlsten et al. 2022) projects, have built up consistently defined, large samples of satellite galaxies around tens of hosts analogous to the Milky Way (MW) for the first time, significantly extending earlier works (e.g.Zaritsky et al. 1993Zaritsky et al. , 1997;;Crnojević et al. 2016;Bennet et al. 2019;Carlsten et al. 2020).These new samples have opened the door to statistical studies of satellites in MW-mass systems other than those in the Local Group (LG), which until now have been almost the sole point of detailed comparison with simulated systems in this mass regime.However, in the case of satellite quenching, the removal of gas and the permanent cessation of star formation (SF), SAGA and ELVES do not necessarily align with findings from the LG, or perhaps even with each other.
In the LG virtually all satellites within the virial radius of either the MW or M 31 (∼300 kpc) are devoid of gas and quenched (Spekkens et al. 2014;Putman et al. 2021), with the most notable exceptions being the LMC and SMC.Initial analyses of the ELVES population (Carlsten et al. 2022) suggests that it follows a similar trend.However, in the case of SAGA, the vast majority of the satellites are still star-forming.Karunakaran et al. (2021) found that over 85% of the SAGA satellites are star-forming, in comparison to less than 40% of equivalently selected satellites from the A Project Of Simulating The Local Environment (APOS-TLE, Fattahi et al. 2016;Sawala et al. 2016) and Auriga (Grand et al. 2017) simulations.
The satellite mass ranges sampled by SAGA and ELVES are overlapping but different, complicating attempts at direct comparisons.Greene et al. (2023) find that the quenched fraction in ELVES systems is comparable to the LG, and Font et al. (2022) and Carlsten et al. (2022) argue that incompleteness in SAGA for low surface brightness (LSB) dwarfs could cause the apparent discrepancy between SAGA and simulations and the LG.However, Karunakaran et al. (2023) explored a number of SF metrics and selection biases in the SAGA and ELVES samples in an attempt to make like-for-like comparisons, concluding that the quenched fraction of satellites is similar for both (after accounting for differences in sample selection) and in tension with the high quenched fraction seen in the LG and simulations.
Quenching is a pivotal process that galaxies undergo when they fall into a group or cluster.Quenching transforms a galaxy from a member of the star-forming population typically found in the field, to a quiescent member of a larger structure (Dressler 1980;Grebel et al. 2003).Theory and simulations indicate that ram pressure stripping (RPS) likely places a significant role in driving satellite quenching all the way from Milky Way-mass groups to the largest clusters (e.g.Gunn & Gott 1972;Vollmer et al. 2001;Tonnesen & Bryan 2009;Simpson et al. 2018;Oman et al. 2021;Boselli et al. 2022;Wright et al. 2022).There are even some claims of relatively isolated, low-mass galaxies undergoing RPS in the field, presumably due to collisions with denserthan-average regions of the intergalactic medium (e.g.Benítez-Llambay et al. 2013, 2017;Beale et al. 2020;Yang et al. 2022).
RPS can occur when a gas-bearing galaxy moves at high speed through the intracluster medium (ICM) in a galaxy cluster or the circumgalactic medium (CGM) of the central galaxy in a galaxy group.If the ram pressure exerted on the gas content of the infalling galaxy exceeds the gravitational attraction between the gas and the other components of the disk, then gas will be stripped (Gunn & Gott 1972).If a galaxy loses all or most of its gas reservoir through RPS, then it will no longer be capable of sustained SF and will quench.
There are numerous examples of RPS in progress in galaxy clusters (e.g.Kenney & Koopmann 1999;Kenney et al. 2004;Koopmann & Kenney 2004;Chung et al. 2009;Poggianti et al. 2017;Boselli et al. 2018;Ramatsoku et al. 2019;Cortese et al. 2021;Jones et al. 2022;Boselli et al. 2023), which are often called "jellyfish" galaxies because of the tendrils of gas (and, frequently, in situ SF, e.g.Jáchym et al. 2014;Boselli et al. 2018) that extend from their disks.The phenomenon is also common in large groups (e.g.Roberts et al. 2021).However, for MW-mass systems there are only a handful of clear examples of ongoing RPS, such as, the LMC in the Local Group (LG) itself (e.g.Luks & Rohlfs 1992;de Boer et al. 1998), Holmberg II in the M 81 system (Bureau & Carignan 2002;Bernard et al. 2012), and ESO 324-G024 in the Cen A system (Fritz et al. 1997;Johnson et al. 2015).It has also been suggested that the LG Phoenix and Pegasus dwarfs might be experiencing RPS (Young et al. 2007;McConnachie et al. 2007).
In this work we approach the issue from a different perspective by trying to understand the details of the process of quenching in the SAGA systems, where most satellites have either yet to quench or are currently in the process of quenching.We have conducted observations of eight (plus one partial) SAGA systems, 42 (45) confirmed satellites, with both the Jansky Very Large Array (VLA) and the Canada-France-Hawaii Telescope (CFHT).The VLA H i observations are to map the gas reservoirs of the satellites and to look for signs that gas is being lost.We identify four cases of likely RPS in progress, effectively doubling the number of known examples in MW-like systems.The CFHT Hα observations trace recent SF, providing a uniform census of star formation rates (SFRs) and a robust means to separate quenched and star-forming satellites.By assessing both the SF, as well as the gas that fuels it, we aim to build a clearer picture of how quenching is proceeding in these systems.
This paper is organized as follows.In the following section we discuss the SAGA sample from which our targets were selected, as well as the VLA and CFHT observations.In §3 we present our results, including H i masses (and limits) and SFR estimates for all targeted satellites.We discuss these results in §4 and present our conclusions in §5.
Throughout this work we assume that all satellites are at the same distances as their SAGA hosts (unless stated otherwise) and we use host distance estimates from Mao et al. (2021), which are reproduced in Table 1.

SAMPLE & OBSERVATIONS
Our target sample was selected from the SAGA-II release (Mao et al. 2021); here we briefly describe the SAGA selection process.The SAGA project (Geha et al. 2017) used the HyperLEDA database (Paturel et al. 2003;Makarov et al. 2014) to select all MW-like host galaxies in the distance range 25 ≲ D/Mpc ≲ 40, based on the distance estimates of the NASA-Sloan Atlas (Blanton et al. 2011).SAGA defines "MW-like" as −23 > M K > −24.6, intended to approximately correspond to the MW's stellar mass, with the K-band magnitudes taken from the 2MASS Redshift Survey catalog (Huchra et al. 2012).The final list of hosts was trimmed to reject those with bright foreground stars and those that are themselves within the virial radius of a much larger system (Mao et al. 2021, Section 2.1.2).
SAGA then selected satellites around each host within 300 kpc (projected), taken as the approximate virial radius of a MW-like galaxy.Potential satellites were first selected from the photometric catalogs of SDSS DR14 (Abolfathi et al. 2018), the Dark Energy Survey DR1 (Abbott et al. 2018), and the DESI legacy imaging surveys DR6 and DR7 (Dey et al. 2019).These catalogs were combined and cleaned of duplicates, Galactic cirrus, and other spurious detections (Mao et al. 2021).The satellite candidates without prior redshift measurements were then targeted with follow-up spectroscopic observations.As not all candidates could be observed, SAGA target primarily those within a specific region of surface brightness-color space (Mao et al. 2021, equation 3), designed to capture all galaxies with z < 0.015.Finally, only the candidates within 275 km s −1 of the host galaxy redshift are considered as satellites and retained in the final SAGA-II catalog.
Starting with the SAGA-II catalog, for our H i and Hα observations we selected all systems with three or more satellites that are North of Dec = -25, and therefore straightforward to observe with the VLA and CFHT, a total of 17 systems and 85 satellites.We then excluded systems in the range 13.6 h < RA < 16.8 h as these were not observable in the CFHT 2021B semester during which our project began.This resulted in an initial sample of 10 hosts and 49 satellites.

VLA observations
These 10 SAGA systems were observed between July and September 2022 as part of the VLA project 22A-023 (PI: M. Jones).Despite the fact that the SAGA systems nominally extend over approximately 1 sq deg, in practice the satellites are usually arranged such that all can be imaged with the VLA in just a few 30 ′ -wide L-band primary beams.Of the 10 systems observed, the maximum number of separate pointings required was three, with the average being two.In general these observations were lightly impacted by radio frequency interference (RFI) and antenna outages, however, NGC 2906, NGC 4454, and UGC 4906 were all more severely affected.For NGC 2906 (a single pointing) no usable data were collected, while for NGC 4454 and UGC 4906 about half of the data were lost to RFI.In addition, due to an error in the scheduling block files, the second pointing for NGC 4454 was not observed, meaning that 3/4 satellites in that system have either poor or no coverage.This reduced the sample to 8 systems (42 satellites) with complete H i coverage and one additional system with partial coverage.
The data were reduced using CASA (McMullin et al. 2007) and the H i pipeline 1 of Jones et al. (2023).Each observation was reduced separately following standard automated and manual flagging and calibration proceedures.Multiple pointings were combined into mosaics for imaging.All fields were imaged with Briggs robust values of 0.5 and 2. The former represents a compromise between sensitivity and spatial resolution, while the latter optimizes sensitivity for sources that are comparable in size to the beam (or larger).Multi-scale CLEANing and auto-masking were used within CASA's tclean task.The images were CLEANed down to approximately 2.5σ within the mask.Spectral resolution was also re-gridded to 5 km s −1 during imaging.The final RMS noise values and synthesized beam sizes for each system are shown in Table 1.Mao et al. (2021) and assuming a mass-to-light ratio of 0.96 and M⊙,K = 3.41, as in (Longhetti & Saracco 2009).Distance estimates are from Mao et al. (2021).The RMS noise and beam size values quoted here are (at the pointing center) for the robust=2 H i cubes.There is no H i mass measurement for PGC 13646 as it was outside of the VLA primary beam.† Note that two pointings are required to image all of the NGC 4454 satellites in SAGA, but only one pointing was observed.
Source masking in the image cubes was performed using SoFiA v1.3.2 (Serra et al. 2015).The standard smooth and clip algorithm was implemented with no spatial smoothing, as well as smoothing at roughly half and equal to the synthesized beam diameter.In addition, the data were smoothed to a spectral resolution of ∼15 km s −1 .A threshold of 4σ was used to clip the smoothed data cube and the candidate detections were then merged.A reliability threshold of 95% was enforced to remove spurious detections.These source masks were then used to generate moment maps and spectra of each detected target using the primary beam-corrected H i data cube.
In the case of non-detections, we returned to the (primary beam-corrected) data cube and extracted spectra over an area equal in size to one synthesized beam (from the robust=2 cube), centered at the optical center of each target.The RMS noise in each of the spectra was measured within ±100 km s −1 of the known velocity of each target.These noise values were then used to estimate 3σ upper limits on the H i masses of each non-detection, assuming a nominal velocity width of 30 km s −1 and a roughly top-hat shaped spectral profile.There is one satellite (LS-429812-2469 of UGC 903) that we marginally detected.However, we still use the upper limit as the peak S/N is only ∼2.
For the detected satellites H i masses were calculated simply by integrating all the flux within their source mask generated by SoFiA.As the H i morphology of a number of the satellites is quite disturbed, and because they mostly have relatively narrow spectral profiles (not exhibiting double-horn features), this simplistic approach was preferred to profile fitting.However, in three cases manual adjustments were made to the SoFiA masks.The H i emission from two satellites were lightly blended with their hosts, while NSA-542307 was blended with a low-level RFI feature that we were unable to remove via flagging.These were all manually separated using the 3-dimensional visualization tool SlicerAstro (Punzo et al. 2016;Punzo et al. 2017) in order to measure their individual H i masses.
We use the same approaches as above to measure the host H i masses and limits in Table 1.However, for the undetected hosts we assumed a fiducial velocity width of 100 km/s.Three of the hosts (UGC 903, NGC 4158, and NGC 7541) are detected in the Arecibo Legacy Fast ALFA (Arecibo L-band Feed Array) survey, or ALFALFA (Giovanelli et al. 2005;Haynes et al. 2018), and after correcting for the different distance estimates all three of our H i mass measurements agree within 1σ with those from ALFALFA.

CFHT observations
All SAGA systems observed with the VLA, except NGC 4454, were also observed in Hα and r-band with the MegaCam instrument on CFHT during programs 21BC05 and 23AC01 (PI: K. Spekkens).The MegaCam imager consists of 40 9.5 megapixel CCDs, covering a full square degree field of view (FoV).This meant that every target system could be observed in its entirety in a single pointing, although in some cases, where the satellite distribution was lopsided, the host was not placed at the center of the FoV.Each target was observed with either a five-or seven-point large dither pattern to fill in gaps between the CCDs.The nominal exposure times for each system were 2000 s in Hα and 300 s in r-band, however, in some cases usable exposures were taken during partial executions of observing blocks, leading to slightly longer total integrations.These images were generally taken during grey time with seeing below 1.2 ′′ .The images were processed with the Elixir-LSB pipeline (Ferrarese et al. 2012) and mosaicked with MegaPipe (Gwyn 2008).Continuum-subtracted Hα images were produced by subtracting the r-band images from the Hα images.
Hα fluxes were measured using the Aperture Photometry Tool (Laher et al. 2012).The Hα morphology of most of the detected satellites is highly irregular and clumpy, thus elliptical apertures (with annuli for sky subtraction) were manually adjusted for each source.In cases of Hα non-detections, a circular aperture of 5.61 ′′ (30 pixels) in radius, and centered on the satellite coordinates from Mao et al. (2021), was used to estimate a 3σ upper limit for the Hα flux.To convert the observed magnitudes, m Hα , to Hα fluxes we used the expression: F Hα = 3.621 × 10 −20 ∆ν10 −0.4mHα erg s −1 cm −2 , where ∆ν is the filter bandwidth in Hz (7.18×10 12 Hz).Fluxes were corrected for Galactic extinction (Fitzpatrick 1999;Schlafly & Finkbeiner 2011), but not internal extinction, and then converted to luminosities (in erg s −1 ) based on the host distance estimates (Table 1) and finally to SFRs following Kennicutt (1998), i.e.SFR Hα = 7.9 × 10 −42 (L Hα /erg s −1 ) M ⊙ yr −1 .We assumed an uncertainty in this conversion of 0.2 dex, as in Kennicutt (1998), which was added in quadrature with the flux measurement uncertainties.

Gas fraction of satellites
The H i mass measurements and limits from §2.1 (Table 2) were combined with the stellar mass estimates of Mao et al. (2021), to calculate gas fractions, M HI /M * , for every satellite in the nine systems targeted (with the exception of LS-321038-4238 in the NGC 4454 system, which was not observed with the VLA).These gas fractions are plotted as a function of projected separation between host and satellite in Figure 1.Detections are plotted as filled blue circles, while upper limits are unfilled circles.The LG satellites from Putman et al. (2021) are shown as filled and unfilled grey squares for comparison (stellar mass estimates from Carlsten et al. 2022).In both cases the symbol sizes are scaled with log M * .A floor for gas fraction upper limits is set at log M HI /M * = −2.0 in order to prevent the vertical axis becoming too stretched.
Focusing first on the detections, we find a pronounced trend of decreasing gas fraction with decreasing hostsatellite projected separation, declining approximately 1 dex over 300 kpc.This is the qualitative expectation if the satellites are being stripped of their gas by interac-tions with the central galaxy and its hot gas halo, however, the situation appears notably distinct from that of the LG where virtually all satellites within 300 kpc are devoid of gas.This is discussed further in §4.1.Related trends in the H i content of satellite galaxies in groups have been found in various works (e.g.Hess & Wilcots 2013;Odekon et al. 2016;Brown et al. 2017;Jones et al. 2020;Zhu & Putman 2023), but to our knowledge this is the first such measurement for a sample focused on satellites around MW-like hosts.
The upper limits are mostly mixed in with the detections and are not ruling out the possibility that the non-detections could be following the same trend as the detections.The exceptions are the three limits that are significantly below all of the detections.These are the mostly likely cases where satellites have lost their H i reservoirs.These are all fairly massive (log M * /M ⊙ ∼ 9) satellites with low specific star formation rates, that have likely been members of their current systems for some time, and as a result have lost the majority of their initial gas reservoir.The other non-detections are likely a mixture of genuinely gas-poor satellites and those with H i masses too low to be detected simply because these satellite are low mass (both stellar and H i).

Star formation rates
SFRs for most SAGA-II satellites have been estimated previously based on GALEX NUV (and some FUV) fluxes (Karunakaran et al. 2021).Hα traces SFR over the past ∼10 Myr, as opposed to SFR over a timescale on the order of ∼100 Myr with NUV (e.g. Lee et al. 2011), thus the two tracers are highly complementary.Here we present Hα-based SFR estimates from the CFHT imaging and compare to the previous NUV-based estimates.
Figure 2 shows the comparison between our SFR estimates and those based on NUV from Karunakaran et al. (2021). 2 In general we see excellent agreement between the two measurements.There is a slight tendency for points at high SFRs to be below the one-to-one line (i.e.SFR Hα < SFR NUV ) and those at low SFRs to be above it.There is one satellite detected in Hα, but undetected in NUV, and with the Hα SFR measurement in tension with the SFR limit from NUV.This is probably a result of the very shallow GALEX observations for this particular object causing the previous upper limit to be in error.We also note that optical spectra from Mao et al. (2021) detect Hα line emission in this object.There is another outlier (LS-595052-1940)  LGS3 LG SAGA Figure 1.Gas fraction of SAGA satellites as a function of projected distance from their hosts (blue points and limits).Values for Local Group satellites are shown in grey for comparison (Putman et al. 2021).In this case the separation is not the projected separation, but the smaller of the physical separation between the satellite and either the MW or M 31.Filled symbols indicate H i detections and unfilled symbols are upper limits.The symbol diameters are scaled by log M * .A floor in gas fraction upper limits is set at log MHI/M * = −2.0.Error bars are omitted for clarity.We remind the reader that many Local Group satellites would be too faint/low-mass to be included in SAGA, hence many LG satellites have symbols that are too small to see (log M * /M⊙ ≲ 6).
Comparison between NUV-and Hα-derived SFR estimates for the SAGA satellites in common between this work and Karunakaran et al. (2021).The solid black line is the line of equality and the dotted black line is a fit to the satellites detected in both NUV and Hα.
above of the main locus of points.In this case we suspect that a background UV source in close proximity to the dwarf led to an overestimate of SFR NUV .Using HyperFit (Robotham & Obreschkow 2015) we fit a straight line to all the satellites with detections in Figure 3. Specific star formation rate of SAGA satellites as a function of projected distance from their hosts.Filled blue points indicate H i detections, while unfilled points show nondetections (only for objects where the H i mass limits are less than the stellar mass estimates, i.e.only those guaranteed to be gas-poor).Symbols are scaled by log M * , as in Figure 1.
both Hα and NUV (excluding the outlier LS-595052-1940).This fit is shown by the dashed black line in Figure 2. As expected the best fit slope is slightly below unity and there is a non-zero intercept.The gradient and intercept values are 0.84 ± 0.06 and −0.36 ± 0.09, respectively.
Figure 3 shows the specific star formation rates (SS-FRs) of the SAGA satellites as a function of (projected) host-satellite separation.There is a drop-off in SSFR for H i-detected satellites that are particularly close to their host (D proj ≲ 100 kpc).For the H i non-detections we only include those for which the upper limit provides a tight constraint on the gas fraction (log M HI /M * < 0), as the other non-detections could be just as gas-rich as the detections (cf. Figure 1).These non-detections all have low SSFRs, indicating that proximity to the host is not the only relevant factor here.We will return to this point in §4.2.

Ram pressure stripping in progress
Upon visual inspection of the H i maps of the detected SAGA satellites it became apparent that four objects display signs of possible ongoing RPS in their H i morphology (Figure 4).As ram pressure exerts a force on the gas in a galaxy, but not its stars, a telltale (necessary, but not sufficient) sign of RPS is a gas tail containing no stars (or only stars that likely formed recently in situ).In addition, RPS tails should occur in only one direction, as a galaxy's gas is swept back into its wake as it moves through the CGM or ICM.In contrast, tidal tails typically form both leading and trailing arms as the tidal field is roughly symmetric about the satellite.Tides are also gravitational and therefore impact all matter equally.Thus, if stars are coincident with gas at the stripping sites, then stellar tails should accompany the gas tails.However, in any real scenario gas stripping is unlikely to be driven purely by tides or ram pressure and they can be difficult to reliably separate in practice.Below we will argue that there are a number of factors that indicate that these four cases (Figure 4) are very likely dominated by RPS, but ultimately the confirmation of this will require further, more detailed observations.
In the cases of NSA-19694 and LS-359104-131 (Figure 4, left panels) we see clear H i tails with no apparent accompanying stellar component.For LS-359104-131, the tail extends only to the SW, with no convincing signs of extensions in any other direction (given the sensitivity and resolution of the data).For NSA-19694 there is a main tail that extends to the NE as well as a small nub on the NW side of the galaxy.Although we do not know the 3-dimensional velocity of the galaxy, and therefore cannot know the direction of any ram pressure exerted on it (we return to this point below), this configuration would make sense if ram pressure is forcing gas to the NW and the two tail components come from opposite sides of the galactic disk.The gas distribution of galax-ies undergoing RPS can often form a bow shape as the outskirts of the disk (on either side) are stripped before the center (where the gravitational restoring force is stronger).This bow shape can be slightly asymmetric either because the original H i distribution is asymmetric (which is exceedingly common, e.g.Espada et al. 2011) and/or because the orientation of the galaxy results in one side of the disk taking the brunt of the ram pressure interaction.
Returning to LS-359104-131, in Figure 5 we show a wider field moment zero map made with the robust=2 H i data cube ( §2.1), which has better column density sensitivity for extended emission (but poorer angular resolution).To produce this map we also repeated the masking process with SoFiA using a lower threshold (3.5σ), accompanied by an even higher reliability criterion (>99%), in an attempt to include low significance H i emission surrounding brighter emission without adding spurious detections.
In Figure 5 we see an H i cloud with no apparent optical counterpart approximately 13 ′ (∼140 kpc in projection) to the south, on the other side of the host (NGC 7541) and its largest satellite (NSA-637123).Although it first appears that this H i is likely from an ongoing interaction between NGC 7541 and NSA-637123, inspection of their H i morphologies and kinematics indicates that neither of these galaxies are significantly disturbed, and have likely not yet begun strongly interacting, despite their small angular and velocity separations.Thus, it is unlikely that the orphaned gas cloud originated from either of these galaxies and we instead speculate that it is an extension of the RPS tail of LS-359104-131.Indeed, we note that its velocity is only ∼60 km s −1 offset from LS-359104-131 and there appears to be low significance indications that LS-359104-131 and the southern H i cloud may connect via a spatially (and kinematically) continuous stream of gas.While a detailed investigation into the nature of this potential connection is beyond the scope of this work, we share the results of our preliminary investigation via an interactive Jupyter notebook hosted on Binder3 .We leverage the 3D rendering capability of IPyvolume and a masked subcube of the NGC7541 field to highlight the low S/N connection between LS-359104-131 and the southern cloud.Deeper H i imaging might reveal a continuous structure at higher S/N or perhaps the stripped gas does not have a high enough density to self-shield and there may be gaps, especially near NGC 7541, where a portion could  already have been photoionized.We note that with a full end-to-end length of ∼15 ′ (∼170 kpc, in projection) and an assumed maximum velocity of 275 km s −1 , if this is one continuous structure then the southern portion must be at least ∼600 Myr old.NSA-343657 and DES-373393030 (Figure 4, right panels) show a slightly different morphology, but are still consistent with RPS.In these cases the VLA data do not have sufficient resolution to clearly resolve tails, but the H i distributions are offset in one direction from the optical centers of these galaxies.DES-373393030 is almost entirely unresolved, but its H i distribution is clearly centered to the NW.If we assume that it has not yet reached pericenter and is still falling towards its host, NGC 14214 , then this would mean that the bulk of its H i gas is trailing in its wake.Finally, for NSA-343657 the H i distribution forms almost a "V" shape, which we suspect is another example of the bow shape typical of RPS, but smeared out by the large and elongated beam shape.However, it should also be noted that there is a slight suggestion of a leading tail, which could make the distribution more of a "Z" shape, which would instead suggest a tidal interaction as the cause.The resolution of the current data notwithstanding, the fact that most of the galaxy's gas appears to be residing in its wake (again assuming it is falling towards its host), still means this object is a strong candidate for RPS in progress.
The white arrows in Figure 4 point in the direction of the host galaxy of each satellite.All four RPS candidates are significantly less than 100 kpc (projected) from their hosts, comparable to the separation between the MW and LMC.However, we do not know if they are pre-or post-pericenter passage in their orbits.Before pericenter passage, a satellite's RPS tail should point roughly away from the host, whereas shortly after pericenter passage it would point towards the host.In all four cases we see that the tails are roughly (anti-)aligned with the direction of the host, strengthening the notion that they are indeed the result of RPS.

A new satellite of UGC 903
In the UGC 903 system we detected an additional satellite that was not in the SAGA-II sample, which we will refer to as dw0122+1735 (Figure 6).This is a LSB dwarf approximately 2.7 ′ to the SE of UGC 903 that we identified via its H i line emission, which is lightly blended with that of UGC 903, but clearly concentrated on the LSB optical counterpart (Figure 6, right).
We used GALFIT (Peng et al. 2002(Peng et al. , 2010) ) and DECaLS images to fit the optical properties of this dwarf galaxy.As it is faint and LSB we used a g + r co-add image to fit the structural properties, obtaining a half-light radius R eff = 10.4±0.3 ′′ , a Sérsic index n = 0.36±0.03,and an axial ratio a/b = 0.66 ± 0.02.We then refit in each band separately to obtain the photometry, holding the parameters above fixed.As the object is slightly irregular, and not perfectly modeled with a Sérsic profile, we also perform manual aperture photometry and use the difference between the two as an estimate of the uncertainty.This gives g = 19.4±0.4,µ 0,g = 25.9±0.4,r = 18.8±0.3,and µ 0,r = 25.2 ± 0.3 (extinction corrected following Schlafly & Finkbeiner 2011).For consistency with other objects in Table 2 we use the stellar mass estimator from Mao et al. (2021), which gives log M * /M ⊙ = 7.56.We separate the H i emission from UGC 903 using SlicerAstro to measure an H i mass of log M HI /M ⊙ = 7.53 ± 0.09.These correspond to a gas fraction of 0.93.We take the flux-weighted average velocity of the H i profile to estimate the radial velocity of dw0122+1735 as 2439 km s −1 (∆cz = −77 km s −1 relative to UGC 903).
Finally, we note that based on the photometry above and the distance estimate to the UGC 903 system (Table 1), dw0122+1735 falls well within the standard crite- ria (van Dokkum et al. 2015) to be considered an ultradiffuse galaxy (UDG), with an effective radius of 1.9 kpc and M V = −13.9.  2 both indicate that unlike satellites in the LG, most SAGA satellites have significant H i gas reservoirs.At least part of the explanation for this difference in gas fractions is the different distribution of stellar masses between SAGA and the LG (cf.Karunakaran et al. 2020a).SAGA is only complete for relatively massive satellites (log M * /M ⊙ ≳ 7.5), while the MW only contains a total of three satellites in that stellar mass range (and meeting the other SAGA selection criteria), while M 31 has six (e.g.McConnachie 2012).This means that SAGA is only seeing the massive end of the satellite population.For example, if we take the SMC stellar mass as approximately 10 8.5 M ⊙ (e.g.Besla 2015) then there are 14 satellites more massive than the SMC in the nine target SAGA systems, which make up about half of the H i detections (three are non-detections), while the remaining half almost all have stellar masses above 10 7.5 M ⊙ .The gas fractions of both the LMC and SMC appear to be roughly consistent with the observed trend in Figure 1, and it is this class of satellites (or slightly lower mass) that we are generally detecting with our H i observations.We also note that although there are no direct equivalents in SAGA to the most the stringent gas fraction limits for LG satellites plotted in Figure 1, this is simply a function of sensitivity and distance.Likely some of the H i non-detections from the SAGA sample are just as gas-poor as many of the LG satellites (of similar stellar mass), but given the distance to the SAGA systems it is not possible to obtain comparable limits on their H i content.
However, despite these differences, meaningful comparisons can be made, especially to simulations that have attempted to reproduce the situation in the LG.For example, Simpson et al. (2018) use the Auriga simulation suite (Grand et al. 2017) to assess quenching in MW-like systems.They find that ram pressure stripping is the dominant mechanism removing gas from satellites and that this acts rapidly (<1 Gyr) to quench them.At a stellar mass of log M * /M ⊙ ≈ 8 they predict that half of satellites (within 300 kpc of host) should be gaspoor.Looking at the satellites in our sample that are 7.75 < log M * /M ⊙ < 8.25, we indeed see that five are gas-poor and five are gas-rich.We consider a gas fraction of unity to be the dividing line and err on the side of over-counting gas-poor satellites when there are only upper limits for the H i masses.We note that this is slightly different to the definition used by Simpson et al. (2018), which was a fixed threshold of M HI = 10 5 M ⊙ at all stellar masses.
At slightly higher stellar masses we see that most satellites have significant gas reservoirs.This means that either ram pressure stripping is no longer effective at removing their gas, or they are on their first infall (but smaller satellites are not).Comparing again to Simpson et al. (2018), we see that this finding is also consistent.As shown in Simpson et al. (2018) figure 4, log M * /M ⊙ = 8 is roughly the threshold at which ram pressure stripping becomes inefficient, with almost no satellites at log M * /M ⊙ = 9 being quenched.
These simple comparisons suggest that, despite the first impression that most SAGA satellites are starforming and most LG satellites are quenched, the fraction and masses of the satellites detected in H i is roughly consistent with the expectations from the Auriga simulations, which were designed to reproduce MW analog systems.However, in Figure 1 we also noted a clear trend between the observed gas fraction and project host-satellite separation.This trend is presumably a signature of the gas removal mechanism.Given that there is expected to be considerable host-to-host scatter in CGM density (e.g.Ramesh et al. 2023;Saeedzadeh et al. 2023) it is unclear if ram pressure stripping should produce such a uniform trend when data are compiled from several different systems.We will revisit this point in more detail §4.5.

Star formation rates
The fit between Hα and NUV-derived SFR estimates in Figure 2 shows that although the two metrics mostly agree, there is a slight tendency for Hα to underpredict the SFR relative to NUV at low SFRs and the reverse at high SFRs.A similar trend was first noted for galaxies with low SF activity (SFR ≲ 0.1 M ⊙ yr −1 ) in Lee et al. (2009).Even though that work used FUV (which is more directly comparable to Hα) rather than NUV, the slope that we fitted is consistent within the uncertainties.Lee et al. (2009) also performed a thorough accounting of all conceivable sources of this discrepancy, but concluded that no one factor could explain it, with the possible exception of a varying stellar initial mass function.We have few detections in the low SFR regime and further investigation of this offset is beyond the scope of this paper.
When the SSFRs of the SAGA satellites were investigated in Figure 3, there were no objects with clearly elevated SF.In particular, we note that none of the four satellites that appear to be currently undergoing RPS have elevated SSFRs, nor do they show any obvious signs of fronts of SF.This suggests that ram pressure experienced by these satellites is insufficient to produce the compression-induced waves of SF sometimes seen in higher mass galaxies undergoing RPS (e.g.Gavazzi et al. 1995;Vulcani et al. 2018;Boselli et al. 2021).This find-ing is consistent with other work that has found that this type of SFR-enhancement is ineffective in galaxy groups (e.g.Vulcani et al. 2021;Roberts et al. 2021), presumably because of the lower density of the intergalactic medium and the lower velocities involved.
In Figure 3 we also saw that there is a downturn in the SSFR of the H i-detected satellites at small (projected) host-satellite separation (D proj ≲ 100 kpc).In Figure 7 we show that this is likely a result of the suppressed gas fractions at small projected separations.At approximately a gas fraction of unity there is a sharp drop in SSFR of around 0.5 dex.Despite H i not being the direct fuel for SF, this transition is so marked that there is only one SAGA object with a gas fraction below unity and log SSFR/yr −1 > −10.5.
To verify that this transition is a feature specific to satellites, we drew a random selection of ALFALFA galaxies from the Durbala et al. (2020) catalog that is 10 times the size of our SAGA sample, but matched in terms of its distribution of stellar masses.5 Unsurprisingly we see that the ALFALFA sample extends to higher SFR and gas fraction values.This is a predominantly a field galaxy population (e.g.Guo et al. 2017;Jones et al. 2020) and its intrinsic H i-selection favors gas-rich galaxies.However, crucially we see a smooth extension of the main relation between gas fraction and SSFR below a log M HI /M * = 1 where there are no SAGA satellite galaxies, suggesting that this disconnect is not merely a feature of the relationship, but rather that it is specific to satellites.The sharpness of this transition also suggests that it is the underlying cause of the trend seen in Figure 3, rather than there being a direct correspondence between host-satellite separation and SSFR, e.g. as a result of strong tidal interactions.
In Figure 7 we also overlay two literature measurements of the relationship between gas fraction and SSFR.In green we show the relation from the extended GALEX Arecibo SDSS sample (xGASS, Catinella et al. 2018), which is another field sample but is not H iselected (unlike ALFALFA).This relation roughly overlaps the ALFALFA points at low gas fractions, but if extrapolated higher would be slightly to the left of the ALFALFA population (a result of the H i-selection of ALFALFA).The orange relation is from the spectral line stacking work of Brown et al. (2017) and shows the average relation for satellites in groups with halo masses in the range 12 < log M halo /M ⊙ < 13.This relation has a slightly lower SSFR at a gas fraction of unity than the xGASS relation, and then follows a much steeper decline in SSFR with decreasing gas fraction.This relation is a remarkable match to the trend in our H i detections and again suggests that we are seeing behavior that is specific to satellites.Brown et al. (2017) argued that their results indicate that the timescale for gas removal from satellites was considerably shorter than that for the subsequent shutdown of SF.Several other works have also made similar suggestions (e.g.Wetzel et al. 2013;Oman et al. 2021).However, if this were the case then we would expect to see a population of H i-poor satellites with elevated SSFRs (for their gas fractions).Instead we see the opposite, a population of satellites with low SSFRs for their gas fractions, relative to field galaxies (i.e.ALFALFA and xGASS).This suggests that the shutdown of SF does operate on a similar timescale to the gas removal mechanism and that satellites are not largely stripped of their gas before SF starts to shutdown, at least not in the mass regime applicable to SAGA.
Finally, in Figure 7 we include the H i values for the ELVES sample from Karunakaran et al. (2022).ELVES (Carlsten et al. 2022) has fewer MW-like hosts than SAGA, but they are significantly closer (D < 12 Mpc instead of 25 < D/Mpc < 40 for SAGA) and thus ELVES is complete to much lower satellite masses.This sample has fewer detections than our SAGA sample and it is difficult to draw conclusions from so few values.However, curiously they do not appear to follow the same trend, with the low gas fraction points mostly lying between the two comparison relations.It may be that the current small samples from both SAGA and ELVES paint a somewhat misleading picture that may not be borne out by future observations.Ultimately, more data are needed to draw robust conclusions from this analysis.
In Appendix A we also explore the trend between gas fraction and SSFR in simulations.However, the comparison to SAGA is complicated by their offset gas fractions (see §4.5).

An overlooked, un-quenched satellite
The newly identified satellite of UGC 903 (dw0122+1735) presented in §3.4 was not included in the SAGA sample.In this section we briefly discuss the significance and the likely reason for this omission.
The (projected) separation between dw0122+1735 and UGC 903 is 2.7 ′ (30 kpc).This means that it was not eliminated as a result of the 10 kpc exclusion radius that SAGA enforced in the vicinity of the hosts.At r = 18.8 it is also significantly brighter than the nominal magnitude limit of SAGA (r = 20.75,Mao et al. 2021).Thus, the LSB of this object is the most likely cause of it being missed during the SAGA search for satellite candidates.There is also considerable Galactic cirrus in the vicinity of UGC 903, which may have contributed to it being overlooked.
There has been some debate in the literature (Karunakaran et al. 2021(Karunakaran et al. , 2023;;Font et al. 2022;Carlsten et al. 2022) about whether or not SAGA is systematically missing LSB satellites that would otherwise meet their nominal selection criteria.This satellite, dw0122+1735, is likely an example of exactly that, but one lone example does not settle this debate, which is a statistical consideration.However, it does raise an important complication.Font et al. (2022) and Carlsten et al. ( 2022) both argued that SAGA under-counts the number of quenched satellites because (they argue) quenched dwarfs tend to have lower surface brightness.But, the possibility of missing un-quenched galaxies (due to the same hypothesized observational bias) was not considered.Although dw0122+1735 is not presently forming stars (it is undetected in UV and Hα), it does contain an appreciable neutral gas reservoir, which could fuel SF in the future, and therefore cannot be considered quenched.If other equivalent objects have been missed then it is even possible that the quenched fraction of satellites calculated from SAGA is underestimated and could disagree even more severely with simulations than indicated in Karunakaran et al. (2021).
Finally, as noted in §3.4,dw0122+1735 meets the usual surface brightness and physical size criteria to be considered a UDG.Gas-bearing UDGs are not uncommon (e.g.Leisman et al. 2017;Jones et al. 2018;Janowiecki et al. 2019;Karunakaran et al. 2020b), however, this particular UDG is noteworthy as it is a satellite, gas-rich, and relatively red (g − r = 0.6), the combination of which makes it somewhat unusual.There is one similar object in the sample of Karunakaran et al. (2020b), but it is also an order of magnitude more massive.

Ram pressure stripping toy model
The H i observations of the four satellites apparently in the process of being ram pressure stripped (Figure 4) are not at sufficiently high resolution to allow us to accurately measure the truncation of the leading edge of their H i disks.The truncation of the H i disk is a direct result of the impact of ram pressure and in the case of the LMC has been used to estimate the CGM density of the MW (Salem et al. 2015).
Although we cannot use the detailed H i properties of the satellites to estimate the CGM densities of the SAGA hosts, we can use their combined global H i properties if we make some simplifying assumptions and approximations.To do this we construct a toy model inspired by that of Wang et al. (2020).First, we assume that all satellites were previously field dwarfs, and when they originally fell into their current system, had H i content typical of field dwarfs of similar stellar mass.We also assume that all subsequent H i mass loss is the result of RPS from passing through the CGM of their present day hosts.Starting with these assumptions we can construct a toy model to estimate the CGM density as described below.
In what follows we will only consider the satellites that were detected in our H i observations.This is because for these satellites it is reasonable to assume that they recently fell into the group and/or are still undergoing RPS.For the satellites devoid of gas, they could potentially have been satellites for many gigayears and their present day gas-poor state might be the result (at least in part) of effects other than RPS.For example, they may have been partially ram pressured stripped and then consumed the remainder of their gas through SF.
To approximate the ram pressure stripping threshold we will use the common Gunn & Gott expression: where ρ cgm is the CGM mass density, v is the velocity of the satellite through the CGM, Σ * is the stellar surface density, Σ gas is the gas surface density, R is the galactocentric radius, and G is the gravitational constant.We note that this form of the expression implicitly assumes a 100% self-gravitational coupling for the gas in the disk.
For the stellar surface density we simply assume an exponential disk with total mass estimated by Mao et al. (2021) for each satellite.The scale radius is estimated via the stellar size-baryonic mass relation of Bradford et al. (2015).Although we could estimate this more accurately by re-fitting the light profiles of each satellite, in practice Σ * is not particularly important as the H i distribution is typically much more extended than the stellar distribution and the Σ 2 gas (R) term dominates over the Σ * (R)Σ gas (R) term.Thus, this simple approximation is sufficient for our toy model.
To estimate the original (at infall) H i mass of each satellite we use the stellar mass-gas mass (where M gas = 1.4MHI ) scaling relation of Bradford et al. (2015) as this focuses on field galaxies in the appropriate stellar mass range.In addition to knowing the total original gas mass, we must also know how it was originally distributed in order to evaluate Equation 1. Fortunately, the radial distribution of H i in galaxies is extraordinarily similar (after normalization of the scale radius and amplitude) over a wide range of masses (Wang et al. 2014(Wang et al. , 2020)).We use the averaged profile from Wang et al. (2016) as the template profile for the original radial distribution of H i in all of the SAGA satellites.The H i radii are calculated using the Wang et al. (2016) H i size-mass relation.
The final piece of information that we need is to know how fast the satellites are traveling through the CGM.In practice the relevant quantity is the maximum value of ρ cgm v 2 (which presumably occurs at pericenter) as this corresponds to the smallest value of R where the gas distribution will be stripped.The only constraint that we have on the value of v at the orbital pericenter is a lower limit from the present day value of ∆cz, the difference in the line-of-sight velocities of the host and satellite (Table 2).For the upper limit we choose 275 km s −1 as this is the limit where the SAGA selection criteria would no longer consider a satellite candidate to be bound to the host.
With a means to estimate all of the relevant quantities we follow a Monte Carlo approach to estimate ρ cgm .The following steps are repeated 1000 times for each satellite in our sample.1.The initial H i mass M HI,init is drawn from the stellar mass-gas mass relation (with scatter), with the requirement that M HI,init ≥ M HI,0 , where the 0 subscript denotes the value today.
2. The H i scale radius is drawn from the H i sizemass relation (with scatter) based on M HI,init .
3. The stellar distribution scale radius is drawn from the stellar size-baryonic mass radius (with scatter) based on M * + 1.4M HI,init .
4. The satellite maximum velocity relative to the CGM is drawn from a uniform distribution over the range |∆cz| ≤ 275 km s −1 .
5. Based on the amount of H i mass that is missing, M HI,init − M HI,0 , the RPS truncation radius of the gas distribution is calculated, along with the corresponding CGM density.
The median values of ρ cgm , with 1σ error bars, are plotted for each of the satellites detected in H i in Figure 8.The horizontal blue dash-dot line is the average value for all satellites and hosts (log ρ cgm /g cm −3 ≈ −27.3).For comparison we consider an estimate of the MW CGM density (Figure 8, horizontal black dashed line).Salem et al. (2015) calculated the peak ram pressure experienced by the LMC based on the truncation of its H i disk along its leading edge, which in combination with its orbital parameters provided an estimate of the MW CGM density at the pericenter of the LMC's orbit (∼50 kpc).The average CGM density value that we calculate is slightly higher than that of Salem et al. (2015), which might indicate that we are overestimating somewhat as it is unlikely that most of the SAGA satellites have passed closer to their host than 50 kpc (although most of their orbital parameters are unknown).However, considering the long list of significant assumptions that we needed to make above our estimate is still relatively close to that of the MW.
Although there are many assumptions, the step that results in the majority of the uncertainty seems to be step 4, randomly drawing the maximum velocity of the satellite through the CGM.Both the mass surface density and the velocity are squared quantities in Equation 1, however, owing to the uniformity of most H i disks, the latter is is much more uncertain.The satellites with the smallest error bars are those with the largest values of |∆cz|, i.e. the tightest constraints on the velocity.Similarly, the average value of ρ cgm is highly sensitive to the upper limit for the velocity that we set above.
For example, increasing this value to ∼350 km s −1 would result in agreement with the range of values for ρ cgm for the MW.Therefore, we cannot draw strong conclusions from this toy model, but we note that the similarity to the MW values is remarkable given the number of assumptions.Very tight constraints via this approach would require detailed knowledge of the orbital parameters of the satellites, which unfortunately will not be feasible for most SAGA systems.However, a considerably larger sample than the one presented here would statistically sample the distribution of orbital parameters, and could be matched with theoretical priors for those parameters in order to construct a tighter constraint.In addition, a larger sample is likely to contain a significant number of satellites with large values of |∆cz|, which would further strengthen the constraint.
In the future we plan to obtain higher-resolution H i imaging of the four satellites apparently undergoing RPS such that we can explicitly measure their truncation radii, thereby eliminating the need for most of the steps above and placing a direct constraint on the maximum ram pressure that they have experienced.The corresponding maximum velocity will still present a challenge for estimating the corresponding CGM density, however, with only four satellites we will attempt to construct tighter constraints on these velocities based on the directions of their ram pressure tails, their current hostsatellite separations, and the values of ∆cz.

Gas fractions and projected separations in simulations
Over the past few years there have been several investigations of the quenching of satellites in MW-like systems in cosmological and zoom-in simulations (e.g. Simpson et al. 2018;Digby et al. 2019;Akins et al. 2021;Wright et al. 2021;Engler et al. 2023).Although there is not a consensus on exactly how satellites quench in simulations, RPS is noted as a physical process that is likely important.Engler et al. (2023) analyzed the satellite populations in MW and M 31-like systems in the TNG50 simulation finding broad agreement between the properties of these simulated satellites and genuine LG satellites.They also measured a decline in H i mass as a function of host-satellite separation for relatively massive (M * > 10 8.5 M ⊙ ) satellites in TNG50, analogous to that shown in Figure 1.They also find clear examples of ongoing RPS, with similar morphologies to the observations in Figure 4. Similar RPS tails have also been noted in the DC Justice League simulations (Akins et al. 2021).Engler et al. (2023) find that only a few percent of satellites should host visible RPS tails at z = 0.As we have identified four strong candidates out of 45 satellites, it seems that the occurrence rate might be a factor of a few times higher.However, they consider satellites down to M * = 5 × 10 6 M ⊙ , and the timescale for gas loss and quenching is expected to scale inversely with satellite mass.
As noted by Crain et al. (2017) for the Evolution and Assembly of GaLaxies and their Environments (EA-GLE) simulations, and as we will discuss further below, it is common for current simulations to systematically under-predict the H i masses of low-mass galaxies.However, as pointed out by Wright et al. (2021) in their study of satellite quenching in EAGLE, when dealing with relative gas masses (e.g. between quenched and star-forming galaxies), or changes in gas mass, this bias is of minimal importance as the physical processes that are removing the gas are likely quite distinct from those driving this bias.In order words, it is the transition that is of most interest here and the simulations can still provide insight even if the absolute H i content is offset.
To make a more direct comparison between the observed properties of satellites and theoretical expectations, we make use of the APOSTLE (Fattahi et al. 2016;Sawala et al. 2016) and TNG50 (Pillepich et al. 2019;Nelson et al. 2019a,b) simulations.APOSTLE is dedicated to simulating systems that are comparable to the MW and LG, and is thus well-matched to our sample.It consists of 12 simulation volumes analogous to the MW and M 31 system, or a total of 24 MW-like satellite systems.For simplicity we use only the intermediate resolution L2 run of APOSTLE and require satellites to have a stellar mass of log M * /M ⊙ > 7.1 to roughly mimic the SAGA selection criteria and to ensure that each consists of at least 100 star particles in the simulation (which have masses of 1.1 × 10 5 M ⊙ ).This gives a sample of 114 satellites.
As TNG50 is not focused on MW-like systems we had to define our own set of hosts, taking all those in the stellar mass range 10.5 < log M * /M ⊙ < 11.5 (cf.Table 1).Potential hosts with a more massive neighbor within 300 kpc were removed, as well as those that were themselves satellites within a massive halo (log M total /M ⊙ > 13).Finally, as TNG50 contains some objects that have significant gas mass, but almost no dark matter, we place a lower limit of 10 8 M ⊙ on the sub-halo mass of satellites to avoid erroneously including any such objects as satellites.This gives a total of 246 TNG50 hosts6 and 2030 satellites (with log M * /M ⊙ > 7.1).With this mass cut each satellite stellar component is made up of at least 100 particles (of mass 8.5 × 10 4 M ⊙ ).Satellite stellar masses are taken as the sum of all stellar mass particles that are gravitationally bound to the sub-halo.We estimate H i masses for gas particles in APOSTLE following the convention for simulations using the EAGLE galaxy formation model established by Bahé et al. (2016, see also Crain et al. 2017 for further discussion), combining the partitioning scheme for neutral-to-ionized gas of Rahmati et al. (2013) and the atomic-to-molecular partitioning scheme of Blitz & Rosolowsky (2006).For TNG50, we instead adopt the H i profiles of Diemer et al. (2019, see also Diemer et al. 2018).We use their "projected" (or "map", i.e. not "volumetric") profiles -see Stevens et al. (2019) for a discussion motivating this choice -and the Gnedin & Draine (2014) atomic-to-molecular partitioning scheme.
For both simulations we extracted the total H i mass within 1, 5, 10, 30, and 50 kpc of the center of mass of each satellite.The H i mass of each galaxy was then taken as the value within the first step that exceeds four times the stellar half mass radius.This was compared to instead selecting the step where the H i surface density dropped below 1 M ⊙ /pc 2 , and was found to be mostly consistent.However, the surface density approach gave spurious results for some galaxies (likely those with central holes in their H i distributions), and thus the simpler approach was favored for being more resilient to outlying cases.
Figure 9 shows the comparison of the gas fractions of the H i-detected SAGA satellites (blue points) and the APOSTLE (top) and TNG50 (bottom) simulations as the background intensity maps.The maps for the simulations were generated by projecting each system on to 10,000 different, randomly chosen lines of sight.Additional selection criteria were enforced to match those from SAGA (Mao et al. 2021).With each different projection only satellites with a projected separation from their host of 10 < D proj /kpc < 300 and a velocity separation |∆cz| < 275 km s −1 were retained.We also required the H i masses of satellites to be above 107.5 M ⊙ in order to roughly approximate our H i detection limit (which in practice varies somewhat between different satellites and systems). 7 Both simulations show a qualitatively similar behavior to the observations of the SAGA sample, with gas fractions declining with projected separation.In the case of TNG50 the distribution of gas fractions appears constant to D proj ≈ 100, after which a significant downturn begins.For APOSTLE there is a gradual change in gas fraction between D proj = 300 kpc and 0 kpc (but the interpretation suffers from the small sample size compared to TNG50).The SAGA sample with H i measurements is also currently quite small and the data could reasonably be fit with either of the functional forms discussed above, a plateau with a downturn around 100 kpc, or a roughly linear decline across the entire range of D proj values.
A clear difference between the observations of the SAGA satellites and those in both simulations is that on average the observations lie at higher gas fractions (by approximately 0.5 dex).Even at D proj = 300 kpc the simulations do not reach the values seen in the observations, which implies that this discrepancy must extend into the field.As noted above, it is a known issue (see also e.g.Oman et al. 2019, figure 1) that simulations currently tend to under-predict the H i content of lowmass galaxies.Thus, the satellites in the simulations likely never had as much H i as the SAGA satellites did when they first fell into their present day systems.
Despite this offset, the existence of analogous trends in gas fraction in our observations and the APOSTLE and TNG50 simulations, as well as the reproduction of many other satellites properties (e.g.Sawala et al. 2016;Engler et al. 2023), indicates that we can use the simulations to investigate the processes that are likely driving the observed trends.For many of the intermediate mass (6 ≲ log M * /M ⊙ ≲ 8) satellites in SAGA the dominant factor is likely RPS (cf.Digby et al. 2019).However, at higher masses there is likely a transition (e.g.Fillingham et al. 2019) where the satellite's potential well becomes deep enough that RPS is no longer an effective means of quenching.Even if these high mass satellites are observed to have RPS tails, RPS likely cannot remove enough gas to quench their SF.Instead starvation could be the dominant mechanism that slowly shuts off SF, leading to much longer quenching timescales.
In future work we hope to expand the sample to permit more detailed comparisons with simulations, such that we can begin to understand where this transition occurs.In particular, with a larger sample (both from simulations and observations) we could attempt to identify the radius at which the suppression of gas fraction and SFR begins for satellites of a fixed stellar mass.This would indicate the point where ram pressure begins to have a meaningful impact on the satellites.If ram pressure is the dominant mechanism then this characteristic radius should get progressively smaller for higher mass satellites, with the quenching mechanism eventually transitioning to starvation (as outlined above).Such findings would be strong evidence for ram pressure as the driving force for satellite quenching, with a transition to starvation at higher masses.

CONCLUSIONS
We have fully imaged eight SAGA systems (42 satellites) in both H i and Hα, providing tracers of both the ongoing star formation and the gas that fuels it.A ninth system (3 satellites) has been partially observed.We observe a clear decline in gas fraction (M HI /M * ) as a function of projected host-satellite separation.Over the full range of separations considered (∼25-300 kpc) the typical gas fraction changes by around an order of magnitude.A similar trend is seen in our comparisons to the APOSTLE and TNG50 simulations, however, there remains a ∼0.5 dex offset in gas fraction for all satellites, preventing a full like-for-like comparison.We also find that the fractions of gas-rich and gas-poor satellites are consistent with those from the Auriga simulations, although a larger sample is needed for detailed comparison.These simulations predict that RPS is the dominant mechanism driving gas removal, and it is likely that we are indirectly seeing its effects through this trend in gas fraction.
Four satellites also exhibit signs of likely ongoing RPS, in the form of clear one-sided H i gas tails (Figure 4).All four are within 75 kpc (projected) of their host galaxy and are likely close to the pericenter of their orbits, where the CGM density and ram pressure are expected to be highest.In one case the tail appears to extend back around to the opposite side of the host (Figure 5).Only a handful of examples of satellites undergoing RPS in MW-mass systems are known, and thus these four represent a significant expansion of that sample.Furthermore, identifying four examples in only eight systems suggests this is a remarkably common phenomenon and that further observations of MW-like systems have a strong chance of uncovering additional examples.We plan to obtain higher-resolution H i imaging in order to confirm and further explore the morphology of these features.
We constructed a RPS toy model to attempt to estimate the 3-dimensional mass density of the CGM based on the current and assumed initial gas content of the satellites.This approach ( §4.4) contains some considerable assumptions and uncertainties, but permits an order of magnitude estimate of the CGM density.We find a value about 0.5 dex higher than estimated for the MW (Salem et al. 2015).The constraining power of any individual galaxy depends strongly on how well aligned its orbital velocity vector is with the line-of-sight.We hope to place better constraints on the CGM density in these systems with an expanded sample as well as higher resolution observations that may permit the H i disk truncation to be spatially resolved.
In the UGC 903 system our H i observations revealed a previously unknown satellite that was not included in the SAGA sample, but does meet the nominal selection criteria.This satellite was likely overlooked because it has a very low surface brightness and, in fact, fits the standard definition of a UDG.As it is gasbearing, it cannot be considered quenched as it may go on to form more stars in the future, but it is undetected in Hα and UV and is therefore not presently forming stars.This finding complicates an active debate about whether SAGA is systematically missing LSB satellites and whether this would lead to a suppression of the quenched fraction observed (Karunakaran et al. 2021(Karunakaran et al. , 2023;;Font et al. 2022;Carlsten et al. 2022).
Our Hα imaging produced SFR estimates that were largely consistent with the NUV estimates from Karunakaran et al. (2021).Thus, in this case the main advantage of Hα over UV as a SF tracer is the angular resolution and we will present an assessment of the Hα morphologies of the satellites in a subsequent paper.In Figure 2 we tentatively identified a trend analogous to that of Lee et al. (2009), where at very low SFRs there is a tendency for Hα to slightly underpredict the SFR relative to UV.However, we currently have few detections in this regime.
In line with the trend in gas fractions, there is also a suppression of SSFR at small projected separations.This appears to be driven by the lower gas fractions at small separations.There is a sharp decline in SSFR for satellites with gas fractions below unity.Although there is currently only a small number of satellites in this regime, they are consistent with the Brown et al. (2017) scaling relation (from H i stacking) between gas fraction and SSFR for satellite galaxies.This suggests that the interplay between gas content and SF may be fundamentally different for the SAGA galaxies in comparison to field galaxies, by virtue of the fact that they are satellites rather than centrals.This behavior does not seem to be reproduced in the simulations that we compared to, but their offset gas fractions make this comparison problematic (see Appendix A).Karunakaran et al. (2021) sparked a continuing debate about the high fraction of SAGA satellites that are starforming versus quenched, and whether this is inconsistent with our understanding of the satellites in LG and models that are designed to match them.Based on our observations of eight SAGA systems we find no significant disagreement in our Hα-based SFR estimates and the UV-based estimates of Karunakaran et al. (2021).However, we also do not find significant deviations in the expected fraction of gas-rich and gas-poor satellites from the Auriga simulations (Simpson et al. 2018), and the general behavior of satellite gas fractions is roughly consistent with the APOSTLE and TNG50 simulations.
Our findings above indicate that H i observations in particular offer a unique insight into the ongoing processes in these systems.Ultimately, we need to expand the available observations of gas in these systems in order to enable more detailed comparisons and to draw more robust conclusions.
We thank the anonymous referee for their thorough reading of the paper and constructive suggestions.et al. 2016;Punzo et al. 2017), GALFIT (Peng et al. 2002(Peng et al. , 2010)), dustmaps (Green 2018), matplotlib (Hunter 2007), numpy (van der Walt et al. 2011), scipy (Oliphant 2007;Millman & Aivazis 2011), pandas (Wes McKinney 2010; pandas development team 2020), HyperFit (Robotham & Obreschkow 2015), IPyvolume, spectral-cube, SoFiA (Serra et al. 2014(Serra et al. , 2015)).7 showed the relationship between gas fraction and SSFR for the SAGA satellites and various comparison samples.The sharp downturn in SSFR around a gas fraction of unity suggests that this relationship may be different for field and satellite galaxies.Here we make the same plot (Figure 10), but for satellites from the TNG50 and APOSTLE simulations.Satellites were selected as those within 300 kpc (physical) and ∆v < |275| km s −1 of their host.Satellites within 10 kpc of their host were removed, in line the with SAGA selection.We also removed any satellites with log M * /M ⊙ < 7.1 to make the sample more directly comparable to SAGA and to ensure all the simulated satellites are well resolved.
Here we see that, unlike the SAGA satellites, the simulated satellites mostly follow the field galaxy (xGASS) relation between gas fraction and SSFR.As expected (cf.§4.5) the simulated satellites are generally more gas-poor than those in the observations.While at first glance it might appear as if the SAGA observations are mostly only sampling the gas-rich tail of the population seen in the simulations, the stellar mass cut mentioned above means that both the observational and simulated samples correspond to a mass regime where SAGA is relatively complete.The H i observations are naturally more likely to detect gas-rich satellites, but even taking the non-detections into account, a significant fraction of the SAGA sample fall at the very extreme of the gas-fraction-SSFR parameter space sampled by the simulations.Thus, what we are mostly seeing here is that the simulations are under-predicting the H i masses, by around 0.5 dex.However, if the simulations merely increased the gas masses of these satellites, then they would likely continue to follow the green relation on which they currently fall, and would therefore probably end up with SSFRs that are higher than the SAGA satellites.This finding indicates that the resolution to this issue is not as straightforward as simply making the simulated satellites more gas-rich.It may be that the relationship between 12 h 29 m 00 s 28 m 50 s 40 s

Figure 4 .Figure 5 .
Figure 4. VLA H i moment zero (robust=0.5)contours overlaid on DECaLS grz images for the four examples of likely ram pressure stripping in progress.From top-left to bottom-right the satellites (hosts) are NSA-19694 (NGC 4454), NSA-343657 (NGC 6278), LS-359104-131 (NGC 7541), and DES-373393030 (NGC 1421).In each case a white arrow points in the direction of the host galaxy.The ellipses in the lower-right corner of each panel indicates the VLA synthesized beam size.The contours begin at 3σ and each subsequent contour is double the previous.For top-left to bottom-right these lowest contour values are 8.4 × 10 18 , 4.1 × 10 18 , 5.7 × 10 18 , and 6.9 × 10 18 cm −2 over 10 km s −1 .

Figure 6 .
Figure 6.Left: DECaLS grz image of dw0122+1735.Right: DECaLS grz image of UGC 903 and dw0122+1735 (left of center) with H i moment zero (robust=2) contours overlaid.The contours begin at 3σ and each subsequent contour is double the previous.The semi-transparent ellipse in the lower-right corner indicates the synthesized beam size.There is a clear concentration of H i centered on the optical source enlargerd in the left panel.
Figure1and Table2both indicate that unlike satellites in the LG, most SAGA satellites have significant H i gas reservoirs.At least part of the explanation for this difference in gas fractions is the different distribution of stellar masses between SAGA and the LG (cf.Karunakaran et al. 2020a).SAGA is only complete for relatively massive satellites (log M * /M ⊙ ≳ 7.5), while the MW only contains a total of three satellites in that stellar mass range (and meeting the other SAGA selection criteria), while M 31 has six (e.g.McConnachie 2012).This means that SAGA is only seeing the massive end of the satellite population.For example, if we take the SMC stellar mass as approximately 10 8.5 M ⊙ (e.g.Besla 2015) then there are 14 satellites more massive than the SMC in the nine target SAGA systems, which make up about half of the H i detections (three are non-detections), while the remaining half almost all have stellar masses above 10 7.5 M ⊙ .The gas fractions of both the LMC and SMC appear to be roughly consistent with the observed trend in Figure1, and it is this class of satellites (or slightly lower mass) that we are generally detecting with our H i observations.

Figure 7 .
Figure 7. Specific star formation rate of SAGA satellites plotted as a function of their gas fractions.H i detections are shown as filled points and non-detections as unfilled points.The size of the points is scaled by log M * , with examples shown on the right side of the plot.ELVES satellites with H i and NUV SFR measurements (Carlsten et al. 2022; Karunakaran et al. 2022) are included for comparison.Grey error bars in the background show a stellar mass-matched comparison field sample from ALFALFA.The orange line shows the Brown et al. (2017) relation from satellites in groups of 12 < log M halo /M⊙ < 13, and the green line shows the xGASS relation for field galaxies in the stellar mass range 9 < log M * /M⊙ < 10 (Catinella et al. 2018).The vertical black dashed line indicates a gas fraction of unity.

Figure 8 .
Figure8.Estimates of CGM density from the satellites in each SAGA system that are detected in H i (blue error bars).The blue horizontal dot-dash line indicates the weighted average of all the points (log ρcgm/g cm −3 ≈ −27.3) and the black dashed line (and dotted lines) indicate the CGM density estimate for the MW fromSalem et al. (2015).

Figure 9 .
Figure 9.Comparison of observed gas fractions of SAGA satellites detected in H i (blue circles) as a function of projected host-satellite separation and the greyscale maps showing the corresponding values averaged over 10,000 random orientations of MW-like satellite systems in the APOSTLE (top) and TNG50 (bottom) simulations.

Figure 10 .
Figure 10.Specific star formation rates of satellites plotted as a function of their gas fractions.Equivalent to Figure 7 showing the APSOTLE and TNG50 simulations in comparison to the SAGA sample observations.The size of the points is scaled log M * , with examples shown on the right side of the plot.The orange line show the Brown et al. (2017) relation from satellites in groups of 12 < log M halo /M⊙ < 13, and the green line shows the xGASS relation for field galaxies in the stellar mass range 9 < log M * /M⊙ < 10 (Catinella et al. 2018).The vertical black dashed line indicates a gas fraction of unity.

Table 1 .
SAGA systems targeted This work used images from the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID 2014B-0404; PIs: David Schlegel and Arjun Dey).Full acknowledgment at https://www.legacysurvey.org/acknowledgment/.This work used data from the Karl G. Jansky Very Large Array.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.The data were observed as part of program 22A-023 (PI: M. Jones).DJS acknowledges support from NSF grant AST-2205863.KS acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC).AK acknowledges financial support from the grant (SEV-2017-0709) funded by MCIN/AEI/ 10.13039/501100011033 and from the grant POSTDOC 21 00845 funded by the Economic Transformation, Industry, Knowledge and Universities Council of the Regional Government of Andalusia.DC acknowledges support from NSF grant AST-1814208.

Table 2 .
Satellite H i and star formation properties