The Satellites Around Galactic Analogs (SAGA) Survey, first presented in Geha et al. (2017, hereafter Paper I) , is an ongoing spectroscopic galaxy survey. Its primary science goal is to characterize the satellite populations down to the luminosity of the Leo I dwarf galaxy (M r,o < −12.3), around 100 Milky Way (MW) analogs outside of the Local Volume (at 25-40.75 Mpc). SAGA Paper I included 27 satellites identified around 8 MW-like systems. In this work, we present the Stage II results of the SAGA Survey, which includes 127 satellites identified around 36 MW-like systems and marks the completion of more than one-third of the full survey. * NASA Einstein Fellow; yymao.astro@gmail.com The SAGA Survey complements the rich observational data sets of MW satellite galaxies (Drlica-Wagner et al. 2020, and references therein) , and provides critical tests of the Lambda-Cold Dark Matter (ΛCDM) paradigm. Apparent discrepancies between the observed MW satellite population and ΛCDM predictions have been considered as "small-scale challenges" to ΛCDM (e.g., Weinberg et al. 2015; Bullock & Boylan-Kolchin 2017) . These challenges include ΛCDM over-predicting both the number and the central mass densities of the MW satellites; however, it has long been unclear if these discrepancies can be explained via a more realistic treatment of baryons (e.g., Brooks et al. 2017; Read & Erkal 2019) , if the ΛCDM paradigm itself is incorrect (e.g., Polisensky & Ricotti 2014) , or if other salient simplifications have impacted the comparison. These small-scale challenges have been largely based on the MW and M31 (e.g., Tollerud et al. 2014) , just two systems with possibly correlated satel-lite populations (Garrison-Kimmel et al. 2014) . The SAGA Survey addresses these limitations by constraining the distribution of satellites around MW-mass galaxies, putting the MW and its bright satellite population into a cosmological context and allowing us to compare observations and theory statistically.
The SAGA Survey also provides a benchmark for understanding low-mass galaxy evolution. The galaxy-halo connection is reasonably well constrained for galaxies whose stellar mass M * > 10 9 -10 10 M , but much uncertainty still exists below this scale (e.g., Somerville & Davé 2015; Wechsler & Tinker 2018) . Unresolved questions include how the interaction between satellite galaxies and their hosts affects satellite star formation histories (e.g., Simpson et al. 2018; Hausammann et al. 2019) . The SAGA Survey focuses on bright and classical dwarf satellite galaxies (M * ∼ 10 6 -10 9 M ), providing crucial information for our understanding of the galaxy-halo connection and host-satellite interactions. For example, a surprising result in Paper I was that the vast majority (26/27) of satellites in 8 SAGA systems were star-forming, as compared to 2 of 5 MW satellites in the same luminosity range.
In recent years, efforts have been undertaken to characterize satellite systems beyond the Local Group, but around nearby systems in the Local Volume (< 20 Mpc). Within this distance, one can use resolved stars or surface brightness, in addition to spectroscopic redshifts, to search for satellites. MW-like satellite systems that have been studied in the Local Volume include NGC 4258 (Spencer et al. 2014 ), M94 (Smercina et al. 2018) , NGC 3175 (Kondapally et al. 2018) , NGC 2950 and NGC 3245 (Tanaka et al. 2018) . A handful of satellite systems around galaxies of different masses have also been studied, such as the systems around M81 (Chiboucas et al. 2013) , M101 (Danieli et al. 2017; Bennet et al. 2019; Carlsten et al. 2019) , and Centaurus A Müller et al. 2019) . Regardless of the host galaxy mass, the number of complete satellite systems in the Local Volume remains limited. It is not until very recently that a statistical study of the Local Volume satellite systems became possible Carlsten et al. 2020b,c) .
Beyond 20 Mpc, dwarf galaxies are difficult to distinguish from the far more numerous background galaxy population using photometry alone (although attempts have been made, e.g., Xi et al. 2018) . Confirming satellite galaxies hence requires extensive spectroscopic follow-up (e.g., Zaritsky et al. 1997) , which is the main challenge of the SAGA Survey. Existing deep spectroscopic surveys either explicitly remove low-redshift galaxies (e.g., DEEP2; Newman et al. 2013) , or have limiting magnitudes which allow study of only the brightest satellites (e.g., SDSS, r < 17.77; Strauss et al. 2002; Sales et al. 2013; and GAMA, r < 19.8; Baldry et al. 2018 ). The SAGA Survey was hence born to tackle this challenge before ongoing, planned, and proposed deep spectroscopic surveys become available (e.g., DESI; DESI Collaboration et al. 2016; 4MOST; de Jong et al. 2012; MSE; McConnachie et al. 2016) .
On the theoretical side, modern hydrodynamical simulations span a wide range of scales, including cosmological volumes that include large samples of MW-mass galaxies but only resolve their most massive satellites (e.g., Schaye et al. 2015; Pillepich et al. 2018) , zoom-in simulations of individual MW analogs that resolve many of their satellites (e.g., Wetzel et al. 2016; Garrison-Kimmel et al. 2019a) , and highresolution simulations that resolve individual star forming regions around individual dwarf galaxies (e.g., Wheeler et al. 2019; Munshi et al. 2019) . SAGA results will serve as observational tests against the physical models implemented in such simulations (e.g., Samuel et al. 2020; Akins et al. 2020) .
Theoretical frameworks that characterize how host halo properties affect satellite populations (Mao et al. 2015; Fielder et al. 2019) , along with empirical and semi-analytical models to produce realistic satellite systems (Nadler et al. 2019b; Jiang et al. 2020) , are also rapidly being developed. These theories and galaxy-halo connection models enable statistical inferences that can marginalize over the uncertainties in baryonic physics, providing a new way to test the ΛCDM paradigm. Nadler et al. (2019a Nadler et al. ( , 2020a have tested this approach with MW satellites to constrain the faint end of the galaxy-halo connection and alternative dark matter models. Incorporating SAGA results in this type of analyses will undoubtedly provide more powerful constraints by adding data in a regime that has been sparse to date.
In this paper, we present the Stage II results of the SAGA Survey, including 127 satellites in 36 MW-like systems. We first review our survey design, including host selection and survey progress ( §2), photometric catalog production ( §3), target selection and spectroscopic instruments ( §4), and the definition of satellite and survey completeness ( §5). We then present our findings, including the properties and quenching fraction of individual satellites ( §6), the luminosity functions, radial distributions, and co-rotating signals of the satellite systems ( §7), and the comparison between our observation and theoretical predictions ( §8). Readers who are familiar with Paper I or want to quickly navigate the main results of this paper can read the summary ( §9) first.
As in Paper I, all distance-dependent parameters are calculated assuming H 0 = 70 km s −1 Mpc −1 and Ω M = 0.27. Magnitudes and colors are extinction corrected (denoted with a subscript 'o,' e.g., r o ; using a combination of Schlegel et al. 1998 and Schlafly & Finkbeiner 2011; see §3) . Absolute magnitudes are K-corrected to z = 0 using Chilingarian et al. (2010 
Galactic Plane Imaging not available SDSS imaging LS DR6/7 imaging DES DR1 imaging New Hosts in This work (28) Paper I Hosts (8) Figure 1 . The sky distribution of 205 Milky Way-analog systems identified in the SAGA Survey volume (circles and squares). There are 115 SAGA systems considered for targeting that have complete imaging over at least 99% of the projected virial radius in SDSS (yellow), Legacy Surveys (LS) DR6+7 (green), and/or Dark Energy Survey (DES) DR1 (blue). The solid black line indicates the Galactic plane. There are 36 hosts with complete spectroscopic coverage: 8 from Paper I (orange squares) and 28 newly presented in this work (red squares).
The primary science goal of the SAGA Survey 1 is to characterize the satellite galaxy populations of 100 Milky Way analogs down to an absolute magnitude of M r,o = −12.3. In this section we describe the survey at a high level: first we describe our criteria to select MW analogs (Section 2.1), then outline the SAGA Survey strategy and timeline to achieving the project goals (Section 2.2).
The selection of our Milky Way analogs ("hosts") is primarily based on K-band luminosity and local environment criteria. Our science questions are best answered by comparing the satellite populations of hosts with similar dark matter halo masses, a property that is impossible to directly measure at this mass scale. Among direct observables, K-band luminosity and local environment are the best available proxies for halo masses given available data.
As described in Paper I, we select MW hosts from a complete nearby galaxy catalog ('Grand List' hereafter). This 'Grand List' was constructed when we began the survey in 2014, mainly based on the HyperLEDA database (Makarov et al. 2014) 2 and supplemented by various literature sources. Since then, the HyperLEDA database has incorporated additional sources and properties from the astronomical literature. We therefore rebuild our 'Grand List' to take advantage 1 sagasurvey.org 2 leda.univ-lyon1.fr of these new additions, while maintaining consistency with our prior results.
We aim to construct our 'Grand List' so that it is complete out to 60 Mpc in distance and down to K = 11.9 (corresponding to M K = −22 at 60 Mpc). We first obtain all objects in the HyperLEDA database that are galaxies (objtype = G), have both heliocentric velocity (v) and distance (modbest), and have a K-band magnitude less than 12 mag (kt < 12). Each object in the HyperLEDA has a unique PGC ID that serves as the main identifier.
To maintain consistency with our prior results, we opt to use the distance measurement from the NASA-Sloan Atlas (NSA; Blanton et al. 2011 ) and the total K-band magnitude from the 2MASS Redshift Survey catalog (2MRS; Huchra et al. 2012 ). More specifically, we spatially match the objects with the NSA v1.0.1 catalog, 3 and choose the NSA distance (ZDIST) when available, unless HyperLEDA has a redshift-independent distance measurement and the final error (e_modbest) is less than 0.04 dex, in which case we choose the HyperLEDA distance measurement. The distance reported in the NSA is corrected for local peculiar velocity with the Willick et al. (1997) model. While this is not necessarily the best distance measurement for each galaxy, the availability of quality distance measurements in our target distance range (between 25 and 40 Mpc) is scarce, and the available distance estimates usually have an uncertainty of 0.1 dex (Leroy et al. 2019 We define "low-z" as 0.003 < z < 0.03; "discovered" indicates that the redshift was first obtained by SAGA.
For K-band magnitude, we use the Extragalactic Distance Database (EDD; Tully et al. 2009 ) 4 to obtain extinctioncorrected, total Kand H-band magnitudes from the 2MRS catalog. We use the prescription in Chilingarian et al. (2010) and K − H color to K-correct the K-band absolute magnitude to z = 0. For objects that do not appear in the 2MRS catalog, we use the K-band magnitude in HyperLEDA but shift it by −0.03 mag. The value of this small shift is derived from the median difference for objects in both catalogs and does not have significant impact on our host list construction.
With the compiled information in our 'Grand List,' we define our Milky Way hosts based on four criteria: (1) an absolute K-band criteria, (2) a set of cuts to reduce stellar foregrounds, (3) cuts based on the local galactic environment and (4) distance cuts for our survey. We detail these criteria below and summarize in Equation 1.
First, we require our Milky Way hosts to have an absolute K-band magnitude in the range −23 > M K > −24.6. We use this magnitude as a proxy for total stellar mass. As shown in Paper I, this range corresponds to the possible magnitude of the central galaxy hosted by a halo of virial mass of 1.6 × 10 12 M , when a simple galaxy-halo connection model is assumed (see Figure 2 of Paper I). The corresponding stellar mass range for the host galaxy is 10 10 -10 11 M . The virial radius corresponding to this halo mass is ∼ 300 kpc.
We next require that hosts are sufficiently far from the Galactic plane to reduce stellar foregrounds (|b| ≥ 25 • ). To avoid saturated regions around very bright stars, we search the Hipparcos-2 catalog (van Leeuwen 2007) within the virial radius. We remove hosts with a very bright foreground star 4 edd.ifa.hawaii.edu H P < 5 (H P is the broad optical Hipparcos magnitude) which can significantly reduce the region of usable photometry.
We then examine the candidate host's local galactic environment. Using our 'Grand List,' we first require there are no bright galaxies inside the projected virial radius with magnitude of K < K host −1.6 or brighter. 5 The magnitude difference of 1.6 mag was chosen to ensure that no two MW host galaxies in our K-band magnitude range could be found inside the same virial volume. This criteria does not exclude MW-M31 like systems, as the Milky Way is slightly beyond two virial radii from M31 (such a configuration is only excluded for ∼ 35% of the lines-of-sight solid angle). We also use the 2MRS group catalog constructed by Lim et al. (2017, obtained from EDD) and remove hosts associated with a group whose mass is greater than M halo < 10 13 M . Finally, we manually flag a small handful of potential host galaxies that are visibly disturbed and/or in merging systems.
We select galaxies in the distance range of 25-40.75 Mpc with a heliocentric velocity greater than 1,400 km s −1 . The outer distance limit corresponds to the absolute magnitude limit of M r,o < −12.3 given our apparent magnitude limit of r o < 20.75. The inner distance and velocity limits exclude 80 hosts in fields that are too big for our instruments and that are more efficiently studied using other techniques (e.g., Danieli et al. 2018) . These nearby (< 25 Mpc) hosts that also satisfy our stellar mass and environments cuts include the wellstudied Milky Way analogs M94 (Smercina et al. 2018) , as well as 3 of the 10 hosts (NGC 1023 , NGC 2903 presented in Carlsten et al. (2020a) . 
The distance and environmental cuts above are slightly more stringent than those used in Paper I but include all previously published hosts. From our 'Grand List', we identify 205 MW analogs that pass these selection criteria. Of these hosts, 115 pass the available imaging requirements described in Section 3.8 to be considered for targeting. We plot the distribution of these hosts on the sky in Figure 1 . For each host, the SAGA Survey searches for satellites within a projected radius of 300 kpc, approximately corresponding to the virial radius of a MW-mass halo. Given the distance range of SAGA hosts (25-40.75 Mpc, or 0.003 < z < 0.013 in redshift), the search radius of 300 kpc corresponds to 25 -42 on the sky.
Among the 115 hosts that meet all our selection criteria, we identify 4 Local-Group-like pairs (two MW-mass hosts within 1 Mpc of each other). One of these pairs is present in this work: NGC 7166 and ESO 288-025 (separated by 668 kpc). We additionally identify 11 hosts that are in a Local-Group-like environment but whose companions do not meet our selection criteria or do not have imaging available; 3 of these hosts are presented in this work: NGC 5602, NGC 5750, NGC 5962. We plan to observe more fields around these Local-Group-like hosts in SAGA Stage III to enable studies of the correlations between satellite systems and the host environment.
The SAGA Survey is a long-term project that we consider in three main stages summarized in Table 1 : (I) Determine the complete luminosity function for a handful of MW analogs using basic criteria; (II) Design targeting strategies which more efficiently identify high probability satellites without compromising completeness; and (III) Apply these criteria to 100 MW-analog systems. In Paper I we presented the results of SAGA Survey Stage I: 8 complete systems with relaxed target selection criteria, over 12,682 newly obtained galaxy spectra around the 8 systems, 16 newly discovered satellites, and 99 newly discovered very low-redshift (z < 0.03) galaxies (see also Table 1 ).
In this work we present results from Stage II of the SAGA Survey. This includes two major improvements over Stage I. First, we significantly improved our target selection strategies. Based on data from SAGA Stage I, we reduced the average number of targets per host by more than a factor of four (from ∼ 1,500 to ∼ 400). We detail these target selection strategies in Section 4.1. Second, we expanded our photometric data to include both the Dark Energy Survey (DES) and the Legacy Surveys (LS) in addition to SDSS (Stage I used SDSS alone). This expansion provides us access to many more hosts, especially in the southern hemisphere (Figure 1) . In expanding our photometric coverage, some hosts have photometric coverage across multiple surveys. These data provide excellent test cases to understand differences among the photometric catalogs (see also Section 3.7). Table 1 provides key numbers regarding the progress of the SAGA Survey. The three column groups represent the three SAGA Survey stages. We show the numbers for two different field definitions in sub-columns of different fonts: within 300 kpc in bold font (corresponding to the virial radius), and within 1 • (corresponding to the field of view) in regular font around the SAGA hosts. The Stage II results presented here are based on 36 complete systems in total, reaching just over one-third of the goal of 100 complete systems (see Section 5.2 for the definition of a "complete system"). During Stage II, we completed 28 hosts with only 11,726 new redshifts (on average ∼ 430 redshifts within 300 kpc around each host). In SAGA Stage III, we will continue to improve our targeting efficiency, especially in the low-surfacebrightness regime, and aim to complete the remaining hosts with an average of ∼ 300 redshifts within 300 kpc around each host.
In the 36 complete systems of SAGA Stage II, we have identified 123 satellites down to M r,o = −12.3 (and 4 more below the threshold). The SAGA Survey provides the first redshift measurements to 66 of these satellites (54%). Among the fainter half of the satellites (62 satellites with −15.77 < M r,o < −12.3), 53/62 (87%) have the first redshift measurement from SAGA. The SAGA Survey has also discovered more than 200 additional very low-redshift (z < 0.03) galaxies around the 36 SAGA hosts.
As of July 2020, we have obtained close to 40,000 galaxy redshifts in total across all SAGA hosts. We are on track to complete spectroscopic follow-up for our Stage III goal by the end of 2021.
We next detail our pipeline to build the photometric object catalog for each SAGA host. Our photometric object catalogs provide the base catalogs for our target selection. In SAGA Stage I, we used only the SDSS DR12 catalogs as our photometric object catalogs. In SAGA Stage II, we expand to include deeper data sets from the Dark Energy Survey (DES) and Legacy Surveys (LS). 6 Our target selection concentrates on galaxies brighter than an extinction-corrected r o < 20.75. This is two or more magnitudes brighter than the DES or LS photometric magnitude limits.
We first describe the source catalogs from the three surveys (SDSS DR14, DES DR1, LS DR6/7; Section 3.1-3.3). We then describe the process used to merge these into a single photometric object catalog for each host (Section 3.4), the cleaning process (Section 3.5), and the validation of our final object catalog (Section 3.7). Our photometric object catalogbuilding pipeline is developed in an open-source repository. 7
We use the SDSS DR14 photometry catalog (Abolfathi et al. 2018 ) and ingest all photometric objects around each host within a one-degree radius. We use modelMag for galaxy magnitudes and PETROR50_R as the effective photometric radius R r,eff in our source catalog. To correct for extinction, we use the value reported in EXTINCTION, which is derived from Schlafly & Finkbeiner (2011) for DR14. We use the PHOTPTYPE flag for star/galaxy separation and select only galaxy targets.
We then apply a set of selection criteria to identify "good" galaxy targets:
where σ g,r,i are the reported one-sigma photometric errors in each band. The last two sets of criteria aim to remove objects that are anomalously bright due to sky over-subtraction near bright stars.
We use the DES DR1 photometry catalog (Abbott et al. 2018 ) and ingest all photometric objects around each host within a one-degree radius. We use MAG_AUTO_ * _DERED for de-reddened magnitudes (derived from Schlegel et al. 1998), 8 and use FLUX_RADIUS as our effective photometric radius R r,eff . We convert DES griz magnitudes to approximate SDSS magnitudes using the following formulas: 8 g SDSS = g DES − 0.0009 + 0.055(g − i) DES , 6 darkenergysurvey.org and legacysurvey.org 7 github.com/sagasurvey/saga 8 des.ncsa.illinois.edu/releases/dr1/dr1-faq
For star-galaxy separation, we use EXTENDED_COADD parameter in the r band, which is defined as: 8
and we require galaxies to have EXTENDED_COADD = 3. However, for very bright objects (r o < 17), the DES pipeline sometimes mis-classifies stars as galaxies. We hence applied additional cuts to bright objects to remove stars. In particular, we relabel objects that satisfy the following criteria as stars: 0.7 · (r o + 10.2) > r o + 2.5 log 2πR 2 ro,eff , g o − r o < 0.6, and
Finally, we apply the following selection to identify "good" galaxy targets: IMAFLAGS_ISO_R = 0, and FLAGS_R < 4.
We use the combined Legacy Surveys DR6 (BASS/MzLS) and DR7 (DECaLS) photometry catalogs (Dey et al. 2019) and ingest all photometric objects around each host within a one-degree radius from the "sweep" files.
We use the TYPE flag to identify galaxies as all objects whose TYPE = PSF. We use FLUX and MW_TRANSMISSION to calculate dereddened magnitudes, where MW_TRANSMISSION is calculated using the coefficients in Schlafly & Finkbeiner (2011) and the map from Schlegel et al. (1998) . 9 We use SHAPEEXP_R (or SHAPEDEV_R) as our effective photometric radius, R r,eff , for objects that are better fitted with Exponential (or De Vaucouleurs profiles). For objects that are better fitted with a composite profile, we linearly interpolate between these two radii according to the value of FRACDEV to set R r,eff .
We approximate the SDSS grz magnitudes as:
g SDSS = g LS + 0.09, r SDSS = r LS + 0.1, and z SDSS = z LS + 0.02. 9 legacysurvey.org/dr7/catalogs/#galactic-extinction-coefficients These photometric shifts were determined in the several degree region of overlapping coverage between the two surveys. We define "good" galaxy targets in the Legacy Surveys catalogs as objects that satisfy all of the following criteria. Note that each of the following criteria is applied to all three bands (grz) unless otherwise noted:
NOBS ≥ 1 (g and r bands),
RCHISQ < 4 (any one band), FRACIN > 0.7 (any one band), and σ(magnitude) < 0.2.
Even with the above selection criteria, a noticeable number of objects have unrealistic model radii. In most cases, the model is attempting to compensate for nearby sources that were not identified as objects and creates an elongated light profile. Hence we further require our galaxy targets to satisfy:
R ro,eff < min 15 , 10 −0.2(ro−23.5) arcsec, R ro,eff > 10 −0.2(ro−17) arcsec, and
The criteria for selecting "good" galaxy targets from LS catalogs that we described here are rather aggressive -they are needed to avoid a large number of spurious objects in our source catalogs. To ensure our source catalogs are still complete down to our magnitude limit, we only use LS catalogs to supplement SDSS/DES catalogs in this stage of the SAGA Survey (see the next section on the merging process); that is, we do not observe hosts that have only LS coverage. These issues have been mitigated in the recently released LS DR8 and we plan to update our LS photometric catalogs in Stage III.
Sky catalogs are traditionally matched using the nearestneighbor method. While this method is straightforward, it has difficulty dealing with some edge cases, such as when a single object is shredded into different numbers of objects in different catalogs. Low-redshift galaxies tend to be several arcseconds or larger in size, and are often shredded in the above catalogs. Thus, we adapt a different approach to combine available photometric catalogs. We concatenate all objects from all catalogs, and then run a friends-of-friends (FoF) group finder. When the linking length is chosen properly, each FoF group represents one single object in the combined catalog. Within each group, we identify one group member as the primary photometry source for that object. The FoF algorithm and its main advantage over the nearestneighbor method are demonstrated in Figure 2 . Our code for FoF sky catalog matching has been extracted as a standalone, publicly available package. 10 To ensure that we both (1) correctly merge objects that have slightly different positions in different catalogs, and (2) do not accidentally combine distinct objects into one, we use a dynamic linking length for the FoF merging algorithm. We start with a linking length set to 3 ; if there is any group in which there is more than one object from all available surveys, we reduce the linking length by 0.5 and repeat an FoF search on that group, and repeat until the linking length reaches 0.5 or until at least one survey has only one object in that group.
As demonstrated in Figure 2 , each FoF group represents a single object in our final catalog, but it may include multiple photometric objects from multiple surveys. We identify one primary photometry source for that object. Among the choices, we follow the order of DES, LS, and SDSS to select the primary source. If within the group there are multiple good objects from the same survey, we select the object with the brightest magnitude. We use the galaxy properties (magnitude, radius) reported in the primary photometric catalog.
After the FoF matching process, we supplement our source catalog using the NASA-Sloan Atlas (NSA) catalog v1.0.1 (Blanton et al. 2011). 3 The NSA catalog is a reprocessing of the SDSS photometry for galaxies with z < 0.15 using an improved background subtraction technique important for large, extended galaxies. If a merged source catalog object is identified in the NSA v1.0.1 catalog, we use the NSA photometry as the primary photometric information. These objects include many of our host galaxies and their brighter satellites.
In the photometric catalogs, spatially extended galaxies are often shredded into multiple objects. Galactic cirrus or diffraction spikes from bright stars can be mis-identified as galaxies. Although our FoF-merging method described above mitigates many of these problems, spurious objects still contaminate our source catalogs. This issue exists in all three photometric surveys.
We clean these remaining spurious objects both using available spectroscopic data (automated removal) and by eye (manual removal). For automated removal, we first identify all galaxies that have a confirmed spectroscopic redshift and have an error on effective photometric radius less than 50%; we then remove objects that are within two times the effective photometric radius of such galaxies. If the confirmed galaxy is in the NSA catalog, we use the elliptical Petrosian Figure 2 . Illustration of the friends-of-friends (FoF) algorithm with variable linking lengths, used to merge multiple photometric catalogs (Section 3.4). The galaxy shown in the leftmost column is identified as a single object in some photometric catalogs (in this example, DES; top panel in the 2 nd column), and as multiple shredded objects in other photometric catalogs (SDSS and LS; middle and lower panels in the 2 nd column). The numerical indices are ordered by the object magnitudes. If we choose only one particular photometric catalog as the "main catalog" and do a nearest-neighbor match, we would inherit the shredded objects (middle and lower panels in the 3 rd column, or mismatch to a shredded object instead of the main object due to slight difference in position (top panel in the 3 rd column). Note that DES is not always the best-performing catalog; sometimes SDSS or LS performs better, and sometimes all of them produce shredded objects. By using the FoF matching method, all these objects are joined together, and the brightest objects in each photometric catalog are chosen to be the main object for that catalog (rightmost panel).
light profile instead of the simple circular radius. Since these ellipses sometimes overlap with other distinct background galaxies, we only remove objects that do not have redshift information or have a redshift within 200 km s −1 of the said galaxy's heliocentric velocity. As a final check, we visually inspect all objects in our primary targeting region (defined in Section 4.1) and manually flag the handful of spurious objects that remain after automated cleaning. We developed an image list tool 11 based on the cutout service of Legacy Surveys Viewer 12 to facilitate the visual check.
For each object in our cleaned photometric catalogs, we calculate the effective surface brightness µ ro,eff defined as: where r o is total magnitude and R r,eff is the effective halflight radius in r-band from the primary photometric survey, in units of arcseconds. Across the photometric surveys described in Section 3, all of our objects have measured gand r-band photometry. We therefore use Chilingarian et al. (2010) and g − r color to K-correct the r-band absolute magnitude to z = 0, and estimate stellar mass based only on these two bands. We begin with the Bell et al. (2003) relation to determine the mass-to-light ratio based on the g − r color. We adopt a Kroupa IMF and subtract 0.15 dex from the Bell et al. value which assumes a diet Salpeter IMF. We compare these stellar masses to recent estimates from Zibetti et al. (2009) and Taylor et al. (2011) . We find an additional −0.15 dex is consistently needed to match to these two data sets. We calculate stellar mass assuming an absolute solar r-band magnitude of 4.65 (Willmer 2018) . The resulting stellar mass conversion is: log[M * /M ] = 1.254 + 1.098 (g − r) o − 0.4 M r,o . We as-sume a systematic error of 0.2 dex which is larger than the random errors due to propagating errors in photometry and distance.
We compare our merged and cleaned photometric catalogs with various sources to ensure their completeness and validity. For SAGA fields that are in the DES footprint, we compare to Tanoglidis et al. (2020) who re-analyzed lowsurface-brightness galaxies using GALFITM. We confirm that all low-surface-brightness galaxies analyzed in Tanoglidis et al. (2020) and in the SAGA footprint are present in the SAGA photometric catalogs. We find that the surface brightness value derived from DES public catalog are consistently 0.25 mag brighter than the surface brightness reported by Tanoglidis et al. (2020) , but for consistency with other catalogs do not correct our values.
We also compare our object catalogs that are in the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) footprint with the HSC-SSP Public Data Release 2 (Aihara et al. 2019 , overlapping with nine SAGA complete fields, two of which we compared here) and the HSC-SSP low-surfacebrightness galaxy catalog compiled by Greco et al. (2018 Greco et al. ( , 2020 , overlapping with three SAGA complete fields, all of which we compared here). We find that our object catalogs indeed include all real galaxy objects in both HSC catalogs down to our magnitude limit of r o < 20.75. In other words, partly thanks to the inclusion of deep DECam data (DES and LS), our object catalogs do not miss low-surface-brightness objects at the catalog level.
For a host to be considered for SAGA targeting, it must have sufficient imaging data within a projected radius of 300 kpc. Our photometric coverage criteria for targeting are (1) greater than 99% imaging coverage in at least one of the photometric surveys (SDSS, DES, LS), and (2) if the best coverage comes from LS, the area must also have at least 85% imaging coverage from SDSS or DES. We apply the second criterion because we have used a set of more aggressive cuts to remove problematic photometric objects in LS (Section 3.3), and the supplemental SDSS or DES coverage helps ensure that we do not systematically remove good objects. Of the 205 MW analogs identified in the SAGA volume, 115 hosts have sufficient photometric coverage for spectroscopic targeting (Figure 1 ).
In SAGA Stage II, we have improved our target selection strategy (Section 4.1), and obtained 12,690 galaxy redshifts (Section 4.2-4.4), bringing the total number of SAGAobtained redshifts within 1 • around 36 SAGA hosts to 25,372. We also compile 21,075 additional redshifts from the literature (Section 4.5). We describe the association of spectroscopic and photometric objects in Section 4.6, and the resulting redshift distribution in Section 4.7.
Improving our target selection efficiency is one of the major goals of Stage II of the SAGA Survey. Our aim is to complete a host with one or two pointings (200-400 targets) with MMT/Hectospec or AAT/2dF while maintaining high satellite completeness. Achieving this goal requires carefully characterizing the photometric properties of very low-redshift galaxies (the redshift range of SAGA hosts is 0.003 < z < 0.013), a regime in which traditional photometric redshifts are not optimized. There exist very few galaxies at these low redshifts with spectroscopic confirmation or distance measurements that are dimmer than r o = 18, and the SAGA Survey aims to be complete down to r o = 20.75. Hence, it is necessary that we develop our own target selection strategy.
In SAGA Survey Stage I, we employed a set of color cuts (the gri cuts; Equations 1 and 2 of Paper I) to exclude objects that are extremely unlikely to be very low-redshift galaxies. These cuts still leave over 1,000 targets per host on average. In that work we attempted to obtain spectroscopic confirmation for all galaxies satisfying these cuts (within 300 kpc to the host and down to r o = 20.75). As described in Paper I, we achieved more than 90% spectroscopic coverage around six hosts, and more than 80% around the remaining two hosts. This effort resulted in an invaluable training set that allows us to develop new strategies to efficiently identify very low-redshift galaxy candidates; we have applied our strategies to this work and to the low-redshift galaxy program of the Southern Stellar Stream Spectroscopic Survey (Li et al. 2019) . This training set also enable development of machine learning-based methods to select low-redshift galaxies based on photometry (e.g., Wu 2019; N. Dalmasso, in prep.).
In SAGA Stage II, we improve our target selection strategy by adapting a two-fold approach, and devote about half of our spectroscopic resource to each part. We first identify a "primary targeting region," using simple cuts in the surface brightness-magnitude and color-magnitude planes:
where µ ro,eff is the effective surface brightness defined in Equation 2, σ µ is the error on µ ro,eff , and σ gr ≡ σ 2 g + σ 2 r is the error on the (g − r) o color. We only use (g − r) o color as the Legacy Surveys photometry does not include i-band. This primary targeting region is designed to encompass all z < 0.015 galaxies, based on all available redshifts, especially those taken during SAGA Stage I where we obtained Left: N(z) split by target samples. The orange (outermost) histogram shows all galaxies brighter than ro = 20.75 that have redshifts from any source; about 60% of the total to that limit. The blue (middle) histogram shows all galaxies that satisfy the Paper I gri cuts. The dark green (inner) histogram shows all galaxies that are in our primary targeting region, as specified in Equation 3. The inset shows a zoom-in version (0 < z < 0.05) of the same histograms. SAGA satellites are in the redshift range (0.003 < z < 0.013) indicated by the dotted vertical lines. Right: N(z) split by redshift source and apparent magnitude. The yellow (outer) histogram shows all galaxies brighter than ro = 20.75 that have redshifts from the literature and pre-existing surveys. The purple (inner) histogram shows those from the SAGA Survey. The hatched (darker) portion denote galaxies brighter than ro = 19. In both panels and inset, the histograms are normalized so that the height shows the number of galaxies per square degree per redshift bin. a large number of redshifts without any photometric selection. Hence, the probability of a photometric object in this region being a satellite is very high, and we devoted about half of our spectroscopic resource to achieve high spectroscopic coverage for all "good" galaxy targets (as defined in Section 3.1-3.3) within this region. Down to our magnitude limit of r o = 20.75, the number of targets that are within the primary targeting region and between 10-300 kpc to each host in projection is 160 per host on average (∼ 200 targets per square degree). We do not attempt to target galaxies within 10 kpc of the host galaxy because that region is dominated by the host galaxy light. The other half of our spectroscopic resources are devoted to address the potential issue of selection bias and to ensure satellite completeness. This is done by obtaining redshifts for objects outside of the primary targeting region, with various selection methods, including random selection, human selection, and selection based on photometric properties (i.e., selecting objects with similar properties as other low-redshift galaxies). For example, for each pointing that we observe, we allocate 50 or more fibers in "discovery mode," where they are randomly assigned to targets within Paper I gri cuts but outside of the primary targeting region defined in Equation 3. To date, among more than 20,000 redshifts outside of the primary targeting region, we have not found any satellites. All galaxies that have confirmed redshifts in 0.003 < z < 0.015 are within the primary targeting region. This result demonstrates the robustness of our targeting strategy.
The left panel of Figure 3 shows the redshift distribution within (green histogram) and outside of (blue + orange histograms) the primary targeting region among all the Stage I + II redshifts, and compares it with the gri cuts in Paper I (blue + green histograms). The gri cuts effectively remove galaxy targets beyond z ∼ 0.5 but include most galaxies up to z ∼ 0.05. In the inset of Figure 3 , the visible bump around redshift 0.005 < z < 0.01 is due to the presence of satellites in SAGA fields. The definition of our primary targeting region further removes galaxy targets beyond z ∼ 0.3, while including all galaxies within z < 0.015 (and hence all SAGA satellites). In the slightly higher redshift range of 0.015 < z < 0.03, 83% of galaxies are still in the primary targeting region. The small fraction of 0.015 < z < 0.03 galaxies that sit outside of the primary targeting region motivates our two-fold effort to obtain redshifts both within and outside of the primary targeting region, ensuring that we are not missing satellite candidates during the target selection phase.
Stage II of our survey includes redshifts taken with the fiber-based MMT/Hectospec between May 2017 and October 2019. MMT/Hectospec deploys 300 fibers over a one-degree diameter field (Fabricant et al. 2005) . We used Hectospec with the 270 l/mm grating resulting in wavelength coverage of 3650-9200 Å and 1.2 Å per pixel. We use the same low spectral resolution set-ups and data reduction pipelines detailed in Paper I.
Stage II of our survey includes redshifts taken with the fiber-based AAT/2dF observations between June 2018 and May 2020. AAT/2dF deploys 400 fibers over a two-degree diameter field. We used the 580V and 385R gratings in the blue and red arms, respectively, both providing a resolution of R = 1300 (between 1 and 1.6 Å per pixel) over a maximum wavelength range of 3700-8700 Å. We use the same low spectral resolution set-ups and data reduction pipelines detailed in Paper I.
The primary goal of our single-slit work is to observe objects for which we failed to measure a redshift using the facilities above due to low signal-to-noise, bad sky subtraction, or other issues. We obtained data on the 200-inch Hale Telescope at Palomar Observatory with the Double Spectrograph (DBSP) between June 2019 to June 2020. DBSP is a singleslit spectrograph that uses a dichroic to split light into separate, simultaneously observed red and blue channels (Oke & Gunn 1982) . For the blue side, we used the 600/4000 grating at an angle of 27.3 • . For the red side, we used the 316/7500 grating at an angle of 24.6 • . This results in wavelength coverage of 3800-5500 Å and 5500-9500 Å, respectively for the blue and red sides.
We include spectroscopic data in the regions of our hosts from a variety of publicly available surveys. We include spectroscopic data from SDSS DR14 (Abolfathi et al. 2018 ), NSA v1.0.1 (Blanton et al. 2011) , and GAMA DR3 (Baldry et al. 2018) ; these sources were included in Paper I, but now updated versions/data releases are used. In this work, we also include data from the following surveys, listed in decreasing order of number of sources in complete SAGA fields: the WiggleZ Dark Energy Survey (Drinkwater et al. 2018 et al. 1996) . Sources from these surveys are included for quality flag q z ≥ 3, and include a few thousand galaxies in the fields around the SAGA hosts presented here. We additionally include a handful of velocities measured in HI by the ALFALFA survey (Haynes et al. 2018) . We re-target galaxies in cases where a literature value in-dicates a satellite galaxy to ensure homogeneity, but do not count this as a SAGA discovery.
We need to robustly join all the available spectroscopic data and associate these data with the appropriate objects in our photometric object catalogs. While a simple sky coordinate match would serve our purposes in most cases, it does not work as well for extended galaxies, as the separation between the spectroscopic object and the corresponding photometric object can be several arcseconds. Hence, for each spectroscopic object, we search within 20 to find the best-matching photometric object based on the separation and radii of the photometric objects. Each spectroscopic object is preferably matched to the nearest point source that is within 1 , then to the brightest extended source within 20 and within 2 times the effective radius of the extended source. If there is no photometric object that satisfies the above criteria, we match to the nearest photometric object within 20 .
A single photometric object may be associated with multiple spectra from different sources. We adopt the redshift value that has the best reported quality, or, when multiple sources have the same quality, the source with largest aperture. When multiple good redshift values are inconsistent with each other (differ by 150 km s −1 or more), they will not be associated with the same photometric object. The best, nearest redshift value will first be associated with a photometric object, and other redshift values will be matched to the remaining photometric objects that have not yet been matched to a redshift value.
Our code to associate spectroscopic and photometric objects is also publicly available as a component of our survey software. 7 4.7. Distribution of SAGA Redshifts Figure 3 shows the redshift distribution within 300 kpc around all 36 SAGA hosts, split by different targeting regions (left panel; described in Section 4.1) and by redshift sources and magnitude (right panel). Overall, we see that the 25,372 SAGA redshifts (purple histogram) peak around z = 0.2. The majority of SAGA-obtained redshifts are for faint galaxies 19 < r o < 20.75. Conversely, almost all redshifts for dim low-redshift galaxies (19 < r o < 20.75 and z < 0.2) were obtained by SAGA. Specifically, the relative fraction of SAGA to literature redshifts is 46% for bright galaxies (r o < 19) and 82% for faint galaxies (r o > 19). For an in-depth discussion of the SAGA redshift distribution, see Appendix A.
In this section we examine two systematic uncertainties of our survey design. We first review our satellite definition and , satellites newly identified in this work (green squares), and non-satellite galaxies in SAGA fields (orange circles), shown in their heliocentric velocity difference and projected radial distance with respect to their hosts. Satellites are defined to lie within 300 kpc and ±275 km s −1 of their host galaxies in projection; this region is indicated by yellow lines. The dark grey band within 10 kpc indicates the region in which we do not attempt to search for satellites. Right panel: Same as the left panel but plotted with projected distances in units of their mean estimated host halo virial radius and projected velocities in units of their mean estimated host halo escape velocity. These halo quantities are estimated using the abundance matching procedure in Paper I. In both panels, the color maps show the fraction of interloping galaxies, estimated by applying our galaxy-halo connection model to a cosmological simulation (see Section 5.1). The yellow bands show ±1σ escape velocity curves calculated based on the NFW profiles of the potential SAGA host halos in our analysis.
estimate the interloper fraction (Section 5.1). We then define survey completeness (Section 5.2) and carefully develop methods to correct for spectroscopic failures or targeting incompleteness (Section 5.3).
We define a "satellite" as a galaxy that is within a projected distance of 300 kpc of its host and within ±275 km s −1 of the host's redshift, as shown in the left panel of Figure 4 (yellow lines). The projected distance cut roughly corresponds to the virial radius for a typical MW-mass host halo. The velocity cut is slightly larger than the value used in Paper I (±250 km s −1 ), and is motivated by the clear boundary of low-velocity objects seen in Figure 4 . This observationally motivated satellite definition selects galaxies that lie within the three-dimensional virial radius of their dark matter host halo, and conservatively excludes unbound galaxies under the assumption that SAGA host galaxies reside in typical MW-mass dark matter halos. We list these satellites and describe their properties in Section 6.
We apply the same satellite definition when comparing our results with simulations; thus, there is no need to correct for interloping galaxies in our comparisons with theoretical predictions. Nevertheless, we still estimate the fraction of our satellites that may be interloping galaxies in threedimensional space. To calculate this interloper fraction, we populate the c125 − 2048 cosmological gravity-only simula-tion presented in Mao et al. (2015) with galaxies using the subhalo abundance matching (SHAM) technique, which has been thoroughly validated for MW-mass galaxies (e.g., Reddick et al. 2013; Wechsler & Tinker 2018; Cao et al. 2019; Macciò et al. 2020 ). This procedure yields both mock host halos and mock satellite systems (including interlopers) for each SAGA host. We describe this procedure in more detail in Section 8.
We use these predictions to compute the fraction of mock galaxies that lie within the three-dimensional splashback radius of each mock host halo divided by the total number of mock galaxies in the bin, applying an absolute magnitude cut of M r,o < −12.3 to mimic SAGA observations. The result of this procedure is shown by the grayscale map in Figure 4 . We plot the interloper fraction in both "observational units" (i.e., projected distance in units of kpc and velocity in units of km s −1 ; left panel) and in "scaled units" (i.e., projected distance in units of the host halo virial radius, and velocity in units of the host halo escape velocity evaluated at a characteristic distance; right panel).
We find that the galaxies observed by SAGA that pass our satellite definition are very likely to lie within the splashback radius of their corresponding dark matter halo. Here we choose the splashback radius to be 1.5 times the virial radius for the MW-mass host halos; this radius is a physically motivated boundary within which satellite properties are potentially influenced relative to field galaxies (e.g., Diemer & More et al. 2015) . In addition, nearly all galaxies that pass our satellite criteria are predicted to be bound to their corresponding halos, assuming NFW host halo density profiles (yellow bands; Navarro et al. 1996 Navarro et al. , 1997 . Within the SAGA satellite definition, the average interloper fraction is 30%, and the fraction varies with projected distance. A handful of the galaxies we observe lie at low projected distance but high relative velocity with respect to their predicted host halos and therefore deserve further study as potential "fly-by" galaxies (Sinha & Holley-Bockelmann 2012; . Follow-up observations of these galaxies may help distinguish their nature as interlopers or satellites.
In Section 4.1 we designated a "primary targeting region" (Equation 3) in which we aim to achieve high spectroscopic coverage. While we have obtained, and continue to obtain, redshifts for targets outside this region, the chance of finding satellites outside of the primary targeting region is extremely low based on our existing data. Hence, we define a SAGA satellite system as a "complete" system if the spectroscopic coverage is higher than 80% for all galaxies that are both within 300 kpc in projection to its galaxy host and within the primary targeting region, and down to r o = 20.75 (corresponding to M r,o = −12.3 at 40.75 Mpc). Figure 5 demonstrates our completeness definition. The left panel shows the number of targets without redshifts and the spectroscopic coverage within the primary targeting re-gion for all of the 115 SAGA systems passing our targeting criteria (Section 3.8). Among 115 systems, 36 have spectroscopic coverage higher than 80% (and 22 have coverage higher than 90%). We have obtained some SAGA redshifts for an additional 19 systems; these systems are not presented in this work except for being shown as the incomplete systems in Figure 5 .
The list of 36 complete SAGA systems is given in Table 2 . Each of our complete system is assigned a "SAGA name" for internal reference. The table also lists, for each SAGA system, the spectroscopic coverage within the SAGA primary targeting region, the number of satellites (defined above in Section 5.1), and the estimated number of missed satellites (defined below in Section 5.3). Table 2 is also available in machine-readable format on the SAGA website. 1 The right panel of Figure 5 shows the average spectroscopic coverage and the number of targets without redshifts as a function of absolute magnitude. For this plot, we compute absolute magnitude as if all targets were at the distance of their host. All 36 complete systems reach 100% spectroscopic coverage within the primary targeting region for galaxies brighter than M r,o = −15.5. For galaxies fainter than M r,o = −15.5, we maintain on average a 90% spectroscopic coverage down to M r,o = −12.3. For comparison, we show the average spectroscopic coverage from the literature and pre-existing surveys before SAGA; the coverage declines rapidly below M r,o = −17.5, and is not higher than 20% below M r,o = −14.5. . SAGA targets and satellites in the surface brightness µr o ,eff vs. magnitude ro plane (left) and the color (g − r)o vs. magnitude ro plane (right). Confirmed SAGA satellites from the 36 hosts are plotted as green squares or, if they were known from Paper I, as green triangles. Nonsatellite galaxies that have confirmed redshifts z < 0.02 are plotted as orange dots. All other galaxies are shown as small blue dots if they have confirmed redshifts, or as small grey dots if they do not. Brown empty circles represent diffuse low-surface-brightness galaxies (dLSBGs) that we targeted but did not obtain a secure redshift for. The number of brown empty circles is representative of the expected number of unconfirmed satellites. The boundaries of the "primary targeting region" (as defined in Section 4.1, Equation 3) are shown as purple dashed lines.
Figure 5 also shows that, on average, there are about 10 galaxies per host with −15.5 < M r,o < −12.3 for which we do not have redshifts. If all targets within the primary targeting region have exactly the same probability of being identified as a satellite, then we expect there would be ∼ 0.2 satellites for every 10 objects. In the next section, we will utilize a more sophisticated method to estimate and to correct for incompleteness.
While quantifying the spectroscopic coverage and the number of incomplete targets provides a robust way to evaluate our survey completeness, it is important to note that not all targets are equally likely to be satellites. Figure 6 clearly demonstrates this statement by showing the distributions of satellites and targets in the the surface brightness-magnitude plane (left panel) and the gr color-magnitude plane (right panel). The boundaries of our primary targeting region are shown as purple dashed lines.
As stated in Section 4.1, all confirmed satellites are in the primary targeting region, where we find 127 satellites among 5,419 redshifts. Outside of the primary targeting region, we find no satellites among 21,308 redshifts. If we assume that the rate of finding a satellite within (or outside of) the primary targeting region is constant, then its value would be ∼ 2% (or < 0.005% for objects outside).
Of course, the chance of finding a satellite should vary smoothly rather than behave like a step function in photometric space. We hence build a model to evaluate "satellite rates," denoted by R sat , and defined as the fraction of satel-lites among the targets with redshifts at a specific point in the 3-dimensional photometric space (r o , µ ro,eff , (g − r) o ). We describe the details of our model construction in Appendix B, and calculate the satellite rate R sat for each object for which we have not obtained a redshift. These satellite rate values can then be used to correct for our incompleteness; the sum of R sat would be the expected number of unconfirmed satellites.
There is an additional concern about the situation in which we failed to measure redshifts for faint, diffuse, low-surfacebrightness objects using multi-fiber spectroscopy facilities; this occurs when the flux incident on the fiber is too low (i.e., for objects with low fiber magnitudes). While we do attempt to follow up these targets with single-slit spectroscopy, we do not always recover redshifts due to our limited single-slit resources. In Figure 6 , we show these "diffuse" Low-Surface-Brightness Galaxies (dLSBGs) for which we failed to obtain redshifts with AAT or MMT as brown empty circles; they cluster in the faint, low-surface-brightness, and slightly redder regime. Note that within the same photometric regime, we did successfully obtain many redshifts (blue dots), but these are mostly for LSBGs that have slightly higher central surface brightness and are at slightly higher redshift than dLSBGs. While the distinction between the diffuse and lessdiffuse LSBGs are often discernible by the human eye (see examples in Figure 7) , they are not easily separable using photometric catalog properties alone (Tanoglidis et al. 2020) .
To ensure that we have a conservative estimate of our completeness, we visually inspected targets in each SAGA field to identify all dLSBGs. Among the 36 SAGA fields, we . Images of low-surface-brightness galaxies (LSBGs) with similar photometric properties. The four columns, from left to right, show (i) diffuse LSBGs (dLSBGs) that are confirmed satellites, (ii) low-redshift dLSBGs that are not satellites, (iii) dLSBGs that we failed to obtain redshifts for due to low fiber flux, and (iv) less-diffuse LSBGs that have similar magnitudes, colors, and surface brightnesses in the photometric catalog and slightly higher redshifts. Images taken from Legacy Surveys DR8. 12 The width and height of each image correspond to 40 .
found 70 of these objects that we have not obtained redshifts for (down to M r,o < −12.3); 13 of them are slightly outside the primary targeting region. The counts of these dLSBGs without redshifts are shown in Figure 5 as error bars on the left panel and as a dark histogram on the right panel. Based on the redshifts that we have obtained for some dLSBGs and on the radial distribution of these dLSBGs (i.e., accounting for the fact that satellites concentrate closer to host galaxies while background low-redshift galaxies should be uniform), we estimate that about 25-30% of these dLSBGs may be satellites.
We incorporate this dLSBG estimate in our satellite rate model, and then calculate R sat for each target for which we have not obtained a redshift. Among them, the average R sat for targets that are within the primary targeting region but not dLSBGs is ∼ 0.005; while the average R sat for dLSBGs is ∼ 0.3. The sum of these R sat values gives us the expected total number of missed satellites. We found that on average we may be missing about 0.7 faint satellites (−15.4 < M r,o < −12.3) per SAGA host (24 in total among 36 hosts). All SAGA hosts are close to 100% complete above M r,o = −15.4.
In the analysis below, we interpret R sat as the probability of each target without a redshift being a satellite, and we incorporate this correction in binned statistics such as the luminosity function and quenched fraction. For scatter plots that show distributions of individual satellites, we include a representative set of unconfirmed potential satellites. This set contains 24 potential satellites in −15.4 < M r,o < −12.3 and 11 below M r,o = −12.3, selected by a random realization based of R sat . Potential satellites are always marked distinctly from confirmed satellites.
We now turn to the main findings from SAGA Stage II. In this section, we investigate properties of individual satellite galaxies. In Section 7, we study these galaxies as satellite systems, and discuss the implications of our findings for the galaxy-halo connection in Section 8.
In the 36 complete systems, we have identified 127 satellites using the definition in Section 5.1. There are 123 within the absolute magnitude limit of the Survey, M r,o < −12.3. The 127 satellites are distributed around 33 hosts; 69 of these satellites did not have pre-existing redshifts. One host, NGC 4348, has only one confirmed satellite that is fainter than the magnitude limit, and the other 2 complete systems (NGC 5896, NGC 5633) have zero confirmed satellites. We find two satellites around ESO 288-025 that are also within 300 kpc in projection to NGC 7166 and have consistent redshifts (recall that these two hosts are a Local Group-like pair). The two satellites sit in the outskirts of NGC 7166, and hence we assign them to ESO 288-025. Figure 8 shows the images of all 127 satellites. The images are taken from the Legacy Surveys Viewer 12 with the DR8 layer (this layer includes DES imaging where available). Quenched satellites (defined in Section 6.2) and satellites fainter than our absolute magnitude limit M r,o = −12.3 are indicated by a red or cyan triangle in the upper corners, respectively. The coordinates and basic properties of these satellites are listed in Table 3 and are available in machinereadable format on the SAGA website. 1 We also plot the positions and velocities of all SAGA satellites in Figure C. 1, Appendix C.
For satellites that were reported in Paper I, some properties may differ slightly in this work, due to our adoption Figure 9 . Effective surface brightness vs. physical radius (left) and color (right) for all satellites in our 36 complete SAGA hosts. Confirmed satellites are color-coded based on their stellar mass, and a representative set of potential satellites are shown as empty black circles to demonstrate the photometric regions where we may be incomplete. MW satellites are shown as stars for comparison. In the left panel, galaxies within the dotted lines are considered"ultra-diffuse" galaxies (UDGs).
of new photometric catalogs and the rebuilding of host list.
In addition, we have discovered two more satellites around two of the Paper I 'complete' hosts: LS-444338-475 in the system of NGC 6181 and LS-311554-3218 in the system of PGC 068743. We excluded one satellite reported in Paper I (1237666408439677694 in the system of NGC 7716), due to its poor redshift quality; however, this object is below our survey magnitude limit. We also obtained an accurate redshift for NGC 5962b, which was reported to be a satellite of NGC 5962 in Zaritsky et al. (1997) , and found its redshift at 0.338, and hence it is not a satellite of NGC 5962.
In Figure 6 , we showed how SAGA satellites distribute in the surface brightness-magnitude and color-magnitude planes. In Figure 9 , we plot the distributions of satellite effective surface brightness against physical effective radius and color, and color-code each satellite point based on its stellar mass. In Figure 10 we plot the distributions of satellite color, effective surface brightness, absolute magnitude and projected distance to the host in four different views (the color code of Figure 10 is explained in Section 6.2). In both Figure 9 and 10, we include a representative set of potential satellites that do not yet have measured redshifts (as empty circles) based on our satellite rate model prediction.
In Figure 9 , we include the MW satellites using quantities from references listed in Paper I. In the left panel, the SAGA satellites show the same general trend (color gradient) as the MW satellites: intrinsically fainter satellites are on average lower surface brightness and have smaller physical radii. In the lower left panel of Figure 10 , the surface brightness relation appears to flatten near the faint end, even when those unconfirmed satellites are considered. While this flattening behavior may be due to poorly measured photometry, it is nevertheless consistent with the satellites in the Local Group (McConnachie 2012, Figure 7) .
Five of the SAGA satellites are formally in the ultra-diffuse regime (defined by van Dokkum et al. 2015 as R eff > 1.5 kpc and µ g,0 > 24 mag arcsec −2 , which translates to µ ro,eff > 24.7 mag arcsec −2 assuming a n=1 Sérsic profile (Graham & Driver 2005) and average (g − r) o = 0.4). This region is in the upper right corner of Figure 9 left panel (shown as dotted lines), and includes the Sagittarius Dwarf Spheroidal Galaxy of the MW. None of the SAGA satellites approach the size or surface brightness of M31's ultra-diffuse satellite And XIX although this galaxy is below our absolute magnitude limit (Collins et al. 2020) . We also find a small number of ultradiffuse galaxies in our larger low-redshift galaxy sample.
The right panel of Figure 9 shows no clear trend between color and surface brightness, especially if one considers the set of unconfirmed potential satellites. We find that the color distribution of MW satellites is similar to that of SAGA satellites. In the upper left panel of Figure 10 we plot the color-magnitude relation of SAGA satellites. Due to the luminosity-metallicity relationship (Kirby et al. 2013) , we might expect a trend of bluer colors with decreasing satellite luminosity. We do not see this trend, especially when considering the unconfirmed satellites, likely due to the competing trend of increased quenching for the fainter satellites. We will discuss this in detail in the next subsection.
We measure the equivalent width (EW) of the Hα line in our spectra for all SAGA satellites; the Hα line probes star formation within ∼ 5 Myr (Flores Velázquez et al. 2020). We define quenched satellites as having EW(Hα) < 2 Å, follow- ing several authors including Geha et al. 2012 . Since the majority of our spectra are fiber-based, and the fiber coverage is usually a small fraction of the overall galaxy, our EWs cannot be directly translated into star formation rates without correcting for aperture effects. Based on our definition, a galaxy labeled as star-forming is secure, while a quenched galaxy could be star forming in the rare cases where our fiber missed a star-forming region. We color-code the satellite points in Figure 10 based on the measured value of EW(Hα). The representative set of potential satellites do not have spectra, and are shown as empty circles. The color bar in Figure 10 transitions from blue to red at EW(Hα) = 2 Å to indicate quenched satellites in red. Down to our magnitude limit of M r,o = −12.3, there are 18 quenched satellites out of 123 confirmed satellites (15%). This is slightly higher than the 1/27 quenched satellites (4%) presented in Paper I; these numbers are statistically consistent given the small number of satellites presented in Paper I. Since the AAT/MMT fibers cover only a fraction of each galaxy, it is possible we missed a star-forming region in a few cases. All quenched satellites are indicated in Figure 8 by a red triangle in the upper left corner.
As mentioned in Section 6.1, we do not see a clear colormagnitude trend, likely due to the competing trend of increased quenching for the fainter satellites (upper left panel of Figure 10 ). Nevertheless, we do see that quenched satellites tend to be redder when compared with satellites of the same absolute magnitude. There is a modest trend between color and projected distance from the host, such that inner satellites have redder colors (upper right panel of Figure 10 ). Quenched satellites also tend to be closer to their host galaxies; they have a median projected distance of 110 kpc, more centrally concentrated than the median projected distance of 180 kpc for star-forming satellites. Figure 11 shows the satellite quenching fraction as a function of stellar mass for SAGA satellites, MW and M31 satellites (Wetzel et al. 2015) , and for isolated SDSS field galaxies (Geha et al. 2012) . For SAGA satellites, we show the correction due to incompleteness and interlopers (as faint green error bars) in addition to the shot noise (as dark green, capped error bars). The incompleteness correction conservatively assumes all potential satellites are quenched. The interloper correction downweights each confirmed satellite based on their interloper probabilities. Figure 11 . Satellite quenching fraction as a function of stellar mass for SAGA satellites (green points with error bars), MW and M31 satellites (orange crosses with band; Wetzel et al. 2015) , and isolated SDSS field galaxies (purple band; Geha et al. 2012) . See Section 6.2 for quenching definition. For SAGA satellites, the darker error bars with caps denote shot noise, while the fainter error bars denote the combined correction for both incompleteness and interlopers (dominated by incompleteness correction). The quenching fraction increases as stellar mass decreases for both SAGA and Local Group satellites, but the increase is steeper for Local Group satellites.
In Figure 11 , the quenching fraction increases as stellar mass decreases for both SAGA and Local Group satellites, but the increase is steeper for Local Group satellites. All but two of our quenched satellites are fainter than M r,o = −16, with stellar masses less than 10 8.2 M . Interestingly, the two bright quenched satellites are both very close to their respective hosts. Around M * = 10 8.5 M , almost all SAGA satellites are still star-forming, while more than half of Local Group satellites (dominated by M31 satellites) are quenched. At M * = 10 7 M , the quenching fraction goes to about 0.6 for SAGA satellites if we consider incompleteness and interloper corrections, while the majority of Local Group satellites are quenched at this stellar mass.
We note that in Figure 11 , quenching is defined as EW(Hα) < 2 Å for SAGA satellites and field galaxies, and as M gas /M * < 0.1 for Local Group satellites. This difference in quenching definition would not reconcile the discrepancy, as a few of our Hα-based quenched galaxies may be star-forming due to fiber coverage, but gas-depleted galaxies would not be star-forming. Related work by Bennet et al. (2019) measured the quenching fraction of the satellite galaxies of four Local Volume galaxies using the presence of bright blue main sequence stars as an indicator of star formation, and found that, down to the SAGA magnitude limit, the quenching fraction is about 0.6-0.7 for M81 and Cen A, but zero for M94 and M101.
These apparently different quenching fractions around different host galaxies may indicate that host environment or other host properties impact the star formation of Local Group satellites. The rapid variation ("burstiness") of star formation in dwarf galaxies may also cause large scatter in quenching fraction from host to host (e.g., Weisz et al. 2012 ). Comparing our results with numerical simulations will shed light on the relation between the observed quenching fraction and quenching timescales (e.g., Wheeler et al. 2014; Fillingham et al. 2015; Garrison-Kimmel et al. 2019b; Akins et al. 2020 ). In particular, Akins et al. (2020) suggest that the satellite quenching timescale can vary significantly depending on satellite gas mass and the level of ram-pressure stripping. On the observational side, we plan to obtain redshifts for unconfirmed satellites to reduce our error bars on the quenching fraction, and to measure star formation rates for each of our satellites including aperture corrections (e.g., Brinchmann et al. 2004 ).
We now investigate satellite population-level statistics of SAGA systems, focusing on satellite luminosity functions (Section 7.1), radial distributions (Section 7.2), and satellite groups and planes (Section 7.3).
We examine the luminosity functions of SAGA satellite systems. The left panel of Figure 12 shows the cumulative number of satellites around each SAGA host galaxy as a function of absolute magnitude, including our assessment of host-by-host incompleteness as described in Section 5.3. This panel demonstrates the diversity of SAGA satellite systems despite our relatively restrictive range of host properties. It further shows that the MW luminosity function is typical among our statistical sample.
To assess this quantitatively, we compute the Poisson likelihood of observing the MW sample given the set of SAGA systems, marginalizing over the unknown Poisson rate in each absolute magnitude bin, following Nadler et al. (2019b) . We find that the likelihood of observing the MW satellite system is nearly identical to that of a model which simply averages the observed count of SAGA satellites across all systems in each absolute magnitude bin (the log-likelihood difference between these scenarios is ∼ 1). Thus, the MW satellite luminosity function is consistent with being drawn from the same underlying distribution as SAGA data.
In contrast, the un-normalized satellite luminosity function of M31 lies in the high-satellite-number tail of the SAGA luminosity function distribution; two SAGA systems have the same number of satellites as M31. Note that M31 is slightly brighter than our upper limit on host luminosity. We also find that 8 out of 36 SAGA systems have no confirmed satellites brighter than M r,o = −12.3 within 150 kpc of their hosts. Thus, the M94 satellite system (see Smercina et al. 2018 ) is also not an outlier in the SAGA context.
In Figure 13 , we explore the correlation between the host galaxy luminosity and the amplitude of the satellite luminosity function by showing the number of satellite galaxies down to our absolute magnitude limit of M r,o = −12.3 as a function of their host's K-band absolute magnitude, including estimates of host-by-host incompleteness based on the model described in Section 5.3. Total satellite counts for the MW and M31 are also included. We observe a significant correlation between host M K and total satellite count, with a Spearman rank correlation coefficient of ρ s = −0.39 (p-value = 0.018); this correlation is a factor of ∼ 10 more significant than the trend reported in Paper I, largely due to the increased sample size.
The left panel of Figure 12 reveals that several SAGA hosts have bright/massive satellites that are of similar luminosities or stellar masses of the Large and Small Magellanic Clouds (LMC/SMC). About 30% of the 36 SAGA systems contain at least one satellite similar to the LMC (more massive than M * = 10 9 M or M r,o < −18.3), and about 25% of SAGA systems have at least two satellites more massive than M * = 3 × 10 8 M (or M r,o < −17; similar to SMC). Considering that the search radius of SAGA is 300 kpc in projection, our result is consistent with the findings of Liu et al. (2011) and Tollerud et al. (2011) . SAGA systems that host bright/massive satellites also tend to have higher total numbers of satellites; the Spearman rank correlation coefficient between the M r,o value of the brightest satellite and the total satellite count (M r,o < −12.3) is ρ s = −0.66 (p-value = 2 × 10 −5 ).
In the right panel of Figure 12 , we inspect the average cumulative satellite luminosity functions in scaled magnitude, ∆m ≡ m (sat) − m (host) . The scaled magnitude is calculated in r-band, and we use NSA photometry for 30 SAGA hosts that are in the NSA catalog. For the remaining 6 hosts we use M K + 2.5 as an approximation. For comparison, we show several cumulative satellite luminosity functions from the literature, including COSMOS (Nierenberg et al. 2012 ; hosts selected from 0.1 < z < 0.4 and 10.5 < log[M * /M ] < 11), SDSS (Sales et al. 2013 ; hosts selected from z < 0.055 and 10 < log(M * /M ) < 11), and MW-like hosts in the Local Volume (Carlsten et al. 2020b) . Note that the SAGA hosts are selected from 0.003 < z < 0.013 and 10 < log(M * /M ) < 11.
We caution that this comparison should be considered exploratory rather than conclusive, especially given the different host selections and satellite definition. Nevertheless, the shapes of these satellite luminosity functions in scaled magnitude are surprisingly consistent among these studies, hinting at a universal satellite luminosity function for MW-mass host galaxies across different redshifts. The SAGA result is the first to probe the satellite luminosity functions outside of the Local Volume to this faint-end limit, and we may be seeing a slight hint of flattening at the faint end of the luminosity function. Note that the expected ∼ 0.7 unconfirmed satellite per SAGA host is already included in the figure as part of the green error bars.
The amplitude of satellite luminosity functions, on the other hand, is very sensitive to the mass of the host galaxy/halo. For example, in the right panel of Figure 12 For the latter, we calculated the projected distance using a random set of sightlines, and the 1σ range is shown as shaded areas around the median lines. Also shown are the average radial distribution of the eight SAGA hosts that have the most concentrated distribution (thin green line), and a spherical distribution of r −1 projected onto 2D (dotted blue line).
between Local Volume systems and the SAGA systems is again most likely due to different host stellar masses. The median stellar mass of SAGA hosts is about 10 10.4 M , while that of the MW-like hosts presented in Carlsten et al. (2020b) is about 10 10.7 M . Recently, Roberts et al. (2020) studied satellite luminosity functions in the COSMOS field using a statistical approach similar to Nierenberg et al. (2012) , and found a much higher amplitude (by about a factor of five) than the luminosity functions of SAGA, Nierenberg et al. (2012) , and Sales et al. (2013) . Hence, a closer examination on how the amplitude of satellite luminosity function depends on host selection is warranted.
The average cumulative projected radial distribution of confirmed SAGA satellites is shown in Figure 14 . Note that we have not attempted to identify satellites within 10 kpc of the hosts as the region is dominated by host light; hence the cumulative function starts at 10 kpc. Since all SAGA satellites are within 300 kpc by definition, the normalized cumulative radial distribution reaches 1 at 300 kpc. The observed distribution is roughly linear in projected distance, similar to the findings of Tollerud et al. (2011) , and is consistent with a hypothetical case where satellites are distributed according to 1/r 3D within a sphere (shown as dotted blue line). This 1/r 3D profile is less concentrated than a typical NFW profile.
For comparison, we show the radial distributions of galaxies around the MW and the M31. We use the galaxies tabulated in McConnachie (2012) , assuming M V = M r,o + 0.4 (Jester et al. 2005 ) and excluding galaxies that are fainter than SAGA magnitude limit (M r,o = −12.3). We then produce projected radial distributions by "mock observing" these Local Group galaxies with 5,000 sightlines that are sampled uniformly at random. We reject any sightline that results in the MW and the M31 being within 300 kpc in projection with each other, and label any Local Group galaxies that are within 300 kpc in projection to the MW/M31 and within ±500 kpc in the light-of-sight direction as "satellites." Both these operations are done to mimic the SAGA observations. The median distributions for MW and M31 are shown as dashed orange and dash-dotted grey lines respectively in Figure 14 , and the corresponding bands show the 1σ spread due to random projection.
We see that the radial distribution of MW satellites is much more concentrated than the average distribution of SAGA satellites, mostly due to the presence of the LMC and SMC. The radial distribution of M31 satellites is also slightly more concentrated than that of SAGA satellites, but they are marginally consistent. We note that this difference is most likely due to host-to-host scatter. In Figure 14 we also show the radial distribution averaged over only the eight (25%) of the SAGA hosts that have the most concentrated radial distributions; the resulting radial distribution is almost as concentrated as the MW satellite radial distribution.
To compare our results with the Local Volume satellite systems presented in Carlsten et al. (2020c) , we calculate d proj,half , the median projected distance (i.e., radius encompassing half of the satellites), for each SAGA system. Among the 36 SAGA systems, the median and the 16 th /84 th percentiles of d proj,half are 169 +69 −52 kpc. If we limit our sample to only satellites within 150 kpc to their hosts to be consistent with Carlsten et al. 2020c , the median and the 16 th /84 th percentiles of d (<150 kpc) proj,half become 83 +48 −36 kpc. 13 This result can then be directly compared with Figure 4 of Carlsten et al. (2020c) ; we find the distribution of d proj,half among SAGA systems is about 1σ larger than that of the Local Volume satellite systems (including) MW and M31, but is in full agreement with simulation predictions (see also Section 8.2). Our finding is also consistent with the results of Samuel et al. (2020) , which showed that the radial distribution of SAGA Paper I satellites is in agreement with the Feedback in Realistic Environments (FIRE) baryonic simulations, and that the host-to-host variation dominates the scatter. 13 Note that in each SAGA system, there are only ∼ 2 satellites within 150 kpc on average; hence the measurement of d Finally, we find no strong dependence between the radial distribution and satellite magnitude. This can been seen in the lower right panel of Figure 10 .
Over the past decade, the "plane of satellite" questioni.e., whether satellite systems around MW-like galaxies preferentially form a co-rotating planar structure-has raised much discussion (see review by Pawlowski 2018) . Detailed studies have been done on three satellite systems: the MW, M31, and Centaurus A (e.g., Pawlowski et al. 2012 Pawlowski et al. , 2013 Tully et al. 2015) . Additionally, statistical analysis using SDSS data has been done by several groups. In this work, we present an analysis similar to that in Ibata et al. (2014) and Phillips et al. (2015) . A more detailed analyses is expected in Stage III of the SAGA Survey. For each SAGA system, we identify all the satellite pairs (a pair is any two satellites of the same host) and calculate the "opening angle" with respect to their host galaxy. An opening angle of 0 • indicates the two satellites in a pair are on diametrically opposite sides, while an opening angle of 180 • indicates they are on the same side of the host and have exactly the same position angle. For each satellite pair, we check if their spectroscopic velocities have a co-rotating signature, that is, whether the two satellites are moving toward the opposite (same) direction with respect to the host when the opening angle is less (greater) than 90 • . We plot the fraction of pairs that have the co-rotating signature as a function of opening angle in Figure 15 . We only consider satellites that have a velocity difference with respect to their respective hosts greater 25 km s −1 .
Among the SAGA satellite systems presented here, we have identified 36 satellite pairs for satellites within 150 kpc to their hosts and with ∆V > 25 km s −1 . The number of satellite pairs increases to 189 for satellites between 150-300 kpc. The green boxes and blue band in Figure 15 show the co-rotating fraction for the inner satellite pairs and the full satellite pairs, respectively, in 10 bins of opening angle. We calculate the fraction using the mean Bayes estimator, (c + 1)/(n + 2), where n is total number of pairs in that bin, and c is the number of co-rotating pairs. The error bars are estimated using the Wald method, p(1 − p)/n, where p is the estimated fraction. The signal is much noisier in this case of inner satellite pairs, but we show it for comparison with the SDSS result, which is computed inside 150 kpc and presented in Figure 3b of Phillips et al. (2015, shown as a grey line).
With this analysis, we do not observe a clear excess of high co-rotating fraction, either for the full satellite sample or using only satellites within 150 kpc. For pairs that have small opening angles (on nearly diametrically opposed sides), the co-rotating fraction is close to, and slight less than 0.5. However, we note that the SAGA signal is noisier than the SDSS signal, due to the limited number of pairs. A possible corotating excess seems to appear in the last bin near 180 • opening angle (on the same side of the host). Those "corotating" pairs in fact are composed of satellites that are close to each other (less than 100 kpc in projection) and have similar velocities (see Figure C .1 for an illustration). Hence, it is possible that this excess near 180 • opening angle is due to satellite groups rather than satellite planes. Further analysis is needed to fully distinguish the origin of this excess.
To place our results in the context of galaxy-halo connection models, we combine the abundance-matching procedure that associates SAGA host galaxies with dark matter halos with the empirical satellite galaxy model from Nadler et al. (2019b Nadler et al. ( , 2020a . This model populates subhalos with satellite galaxies by extrapolating an abundance-matching relation between luminosity and peak maximum circular velocity (calibrated to the GAMA luminosity function from Loveday et al. 2015 for M r,o < −13) and a satellite size relation (calibrated to Kravtsov 2013 for galaxies with half-light radii r 1/2 1 kpc) into the regime of faint satellites. In addition, it includes a model for satellite disruption due to a central disk potential, which is calibrated to the FIRE baryonic simulations. We update the disruption prescription from Nadler As described in Nadler et al. (2019b Nadler et al. ( , 2020a , this model accurately describes the luminosity, size, and radial distribution of MW satellite galaxies, including the population of satellites accreted with the LMC, when combined with recent observational MW satellite selection functions . To incorporate theoretical uncertainties in our prediction for SAGA satellite populations, we sample from the posterior distribution over galaxy-halo connection and baryonic disruption model parameters presented in Nadler et al. (2020a) , which is derived by fitting the model to most of the known MW satellite population.
To predict the satellite population for each SAGA host, we first follow the procedure in Paper I to map each SAGA host galaxy to a set of potential dark matter host halos in the c125 − 2048 cosmological simulation (a higher resolution version of the box used in Mao et al. 2015) using abundance matching. We assume that the scatter in central luminosity at fixed halo properties is 0.15 dex, and we select all halos mapped to absolute magnitudes within 0.15 mag of the SAGA host galaxy in question, which roughly corresponds to the quadrature sum of the estimated photometric and distance errors for our hosts. We also impose the host M K Figure 17 . SAGA satellite luminosity functions and incompleteness corrections (colored lines and bands) compared to predictions from a cosmological dark matter-only simulation populated with galaxies using the empirical satellite model in Nadler et al. (2019b Nadler et al. ( , 2020a , which has been fit to the MW satellite population. Dark blue lines indicate the mean prediction for each satellite population, and blue (light blue) contours indicate 68% (95%) confidence intervals, which include the effects of host galaxy-halo abundance matching scatter, uncertainty in our galaxy-halo connection model (see Section 8.1), and projection effects.
) and the environment cut (Equation 1d) described in Section 2.1. 14 This procedure yields ∼ 2,000 halos per SAGA host galaxy on average, from which we select a random subset of 300 matched halos for each model realization to capture the probabilistic relation between host galaxy and host halo properties. In particular, these mock host halos have a cosmologically representative range of masses and formation histories that is expected for each host galaxy lu- 14 We do not implement the stellar foreground cuts in this mock selection because our simulation does not include stars or a Galactic disk. We do not implement the maximum halo mass cut because the abundance matching procedure yields very few potential SAGA host halos with present-day virial mass > 10 13 M . minosity, which allows us to estimate the uncertainty in satellite populations resulting from the scatter in these quantities at fixed central luminosity.
For each potential host halo in a given model realization, we then select all other halos above a resolution threshold of V peak = 40 km s −1 that satisfy the ±300 kpc and ±275 km s −1 projected distance and velocity criteria in various projections of our simulation to mimic the SAGA satellite definition; thus, our model self-consistently includes interloping galaxies. We populate (sub)halos with (satellite) galaxies using the model described above, and we measure the "observed" satellite population for each host halo matched to a given SAGA host galaxy for several draws of satellite model parameters from the Nadler et al. (2020a) posterior. We perform 0 1 N( Figure 18 . Comparison of predicted and observed radial satellite distributions, normalized to the number of satellites within a projected distance of 300 kpc. The predicted mean and confidence interval is identical to Figure 17 , and the observed radial distributions are computed using all satellites above our Mr,o < −12.3 absolute magnitude limit. Hence, no observed data are shown for systems that do not have any satellites with Mr,o < −12.3. We do not correct the SAGA radial distributions because our incompleteness model predicts that missing satellites should have the same distribution as confirmed ones.
three sets of mock observations of galaxies in the z = 0 snapshot of the simulation by projecting the simulation along perpendicular axes, and imposing the M r,o < −12.3 SAGA completeness threshold and a surface brightness completeness limit of µ ro,eff < 26 mag arcsec −2 . 15 Mock-observed satellites are weighted by their survival probability according to our baryonic subhalo disruption model, and mock-observed interlopers are assumed to have survival probabilities of unity. 15 The model predicts a very small number of satellites of µ ro,eff > 26 mag arcsec −2 ; however, there are no satellites nor satellite candidates at this low surface brightness in the SAGA object catalog or in the matched HSC catalog.
With these procedures, the model predicts that, at 68% confidence, SAGA hosts inhabit dark matter halos with virial masses between 7 × 10 11 < M vir /M < 2 × 10 12 and that SAGA satellites are hosted by halos with peak virial masses between 2 × 10 10 < M peak /M < 10 11 , and typical presentday virial masses a factor of ∼ 1-6 lower. Figure 16 shows the comparison between the observed and predicted distributions of the total number of satellites (M r,o < −12.3) among the 36 complete SAGA systems. The error bars on the observed number are generated using the procedure described in Section 5.3. To generate the predicted distribution, we repeatedly draw sets of mock host halos corresponding to the 36 SAGA systems; thus, the "predicted" contour includes both statistical uncertainties (due to the limited host sample size) and systematic uncertainties (due to scatter in the host galaxy-halo connection and our satellite model). Since this plot shows predicted and incompletenesscorrected realizations of the SAGA sample, adjacent data points are anti-correlated. We find good agreement between the total predicted N sat distribution and SAGA observations.
To compare our predictions to the data in more detail, Figure 17 shows predicted luminosity functions and Figure 18 shows predicted radial distributions (normalized to the number of satellites with M r,o < −12.3 within a projected distance of 300 kpc) compared to the observed radial distribution for each complete SAGA system. Again, contours indicate 68% and 95% confidence intervals due to scatter in the host galaxy-halo connection, draws from our satellite galaxy model parameters, and projections of our simulation. Our predictions are largely consistent with the observed luminosity functions, and they are in excellent agreement with the observed normalized radial distributions. However, for SAGA hosts with the largest numbers of observed satellites, our model under-predicts both the total number of satellites and particularly the number of bright (M r,o < −15) satellites, for which SAGA observations are highly complete (e.g., see the right panel of Figure 5 ).
To quantify this bright-end tension, we estimate that 0.13 (0.09) predicted satellites with M r,o < −15 must be added to our fiducial prediction per host to bring it into agreement with the data at 68% (95%) confidence assuming Poisson errors on the observed counts, which corresponds to a ∼ 1-2σ discrepancy. This is consistent with Figure 17 , which shows that ∼ 5 SAGA hosts have ∼ 1 additional bright satellite relative to our predicted 95% confidence interval. There is also a hint that the model over-predicts the number of dim satellites, although it is formally consistent with the data given our current incompleteness estimates.
Overall, it is encouraging that the model predictions are in broad agreement with SAGA data, reinforcing our finding that the MW satellite population is not highly atypical.
Before we discuss the implications of the bright-end tension noted above, we emphasize several caveats associated with our current predictions. The satellite model we employed has specifically been fit to the MW satellite population; thus, our predictions for subhalo and satellite galaxy disruption do not account for the varying masses, morphological properties, and host halo density profiles of SAGA host galaxies relative to the MW. In addition, we have not included "orphan" satellites in the model, nor have we explored exactly how our predictions depend on galaxy-halo connection parameters within the region of parameter space allowed by the MW satellite population.
Nonetheless, the potential bright-end discrepancy is reminiscent of similar tensions noted for Local Volume field (Neuzil et al. 2020 ) and satellite (Carlsten et al. 2020b ) galaxies when compared with models that differ from ours in detail, hinting at a more systematic issue that may exist in various theoretical predictions, including those from hydrodynamical simulations. The flexibility of our model allows us to quantify the possible sources of the tension and to study potential solutions. We therefore briefly describe possible solutions in the context of our model, leaving a more thorough investigation to future work:
1. Stellar mass-halo mass relation: Forcing the halos that host the richest SAGA satellite systems to be a factor of ∼ 3 more massive resolves the bright-end discrepancy. However, this would significantly alter the well-validated stellar mass-halo mass relation in the MW-mass regime (Wechsler & Tinker 2018) . 16 Alternatively, this halo mass shift can be achieved if SAGA host magnitudes are systematically biased dim by ∼ 0.5 mag, which we also regard as unlikely; 2. Disruption model: Removing subhalo disruption from our model reduces the bright-end tension to less than ∼ 1σ. It is unlikely that subhalos which host bright satellites undergo no disruption; however, disruption prescriptions calibrated to hydrodynamical zoom-in simulations are dominated by low-mass subhalos and therefore might over-predict disruption efficiencies for more massive subhalos;
3. Global luminosity function: Our abundance matching prediction is calibrated to the GAMA luminosity function down to M r,o = −13; this luminosity function contains both statistical errors (captured by uncertainties in the Schechter function fit to GAMA data) and potential systematic errors (e.g., due to survey incompleteness). While there is no evidence that the GAMA Survey is incomplete down to r = −19.8, we find that varying the GAMA luminosity function amplitude from Loveday et al. (2015) within its quoted 2σ error can fully resolve the bright-end tension.
In future work, we plan to address these questions in detail and to explicitly test for consistency between the MW and SAGA satellite populations by re-fitting the satellite model to SAGA data. We also plan to compare SAGA results with additional models and hydrodynamical simulations.
In this work we present the Stage II results from the ongoing SAGA Survey, including 36 spectroscopically complete satellite systems around Milky Way analogs. This release marks the completion of just over one third of our planned total (Table 1 ). The full redshift data set around these hosts includes 25,372 SAGA-obtained redshifts, peaking around z = 0.2, and 28,425 redshifts from the literature and preexisting surveys. For galaxies in 17.5 < r o < 20.75, 77% of those redshifts are first obtained by SAGA.
Several improvements to the observing strategy and analysis have been made compared to SAGA Stage I (Paper I).
We discussed these survey improvements in Section 2-5 and summarize selected highlights here:
1. We have updated our list of Milky Way analog hosts, and the 'Grand List' from which these hosts are selected, incorporating updates in HyperLEDA and EDD. We now identify 205 hosts that meet our criteria, 115 of which have full photometric coverage at present (Figure 1 ).
2. Our photometric sources now include SDSS DR14, DES DR1, and Legacy Surveys DR6/DR7. We have developed a friends-of-friends merging procedure to combine overlapping photometric coverage (Figure 2 ). Checking against the HSC-SSP, we have verified that we do not miss low-surface-brightness objects in our photometric catalogs.
3. Using redshifts collected from Paper I, we have improved our target selection strategy to efficiently select very low-redshift galaxy candidates (Figure 3, 6) . 4. We have developed a robust method to evaluate our survey completeness ( Figure 5 ), based on both spectroscopic coverage and an accurate model that characterizes of the rate of satellite discovery as a function of photometric properties (Figure B.1, B.2) . This satellite rate model is used to correct for survey incompleteness.
Our primary science results are presented in Section 6-8, focused on the properties of satellites and satellite systems. We summarize them as follows:
1. We identify 127 satellites around 36 SAGA hosts (Figure 8;  2. SAGA satellites follow a surface brightness vs. magnitude relation that is consistent with that of the MW (Figure 9 ). We find a small number of "ultra-diffuse" galaxies in our satellite sample.
3. Consistent with Paper I, the satellite quenching fraction among SAGA systems is lower than that in the Local Group ( Figure 11 ). The quenched fraction increases with decreasing stellar mass: 2/50 bright satellites are quenched (M r,o < −16, similar to LMC/SMC) while 16/55 faint satellites (−16 < M r,o < −12.3) are quenched. We also see a slight trend of increasing quenching fraction as projected radius decreases (Figure 10 ).
4. The satellite luminosity function of the MW is consistent with being drawn from the same distribution as the SAGA systems ( Figure 12 ). The average SAGA luminosity function is in agreement with past studies but probing into fainter regime. The amplitude of the satellite luminosity function is sensitive to host mass ( Figure 12 ).
5. The total number of satellites (M r,o < −12.3) per host ranges from 0 to 9. This number modestly correlates with host galaxy K-band luminosity (ρ s = −0.39, p-value = 0.018; Figure 13 ), and more strongly correlates with the brightest satellite magnitude M r,o (ρ s = −0.66, p-value = 2 × 10 −5 ; Figure 12 ).
6. The satellite radial distribution of SAGA systems is much less concentrated than that of the MW and to a lesser extent M31, but the difference can be largely accounted for by host-to-host scatter ( Figure 14) . Among SAGA satellite systems, we find no evidence for corotating planes of satellites ( Figure 15 ). 7. Our measured total satellite number, satellite luminosity functions, and radial distributions are all largely consistent with predictions based on a ΛCDM+SHAM fit to the Milky Way satellite population, and in particular inform the galaxy-halo connection in the halo mass regime of 10 9 -10 10 M ( Figure 16 , 17, 18).
As more satellite systems are observed and characterized, both by the SAGA Survey and by others, we will be equipped with a powerful tool to test the ΛCDM model and galaxy formation theory. In this work we have touched on several topics, and already found intriguing results such as the wide range of radial distributions and the low quenching fraction of faint satellites, but much remains to be explored. As we make concrete progress towards the completion of SAGA Survey Stage III, we plan to follow up on these questions. In addition, we are publishing our current satellite data in machine-readable format, 1 hoping to foster independent studies of the SAGA satellite systems.
Perhaps the most remarkable aspect of the SAGA results is how this exhibit of 36 satellite systems around Milky Way analogs solidifies the idea that our very own satellite system of the Milky Way is just one "realization" from a diverse distribution. While we can now say that the Milky Way seems to sit comfortably among these SAGA systems in terms of number of satellites, it is interesting to imagine if the Milky Way were like 'Odyssey' (NGC 6181; nine satellites) or 'Don-Quixote' (NGC 5633; no satellites) -the development of galaxy formation theory might have proceeded quite differently than it has in the past few decades. The fascination of seeing how our knowledge evolves as we build a statistical sample is at the heart of the SAGA Survey.
respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory (NOAO); the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOAO. The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.
NOAO is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.
BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program "The Emergence of Cosmological Structures" Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 11433005).
The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration.
The The public archival data from the DES is based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.
This project used public data from the GAMA Survey. GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is www.gama-survey.org.
This work made use of public data collected at the Subaru Telescope and retrieved from the Hyper Suprime-Cam (HSC) data archive system, which is operated by Subaru Telescope and Astronomy Data Center (ADC) at National Astronomical Observatory of Japan. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica In-stitute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.
This work made use of images from the Digitized Sky Surveys (DSS). The Digitized Sky Surveys were produced at the Space Telescope Science Institute (STScI) under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions. The National Geographic Society -Palomar Observatory Sky Atlas (POSS-I) was made by the California Institute of Technology with grants from the National Geographic Society. The Second Palomar Observatory Sky Survey (POSS-II) was made by the California Institute of Technology with funds from the National Science Foundation, the National Geographic Society, the Sloan Foundation, the Samuel Oschin Foundation, and the Eastman Kodak Corporation. The Oschin Schmidt Telescope is operated by the California Institute of Technology and Palomar Observatory. The UK Schmidt Telescope was operated by the Royal Observatory Edinburgh, with funding from the UK Science and Engineering Research Council (later the UK Particle Physics and Astronomy Research Council), until 1988 June, and thereafter by the Anglo-Australian Observatory. The blue plates of the southern Sky Atlas and its Equatorial Extension (together known as the SERC-J), as well as the Equatorial Red (ER), and the Second Epoch [red] Survey (SES) were all taken with the UK Schmidt. All DSS data are subject to the copyright given in the copyright summary. Copyright information specific to individual plates is provided in the downloaded FITS headers. Supplemental funding for sky-survey work at the STScI is provided by the European Southern Observatory.
We acknowledge the usage of the HyperLeda database (leda.univ-lyon1.fr) and the Extragalactic Distance Database (EDD; edd.ifa.hawaii.edu). Support for the development of content for the EDD is provided by the National Science Here we examine in detail how SAGA redshifts and targets distribute in the photometric space. Figure A. 1 provides a close look at the numbers of galaxies that are within the footprints of both SAGA and SDSS, and whether they have only SAGA redshifts, have pre-existing redshifts (e.g., from SDSS), or have no redshift information. We plot these numbers with respect to three photometric quantities: apparent magnitude (r o ), effective surface brightness, (µ ro,eff ), and color ((g − r) o ).
In the left panel of Figure A .1, galaxies within the SAGA footprint and brighter than r o = 16 basically all have preexisting redshifts. Between the magnitude range of 16 < r o < 17.5, about half of the galaxies in our primary targeting region (bars with dark colors in the figure) do not have pre-existing redshifts, and SAGA have obtained redshifts for these galaxies. This result implies that even in the magnitude range of 16 < r o < 17.5, very low-redshift galaxies (z < 0.03) do not have complete spectroscopic coverage. The incomplete spectroscopic coverage in the very low-redshift universe is consistent with the findings of Lazo et al. (2018) and Fremling et al. (2020) . Below the magnitude of r o = 18, most redshifts were, unsurprisingly, obtained by SAGA; similar information has also been shown in the right panel of Figure 3 .
In the middle and right panels of Figure A .1, we see more clearly that SAGA has obtained thousands of redshifts for galaxies outside of the primary targeting region, covering about 10-20% of total objects. These redshifts provide us with the quantitative evidence that our primary targeting region does not exclude very low-redshift galaxies or potential satellites. In these panels we also see how the satellites (orange histogram) distribute in these photometric quantities. We note that in the photometric region where satellites are more populated, the majority of redshifts were obtained by SAGA. . Color indicates whether the galaxy has no redshift (grey), has a pre-existing redshift (blue), has a SAGA redshift (green), or is a satellite (orange). For each of these categories/colors, brightness indicates whether the galaxy is in our primary targeting region (dark) or not (light); the primary targeting region is defined in Equation 3. Note that the y-axis scale is linear between 0 and 10, and logarithmic above 10. Hence, above 10, the visual height difference between histograms should be interpreted as multiples; the vertical bars on the left panel demonstrate the height differences corresponding to a few common multiples.
In Section 5.3, we calculated the "satellite rate," the fraction of satellites among the target with redshifts, as a function of photometric properties. We detail our construction of this satellite rate here. The construction involves a modified logistic regression method; however, despite the name, the core ideas behind the model are straightforward.
Our goal is to build a smooth function in photometric space where the value of the function represents the satellite rate. We first choose to work in the photometric space of (r o , µ ro,eff , (g − r) o ). We then assume that (1) there is a latent variable which takes the form of a linear combination of the photometric properties:
and (2) the satellite rate (probability) is a monotonic function of the latent variable . In other words, we aim to identify a specific direction in the photometric space, determined by (β 1 , β 2 , β 3 ), for which the satellite rate varies monotonically only along that direction. The monotonic function we choose to model the satellite rate is a logistic function,
where the satellite rate R sat ( → −∞) = 0 and R sat ( → ∞) = R max . With this setup, the remaining task is to fit for the parameters (β 0 , β 1 , β 2 , β 3 , R max ). We use the maximumlikelihood method, and the log-likelihood can be calculated exactly as:
This model construction is almost identical to logistic regression, except that we also fit for R max , rather than fixing it to unity. This is because it is impossible to find a region in the photometric space where one would only find satellites, and hence the maximum satellite rate should be capped at some R max < 1. Operationally, we modify the Logit model in the STATSMODELS PYTHON package 17 (Seabold & Perktold 2010) to carry out the maximum-likelihood fit. We fit this model to all galaxy redshifts we have collected among 36 complete SAGA systems. The best-fit parameter values are (β 0 = 0.303, β 1 = −1.96, β 2 = 1.507, β 3 = −5.498, R max = 0.487). To estimate the errors, we bootstrap the satellite population and re-fit the parameters 1000 times. We choose to not bootstrap the full population because the satellites are too rare in the full population. We also fix the R max value in the error estimation process.
We note that the model we built here aims to be descriptive rather than predictive, as there are multiple complications if one desires to build a predictive model for rare events (e.g., King & Zeng 2001) . The main goal here is to construct a model to reproduce the "satellite rate" in different regions of Targets ( sat < 10 3 ) Targets ( sat 10 3 ) Figure B .2. Distribution of galaxies in the surface brightness-magnitude (left) and color-magnitude (right) planes, color-coded by its "satellite rate" (Rsat) produced by our model. The Rsat value can be read off from the right-most color bar. Galaxies whose Rsat < 10 −3 are shown as small grey points to avoid crowding. The boundaries of our primary targeting region are shown as purple dashed lines.
the photometric space, and to smooth out the noise due to limited data points. Figure B .1 demonstrates that our model can indeed reproduce the input data set. The satellite rates our model produces (cyan bands) agree very well with the data in different slices of the photometric space. In the left panel, the satellite rate decreases with magnitude, which is a simple result of the relation between distance, volume, and magnitude. In the middle and left panels, we observe that the model can accurately capture the satellite rate as a function of surface brightness and color, even if the function is not monotonic.
With this model, we can then calculate satellite rate (R sat ) at any given point in the 3-dimensional photometric space.
We demonstrate our calculation in Figure B .2, where each point represents a target, and its green-blue color corresponds to our estimate of its satellite rate R sat . These values are then used to correct for incompleteness, as discussed in Section 5.3.
For completeness, we plot the positions and velocities (with respect to their respective hosts) of all SAGA satellites in Figure C .1. It is noticeable that some close satellite pairs also have similar velocities with each other, and contribute to the excess near 180 • in Figure 15 . Table 2 , with background image from the Digitized Sky Survey (DSS-1) and 1 • on a side. The xand y-axes show Right Ascension and Declination, respectively. The large grey circle indicates a radius of 300 kpc (virial radius) to each host. Satellites are plotted as small circles, color coded by the difference in heliocentric velocity with respect to their hosts. 