© ESO 2022

1. Introduction

According to the standard cosmological cold-dark-matter-dominated Universe with a cosmological constant (ΛCDM model), the structures formed hierarchically, resulting in massive galaxies that had grown through accretion of or merging with smaller ones. Today, massive galaxies should be surrounded by numerous smaller dark matter halos that harbor satellites. A typical Milky Way galaxy should be able to gravitationally bind any low-mass object interior to 250−300 kpc. A direct expectation would be a positive correlation between the number of satellites and the total mass of the host galaxy, because more massive host halos harbor more subhalos (Ishiyama et al. 2013).

Although the ΛCDM model is very successful in explaining the large-scale structure distribution (Eisenstein et al. 2005), on smaller, more local scales, there are many challenges revealed in the Local Group, such as the missing satellites problem (Klypin et al. 1999), the too-big-to-fail problem (Boylan-Kolchin et al. 2011), and the plane of satellite problem (Pawlowski & Kroupa 2013; Müller et al. 2018). Galaxies in the Local Group are close enough to enable deep observations of faint structures in detail. Another recent challenge to the ΛCDM was posed when empirical correlations were found between the number of satellites and: (1) the mass of the bulge (Javanmardi et al. 2019), and (2) the bulge-to-total ratio (Javanmardi & Kroupa 2020) for several nearby galaxies. These two correlations are highly unlikely in the ΛCDM model, as the number of satellites should only correlate with the mass of the dark matter halo and therefore with the rotation velocity, which is closely related to dynamical mass. Javanmardi et al. (2019) used data from Milennium-II simulations (Boylan-Kolchin et al. 2009) to test the significance of these empirical correlations with the predictions of the ΛCDM model, and found no significant correlation. For the former correlation, the linear correlation coefficient from simulations was measured to be 0.32, and for the latter, the correlation coefficient measured was 0.13. Empirical correlations are established for only a few nearest galaxies, and therefore a deep survey is needed to complete satellite population around Milky Way-like galaxies. There is also the question of whether the number of satellites is correlated to some other fundamental galactic property. This information would suggest additional mechanisms guiding galaxy evolution.

Great effort has been made to complete a satellite survey of a large sample of Milky Way-like galaxies outside the Local Group. The Satellites Around Galactic Analogs (SAGA) survey¹ (Geha et al. 2017; Mao et al. 2021) was conceived in order to discover satellite galaxies around 100 Milky Way analogs brighter than M_r, 0 = −12.3 (the SAGA limit, hereafter). At the present time, at Stage II of the project, 36 SAGA galaxies have been surveyed within their virial radii (∼300 kpc) where all plausible candidate satellites have been spectroscopically confirmed. Host galaxies (SAGA galaxies, hereafter) are morphologically diverse, ranging from early-type (ellipticals and lenticulars) to late-type galaxies (Sb and Sc).

The sample is large enough to enable the characterization of a luminosity function of satellites, but also allows us to test correlations between the number of satellites and the physical properties of their host galaxies. This sample can also be used to test predictions of the ΛCDM model on small scales outside of the Local Group.

In Sect. 2, we analyze the dependence between the baryonic mass of galaxies and their near-infrared (NIR) flux, discuss the dynamical mass estimation, and present measurements of angular momentum based on the contribution of individual galactic components (bulge and disk). In Sect. 3, we give a description of observational data used for the study: NIR images that were used to infer stellar mass and data retrieved from various databases on gas mass and galaxy kinematics. In Sect. 4, we present the correlations between the number of satellites and fundamental galactic properties of their host galaxies, and finally in Sect. 5, we discuss the implications of the correlations found on galaxy formation and evolution and present our conclusions.

2. Fundamental galactic properties

There are several fundamental galactic properties that are nontrivial to estimate: baryonic mass, dynamical mass, and total specific angular momentum. Baryonic mass, as a sum of the stellar mass and gas mass, has been shown to be more important than the stellar mass itself (McGaugh et al. 2000). Moreover, it has been shown to correlate, through the bulge-to-total baryonic mass ratio (B/T), with the total number of satellites (Javanmardi & Kroupa 2020). The most difficult to estimate is the stellar mass, because it is measured from luminosity that cannot be converted to mass in a straightforward manner.

Dynamical mass reflects the overall gravitational potential of galaxies, including the contributions from baryonic mass and dark matter. Assuming the rotation velocity dominates velocity dispersion in the outer parts of the galaxy, then the gravitational force balances the centrifugal force, and we can estimate an enclosed mass within the radius at which rotation velocity is measured. Single-dish HI spectra are available for a large number of nearby galaxies, and therefore the width of the HI spectral line is used as a proxy for rotation velocity. The radius of the HI disk is often several times larger than the optical radius and therefore gives a very good estimate of circular motions at large distances from the galaxy center.

The angular momentum of galaxies is conserved during the collapse of material within a hierarchical framework of galaxy formation, but is systematically lower in early-type galaxies than in spirals (Romanowsky & Fall 2012). This can be explained by mergers, whereby the angular momentum from central regions of galaxies is redistributed to their outer parts by dynamical friction. Therefore, in early-type galaxies, angular momentum could be locked up in the outer envelope populated with satellites.

2.1. Baryonic mass

Gas mass is mostly neutral hydrogen (ℳ_HI) with small contributions from heavier elements. Baryonic mass is the sum of stellar mass (ℳ_*) and the gas mass (ℳ_HI):

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{bar}}={\mathcal{M}}_{\ast }+1.33{\mathcal{M}}_{\mathrm{HI}},\end{array}$$(1)

where the contribution of heavier elements is accounted for by the factor 1.33 (Courteau et al. 2014). The mass of neutral hydrogen is measured as an integral under the HI line spectrum; that is, the total flux (F_HI) multiplied by the distance (D) squared:

$$\begin{array}{c}\hfill \frac{{\mathcal{M}}_{\mathrm{HI}}}{{\mathcal{M}}_{\odot }}=2.33\times {10}^5\left[\frac{F_{\mathrm{HI}}}{{\mathrm{Jy}\mathrm{km}\mathrm{s}}^{-1}}\right]{\left[\frac{D}{\mathrm{Mpc}}\right]}^2.\end{array}$$(2)

Estimating the stellar mass of a galaxy is a complex and widely debated problem. There are two main approaches, which rely on extraction of the visible (baryonic) mass from the dynamical mass estimate, or on the stellar population models that correlate the observable properties, such as luminosity and color, to the galaxy stellar mass. Colors can be measured from the optical part of the spectrum or, more favorably, the NIR part, which is dominated by the light of an old stellar population that makes the bulk of the stellar mass, and is least susceptible to extinction.

In distant galaxies, where individual stars cannot be resolved, the integral light coming from both young and old stars is dominated by that from the young stars in the visible part of the spectrum, while old stars dominate the stellar mass. This problem can be solved by using the NIR light, which favors the old stellar population. Therefore, from the luminosity in the NIR (Spitzer 3.6 μm), stellar mass can be estimated as ${\mathcal{M}}_{\ast }={\mathrm{{\rm Y}}}_{\ast }^{[3.6]}{L}_{[3.6]}$. The NIR mass-to-light ratio ${\mathrm{{\rm Y}}}_{\ast }^{[3.6]}$ needs to be constant with a typical value of 0.47 ℳ_⊙/L_⊙ to achieve self-consistency with optical observations (McGaugh & Schombert 2014). The mass-to-light ratio can vary throughout the galaxy and thus affect the mass estimation of individual galaxy components, such as the bulge and disk. Eskew et al. (2012) calibrated 3.6 and 4.5 μm fluxes to measure the stellar mass using the spatially resolved star formation history of the Large Magellanic Cloud (Harris & Zaritsky 2009):

$$\begin{array}{c}\hfill {\mathcal{M}}_{\ast }={10}^{5.65}{F}_{[3.6]}^{2.85}{F}_{[4.5]}^{-1.85}{\left[\frac{D}{0.05}\right]}^2,\end{array}$$(3)

where F_3.6 and F_4.5 are fluxes in 3.6 and 4.5 μm measured in Jy, ℳ_* is in solar masses, and D is a distance in Mpc. The scatter they obtained corresponds to 30% uncertainty in the mass estimate.

Their correlation is important because it provides a way to estimate stellar mass for a large sample of galaxies with NIR photometry taken by Spitzer (Sheth et al. 2010) or WISE (Wright et al. 2010) telescopes. In the NIR, there is only marginal sensitivity to young stellar populations, and measurements can be only slightly effected by dust. We therefore use the above empirical formula to measure stellar masses of SAGA galaxies from the Spitzer NIR images.

2.2. Dynamical mass

The dynamical mass of a galaxy with a rotationally supported HI disk is calculated as (Yu et al. 2020):

$$\begin{array}{c}\hfill \frac{{\mathcal{M}}_{\mathrm{dyn}}}{{\mathcal{M}}_{\odot }}=2.31\times {10}^5{\left[\frac{V_{\mathrm{rot}}}{{\mathrm{km}\mathrm{s}}^{-1}}\right]}^2\left[\frac{R_{\mathrm{HI}}}{\mathrm{kpc}}\right],\end{array}$$(4)

where R_HI is the radius of HI disk. As we have no information on HI disk size from single-dish HI spectra, we use an empirical correlation (Wang et al. 2016) that holds for gas masses over a large range (10^5.5 < M_HI/M_⊙ < 10^10.5):

$$\begin{array}{c}\hfill log\left(\frac{R_{\mathrm{HI}}}{\mathrm{kpc}}\right)=0.51log\left(\frac{{\mathcal{M}}_{\mathrm{HI}}}{{\mathcal{M}}_{\odot }}\right)-3.59.\end{array}$$(5)

A better estimate of dynamical mass comes from detailed modeling of rotation curves. However, from single-dish HI spectra, only the HI line width can be measured and converted to rotation velocity applying the inclination correction.

2.3. Specific angular momentum

An expression for the specific stellar angular momentum valid for spiral galaxies was revised by Romanowsky & Fall (2012) and in the case of an infinitely thin disk with an exponential radial profile is given by:

$$\begin{array}{c}\hfill j_{\mathrm{d}}=2\left[\frac{V_{\mathrm{rot}}}{{\mathrm{km}\mathrm{s}}^{-1}}\right]\left[\frac{R_{\mathrm{h}}}{\mathrm{kpc}}\right],\end{array}$$(6)

where R_h is disk scale length, and V_rot is rotation velocity. SAGA galaxies have different morphologies, but exponential disks dominate their outer radial profiles. For the inner parts of galaxies, where the bulge is present, the expression from Romanowsky & Fall (2012) is valid for an arbitrary value of Sersic index n:

$$\begin{array}{c}\hfill j_{\mathrm{b}}=k_n\left[\frac{v_{\mathrm{s}}}{{\mathrm{km}\mathrm{s}}^{-1}}\right]\left[\frac{R_{\mathrm{eff}}}{\mathrm{kpc}}\right],\end{array}$$(7)

where k_n ≃ 1.15 + 0.029 n + 0.062 n², and R_eff is the effective radius of the bulge. The rotation velocity v_s of the bulge can be calculated using the following prescription from Romanowsky & Fall (2012):

$$\begin{array}{c}\hfill v_{\mathrm{s}}={\left(\frac{v}{\sigma }\right)}^{\ast }{\sigma }_0\sqrt{\frac{ϵ}{1-ϵ}},\end{array}$$(8)

where (v/σ)^* stands for the contribution of rotational over pressure support in the galaxy dynamics. Romanowsky & Fall (2012) derived a median value of (v/σ)^* = 0.7 ± 0.4 which can be applied as a single value for large sample of galaxies. Another two parameters are the central velocity dispersion in km s⁻¹ (σ₀) and ellipticity (ϵ).

Knowing the ratio of the bulge mass to the total stellar mass, that is, the stellar B/T ratio, the total specific angular momentum (Romanowsky & Fall 2012) can be expressed as:

$$\begin{array}{c}\hfill j=B/T·j_{\mathrm{b}}+(1-B/T)·j_{\mathrm{d}},\end{array}$$(9)

where j_b is the bulge angular momentum estimated from Eq. (7) and j_d is the disk angular momentum given by Eq. (6). The stellar B/T ratio in the previous equation reflects the relative contribution from a particular structural component to the total angular momentum. In the case where the stellar B/T ratio is small (B/T < 0.1), the disk component dominates, and the bulge contribution can be neglected (Romanowsky & Fall 2012).

3. Data and sample selection

The Spitzer Heritage Archive² contains 32 out of 36 SAGA galaxies in both 3.6 and 4.5 μm-bands. Archived images are fully reduced but not background subtracted. We estimated and subtracted an average of the median sky values inside a dozen boxes of 10 by 10 arcsec in size far away from the galaxy. We then ran the IRAF ellipse task to extract radial surface brightness profiles for all the galaxies. The majority of them were part of the Spitzer Survey of Stellar Structure in Galaxies (S4G; Sheth et al. 2010; Muñoz-Mateos et al. 2013; Querejeta et al. 2015). The S4G sample includes all but three of the SAGA galaxies analyzed here (NGC 6909, NGC 7029, and NGC 5869).

In order to make a direct comparison of the results from the literature with our results, we list here the choices of selections and cuts used. Our sample is based on the SAGA survey and all SAGA galaxies were selected as isolated galaxies, meaning that there are “no galaxies brighter than M_K + 1³ within 1° of the host, and such that the host is not within two virial radii of a massive (5 × 10¹² M_⊙) galaxy in the 2MASS group catalog (Lavaux & Hudson 2011)” (Geha et al. 2017). As Javanmardi et al. (2019) and Javanmardi & Kroupa (2020) compared nearby galaxies including the most massive Local Group galaxies to the Milennium-II simulation results, their selection criteria favor halos with virial masses of between 10¹⁰ and 2 × 10¹² M_⊙ that are isolated and harbor central disk galaxies. The sample of nearby galaxies to which the simulations are compared in this work comprises: the Milky Way, M 101, M 31, M 33, M 94, M 81, and Cen A. These span the following range of stellar mass 1 × 10¹⁰ − 11.2 × 10¹⁰ M_⊙, and a full range in bulge-to-total ratio (0.01 < B/T < 1). The stellar mass of SAGA galaxies goes from 1.8 × 10¹⁰ to 11.8 × 10¹⁰ M_⊙, and bulge-to-total ratios cover values from 0.01 to 0.8. The only caveat to keep in mind is the brightness limit of the satellite populations. Nearby galaxies from the studies of Javanmardi et al. (2019) and Javanmardi & Kroupa (2020) have the satellites brightness limit of M_V ≤ −8.2 mag, while the satellites of the SAGA galaxies are all brighter than M_r, 0 = −12.3 mag, which is the SAGA limit. In the remainder of the paper, we remove all the satellites fainter than the SAGA limit from the sample of nearby galaxies (Javanmardi et al. 2019; Javanmardi & Kroupa 2020) for direct comparison, and use an empirical correlation found between the number of satellites and the bulge-to-total ratio to recover the faint end of the satellite population for SAGA galaxies to enable direct comparison with the total number of satellites.

First, we inspected the inclinations of all SAGA galaxies and removed the three found to be edge-on (NGC 4348, PGC 13646, and UGC 903). The remaining 29 galaxies were used to infer the relation between bulge-to-total ratio and the number of satellites. Hereafter in our analysis, as it is strongly affected by the accuracy of the distance and/or rotation velocity measurements, we applied two empirical relations (stellar mass–size and the baryonic Tully Fisher relation) to remove seven outliers (NGC 1309, NGC 2962, NGC 4454, NGC 5347, NGC 5602, NGC 5869 and PGC 68743). Another two galaxies were discarded from the sample because we were not able to find the data on either their HI mass or rotation velocity, which was needed to estimate dynamical mass. For the remaining 20 galaxies dynamical mass was estimated and correlated with the number of satellites. Finally, for specific angular momentum measurement velocity dispersion was needed to estimate the contribution of the bulge to the total angular momentum, and for two galaxies we lack these data (NGC 5792 and NGC 7328). These two were therefore removed from the correlation between the number of satellites and the specific total angular momentum. For this last correlation, the total number of SAGA galaxies used is 18. These 18 galaxies have passed all empirical tests and have all the data needed. More details on these selection cuts are given in the following section, and the number of galaxies used for each of the correlations is given in Table 1.

Background-subtracted images are modeled with the GALFIT code (Peng et al. 2010) using multiple components including the disk and bulge with the point spread function (PSF) often needed to describe the nucleus, and additional components like bars and rings where needed. GALFIT has become a standard tool for modeling galactic structure. In Fig. B.1, we plot all the individual modeled subcomponents, along with the composite models. In addition to the radial profiles in 3.6 and 4.5 μm-bands, Fig. 1 presents the sum of all individual subcomponents modeled with GALFIT for each galaxy, and is referred to here as a composite model. Individual components are the nucleus described with the point spread function (PSF), and the bulge and disk. Comparison with radial profile of S4G:P3 is indicated for the 3.6 μm band only (4.5 μm omitted) for clarity. Our radial profiles in the 3.6 μm band show excellent agreement with the radial profiles of S4G galaxies in common with the S4G survey extracted from their Pipeline 3 (Muñoz-Mateos et al. 2015), as can be seen for all SAGA galaxies in Fig. B.1.

Fig. 1. Radial profile of the NGC 1015 galaxy in 3.6 and 4.5 μm-bands modeled with the IRAF ellipse task, denoted as blue and red points, respectively. Overplotted are orange crosses from Muñoz-Mateos et al. (2015) for the 3.6 μm-band. For clarity, only composite models are given for 3.6 μm (black) and 4.5 μm bands (gray), that correspond to the sum of all individual subcomponents. Individual subcomponents are presented only for the 3.6 μm band (point spread function, bulge and disk).

The Spitzer images are in units of M Jy sr⁻¹ given in the AB system. To convert total (integral) magnitudes measured by the GALFIT code into fluxes in Jy that are required for stellar mass estimation (Eq. (3)), we used the following conversion (Sheth et al. 2010):

$$\begin{array}{c}\hfill F={10}^{-0.4(m_{\mathrm{AB}}-8.9)},\end{array}$$(10)

where F is the flux in Jy and m_AB is the magnitude in the AB system. With this conversion, Eq. (3) now reads:

$$\begin{array}{c}\hfill \frac{{\mathcal{M}}_{\ast }}{{\mathcal{M}}_{\odot }}=64.78\times {10}^{10}{\left[{10}^{-0.4(m_{[3.6]}-m_{[4.5]})}\right]}^{1.85}{10}^{-0.4m_{[3.6]}}{\left[\frac{D}{\mathrm{Mpc}}\right]}^2.\end{array}$$(11)

Here, distance D to the source should be provided in Mpc, and total (integral) magnitudes m_[3.6] and m_[4.5] are in the AB system. The calculated stellar mass ℳ_* is in solar masses.

As one of the correlations that has been found for nearby galaxies by Javanmardi & Kroupa (2020), namely that between baryonic B/T ratio and the number of satellites, involves gas mass in addition to stellar mass (i.e., baryonic), we searched the Extragalactic Distance Database (EDD)⁴ (Tully et al. 2009) to complement our stellar mass measurements with gas mass estimates. Most of the galaxies are cataloged in the EDD catalog “All Digital HI” with HI flux measured. The rest of the galaxies are listed in the EDD catalog “Pre-Digital HI”⁵, where the measurements of log F_HI flux are given. Three galaxies are also part of the ATLAS3D project (Cappellari et al. 2011): NGC 2962, NGC 5869, and NGC 6278. For these three galaxies, we retrieved data on H₂ mass. Combining stellar mass with gas mass (atomic and molecular) provides us with the total baryonic mass.

From the EDD catalog “All Digital HI” we retrieved the averaged width of the HI line (Wmx_av). The width of the HI line can be corrected to give maximum rotation velocity, which can be used to infer the specific stellar angular momentum and dynamical mass of galaxies. To convert Wmx_av into rotation velocity V_rot = Wmx_av/2sin(i), we need to correct for inclination, which has been estimated from the ratio of the minor to the major galaxy disk axis (b/a): ${sin}^2(i)=1-[{(b/a)}^2-q_0^2]/[1-q_0^2]$, where q₀ = 0.11 (Guthrie 1992). For three galaxies that were part of the ATLAS3D project (Cappellari et al. 2011), namely NGC 2962, NGC 5869, and NGC 6278, rotation velocities were taken from Krajnović et al. (2011).

4. Results

We modeled all SAGA galaxies in the 3.6 μm band with multiple components using the GALFIT code. The results are presented in Appendices A and B (table and surface brightness profiles, respectfully). The baryonic mass is measured as a sum of stellar and gas mass: ℳ = ℳ_* + 1.33ℳ_HI. After checking galaxy inclinations, which were estimated from the ratio of the disk minor to major axis from GALFIT disk models, we excluded three galaxies as edge-on, because this inclination is unfavorable for bulge mass estimation (NGC 4348, PGC 13646, and UGC 903).

For each galaxy component (bulge and disk), total mass depends on the distance squared (Eqs. (2) and (11)), while B/T ratio as the ratio of masses is not affected by distance uncertainties. To enable comparison with Javanmardi & Kroupa (2020), we excluded all the satellites fainter than M_r, 0 < −12.3 for all the host galaxies from Javanmardi & Kroupa (2020) (and references therein). In the paper by McConnachie (2012), we count five satellites of the Milky Way and eight satellites of the Andromeda galaxy brighter than the SAGA limit and up to 250 kpc. For the M 81 galaxy, we refer to Chiboucas et al. (2013), where galaxies from the M 81 group are listed with their physical separation and absolute magnitudes. Within 250 kpc from M 81 galaxy, we count 20 satellites brighter than the SAGA limit. Crnojević et al. (2019) list all satellites of Cen A galaxy within 250 kpc, and we count 10 satellites brighter than the SAGA limit. M 94 galaxy has been observed by Smercina et al. (2018) but only up to 150 kpc distance, and only two satellites had been detected down to M_V = −9.1 mag. As they are both fainter than our limit, we consider that M 94 has no satellites. M 33 galaxy has one possible but unconfirmed satellite fainter than M_V = −6.5 mag, as reported by Martin et al. (2009), and therefore we consider it has no satellites. M 101 galaxy has eight confirmed satellites in total (Carlsten et al. 2019; Bennet et al. 2019), but they are all fainter than the SAGA limit.

In Fig. 2, SAGA galaxies are given in gray along with their error bars (Table 2), while galaxies from Javanmardi & Kroupa (2020) are given in red. As discussed above, for galaxies from Javanmardi & Kroupa (2020) (red circles in the Fig. 2) the number of satellites is limited to those brighter than the SAGA limit to enable direct comparison. If we assume that the relation holds, then all the galaxies that lie below the green area have numerous satellites that were missed by the SAGA survey. We cannot confirm the validity of the relation with these data, because we only include satellites brighter than a certain limit. However, knowing that, when the total number of satellites is taken into account, all red points from Javanmardi & Kroupa (2020) enter the green shaded area, there is a possibility that the relation truly holds. Another sign that the relation may hold is that there are no points above the green shaded area, the presence of which would question the correlation. Simply, there can be only more satellites than detected, not less. All the points can only go up, never down, if the correlation holds, and this is satisfied with the existing data. We label six SAGA hosts that enter the green shaded area with their NGC numbers, and we consider they have already a total number of satellites detected, and thus are directly comparable to the sample of nearby galaxies with complete satellite census. For these six SAGA galaxies, and those from Javanmardi & Kroupa (2020) with total number of satellites, we measure a Pearson’s linear correlation coefficient of 0.73.

Fig. 2. Relation between the B/T ratio and the number of satellites brighter than Mr, 0 < −12.3 (SAGA limit). Red points are taken from Javanmardi & Kroupa (2020) and references therein. Their values on the y-axis are corrected to represent the number of satellites brighter than the SAGA limit. Gray points are measured in this work. The green shaded area is the correlation between bulge-to-total ratio and the number of satellites from Javanmardi & Kroupa (2020), which is found for the total number of satellites for galaxies MW, M 101, M 33, M 31, Cen A, M 81, and M 94.

Another correlation that is linked to the relation between B/T ratio and the number of satellites is the correlation between the mass of the bulge and the number of satellites (Javanmardi et al. 2019). For each SAGA galaxy, we only have information on the number of satellites brighter than the SAGA limit and, on the other hand, Javanmardi et al. (2019) refers to the total number of satellites, and so we estimated the total number of satellites for SAGA galaxies in order to compare the results. The total number of satellites, if the relation between the number of satellites and B/T ratio (Javanmardi et al. 2019) holds, is: N = 33.6 (±6.5) B/T + 1.7 (±1.3), as given in Javanmardi & Kroupa (2020). This correlation, along with its uncertainties, is represented as the green shaded area in Fig. 2. We calculated the average number of satellites for SAGA galaxies according to this relation, and refer to it as the expected number of satellites. This expected number of satellites is marked with dark red crosses with black arrows for the error bars in Fig. 3. Overplotted are gray circles with errors referring to the number of satellites brighter than the SAGA limit (M_r, 0 < −12.3). However, blue asterisks mark six SAGA galaxies which, according to Fig. 2, have a complete census of satellites (they fall into the green shaded area), and for which a direct comparison can be made with the study of Javanmardi et al. (2019) and Javanmardi & Kroupa (2020). Unlike Fig. 2, where fainter satellites were excluded, the same galaxies (red points) in Fig. 3 harbor the total number of satellites. There is a positive trend indicating that galaxies with larger (more massive) bulges tend to harbor more satellites. Galaxies that have a complete census of satellites and thus are comparable (Fig. 3), labeled as blue asterisks from SAGA sample and red circles from Javanmardi et al. (2019), follow this positive correlation between the mass of the bulge and the (total) number of satellites. For these galaxies with a complete satellite census (13 in total), the Pearson’s linear correlation coefficient is 0.88.

Fig. 3. Bulge mass vs. expected number of satellites for SAGA galaxies (dark red crosses and blue asterisks) and the number of satellites brighter than the SAGA limit (gray circles). Overplotted are the galaxies from Javanmardi & Kroupa (2020) with their total number of satellites as red circles. Blue asterisks indicate SAGA galaxies that according to the relation between the B/T ratio and the number of satellites found in Javanmardi & Kroupa (2020) have a complete (or nearly) complete census of satellites.

4.1. Constraints from the empirical relations

For the other two relations we are investigating, namely number of satellites–specific angular momentum and number of satellites–galactic dynamical mass, we need to be certain of the quality of our measurements as well as that of the measurements we acquire from various resources.

First, we checked our stellar mass measurements using the empirical relation between the radius at which the surface brightness profile of the galaxy at 3.6 μm drops to 25.5 AB mag arcsec⁻² and its stellar mass (Muñoz-Mateos et al. 2015), the so-called stellar mass–size relation. For our sample of galaxies, the relation is given in Fig. 4. Stellar mass is plotted in gray circles with corresponding error bars, while the fit is given as a black solid line. The intention was not to make another correlation, but simply to find outliers. There are three 2-sigma outliers and we shall exclude them on the basis of inaccurate mass measurements. For NGC 5347 and NGC 5602, the composite fit (Fig. B.1) departs from the radial profile at the 25.5 AB mag arcsec⁻² level. This could be why these two galaxies fall off the relation. We cannot identify why galaxy NGC 4454 deviates from the relation, because the fit follows the radial profile. However, this galaxy along with NGC 5347 falls off the baryonic Tully–Fisher relation, as can be seen in Fig. 5. This is an additional reason to exclude these two galaxies from the sample.

Fig. 4. Stellar mass–size relation for SAGA galaxies: gray circles denote stellar mass with their corresponding errors. Stellar mass is correlated with the radius at which surface brightness at 3.6 μm drops to 25.5 AB mag arcsec−2 (R3.6). Galaxies labeled with their NGC numbers are 2-sigma outliers from the fitted relation (black line).

Fig. 5. Baryonic Tully–Fisher relation from McGaugh (2005): the solid line is simply Eq. (5) from McGaugh (2005), not the fit; 1-sigma outliers are labeled with their NGC names, omitting the NGC prefix, with the exception of five-digit number 68743, which corresponds to the PGC prefix. Other points are not labeled for clarity.

Additional parameters of huge importance to our further analysis are the rotation velocity and galaxy distance. Both parameters are equally important in deriving the dynamical mass and total specific angular momentum. The rotation velocity is measured from the width of the HI spectral line (Wmx_ave taken from EDD database), which is corrected for inclination as explained in detail in the previous section. The radius of the HI disk is estimated from ℳ_HI, which depends on the distance squared (Eq. (2)). Another galactic property of interest to us is the total specific angular momentum, for which again the rotation velocity and an accurate distance measurement are needed. We find no information on the HI line width (and/or rotation velocity) for three of the galaxies of our sample: NGC 7079, NGC 6909 and NGC 7029. To check both rotation velocity and distance accuracy, we employ the baryonic Tully–Fisher relation (McGaugh et al. 2000; McGaugh 2005). The baryonic Tully–Fisher relation correlates the total baryonic mass and the flat part of the rotation curve, which is often very close to the maximum rotation velocity that we have estimated using the width of HI line from the EDD database. This is the linear relation that reads: $\mathcal{M}/M_{\odot }=50\cdot V_{\text{rot}}^4$. Inserting this relation into Fig. 5 we calculated 1-sigma deviations and labeled the galaxies that deviate more than 1-sigma with their NGC names omitting the NGC prefix, and for the single case with five digits 68743 omitting the PGC prefix (PGC 68743). All these galaxies have either questionable distance measurements or inclination corrections. In either case, we exclude them from the sample for the subsequent analysis. Two of them are already 1-sigma outliers from the stellar mass–size relation (NGC 5347 and NGC 4454). In summary, we exclude seven more galaxies (apart from three edge-on galaxies) as outliers either from the stellar mass–size relation or the baryonic Tully–Fisher relation. To independently test our rotation velocity measurements, we searched the Hyperleda database through the EDD website for all SAGA galaxies and retrieved rotation velocity. Comparing our rotation velocities with values from the HyperLeda database⁶, we find three 1-sigma outliers: NGC 1309, NGC 2967, and NGC 4158. This is in line with our previous findings, except for the galaxy NGC 2967, which we decided to keep because it is one of the six galaxies that according to Javanmardi & Kroupa (2020) have a complete census of satellites.

4.2. Correlations with fundamental galactic properties

For dynamical mass measurements (Eq. (4)), both the rotation velocity and the radius of the HI disk are needed. For the galaxy NGC 5869, there are no measurements on its HI mass, and so we were not able to use Eq. (4) to infer its dynamical mass from the neutral gas. This galaxy is omitted from Fig. 6, where we plot dynamical mass versus the average number of satellites expected from the relation between the B/T ratio and the number of satellites found in Javanmardi & Kroupa (2020). The dark red crosses in Fig. 6 show the expected number of satellites given their B/T ratios if this latter relation holds for our sample of galaxies. Uncertainties are denoted with arrows. The number of bright satellites detected in the SAGA survey are presented as gray crosses which stand for error bars. The blue asterisks labeled with their NGC names are the only six galaxies in our sample that are considered to have a complete satellite census. Our analysis indicates that these six galaxies should stay in place, while all other galaxies (gray points) should migrate upwards toward red crosses, implying they have more satellites than detected. The inspected range of dynamical mass is rather narrow, but it shows no correlation, unless it is buried in the large scatter.

Fig. 6. Relation between dynamical mass and the number of satellites for SAGA galaxies, for which HI line measurements exist in the EDD database: gray circles with error bars denote the number of satellites detected up to the limit of Mr, 0 < −12.3; red crosses give the expected number of satellites if the relation from Javanmardi et al. (2019) holds. Error bars on the y-axis are simply the square root of the number of satellites. Dynamical mass is calculated from Eq. (4), with error bars corresponding to the errors on the rotation velocity only, meaning they are underestimated.

From each bulge component, we extracted information on the effective radius, Sersic index, and the ratio of minor to major axis (as ϵ = 1 − b/a) from the best fit model (Table A.1) to measure its specific angular momentum (Table 3) using Eqs. (7) and (8). Velocity dispersion was mainly taken from the SDSS SkyServer CasJobs Service⁷ where galaxies can be cross-matched against the SpecObjAll catalog with measured velocity dispersion based on their celestial coordinates. We find 18 matches within 0.5 arcmin radius from the galaxy center. However, the fiber is placed on the galaxy center for only ten of these matches. Three galaxies (NGC 2962, NGC 5869, and NGC 6278) were found in the catalog of Dabringhausen & Fellhauer (2016). We searched the HyperLeda database for the rest of the sample. As a result, data on the central velocity dispersion were gathered for all but six galaxies: NGC 7328, NGC 3689, NGC 5604, NGC 5792, NGC 4454, and NGC 5602. NGC 4454 and NGC 5602 are already excluded as outliers from empirical relations. Of the four remaining galaxies missing information on the central velocity dispersion, two of them (NGC 3689 and NGC 5604) have a stellar B/T ratio of less than 0.1, which makes the bulge angular momentum contribution negligible compared to the disk component, and these two galaxies are kept for further study. Galaxies NGC 7328 and NGC 5792 are excluded because the contribution of the bulge component to the total angular momentum cannot be neglected (stellar B/T ≥ 0.2), and thus the specific angular momentum of the disk cannot well represent the total angular momentum, contrary to the previous two cases of NGC 3689 and NGC 5604 galaxies. All the structural data measured from Spitzer’s images and those additional data retrieved from databases used for angular momentum estimation are listed in Table A.1.

From GALFIT exponential disk models, we extracted disk scale lengths and converted them into physical units (kpc) using galaxy distances from the SAGA dataset⁸ to measure the specific stellar angular momentum of the disks (Eq. (6)). With additional data on maximum rotation velocity and stellar B/T ratio, we measured specific angular momentum for 18 galaxies in total. The correlation between total specific angular momentum and the expected number of satellites is shown in Fig. 7. The expected number of satellites for each galaxy is marked with a dark red cross in Fig. 7 (Table 3), with the vertical arrow representing the error bar. The number of satellites detected in the SAGA survey is denoted with a gray circle with error bars holding the errors on both axes. The Pearson coefficient for the correlation is significant: −0.74. A single outlier in the Fig. 6 is the galaxy NGC 4158, which is also the main driver of the correlation with angular momentum (Fig. 7). Thus, deeper spectroscopic observations of fainter satellites would help us to decide whether this correlation truly holds. This galaxy should have 23 satellites fainter than the SAGA limit.

Fig. 7. Relation between the total specific angular momentum (Eq. (9)) and the number of satellites for SAGA galaxies, for which HI line measurements exist in the EDD database: gray circles with error bars denote the number of satellites detected up to the SAGA limit (Mr, 0 < −12.3); red crosses give the expected number of satellites if the relation from Javanmardi et al. (2019) holds.

5. Discussion and conclusions

All galaxies studied here are part of the SAGA survey as Milky Way analogs. We searched for previous and new correlations between the number of satellites and fundamental galactic properties to understand why some galaxies that are very similar regarding their stellar mass and environment have quite different numbers of satellites. Environment plays no role in their satellite distribution function because they are all selected as isolated galaxies.

The ΛCDM model predicts an excess of bright satellites compared to observations (the “too-big-to-fail” problem; Boylan-Kolchin et al. 2011), and predicts no significant correlations between the number of satellites and the mass of the bulge or the bulge-to-total ratio (Javanmardi et al. 2019; Javanmardi & Kroupa 2020) for the Milky Way-like galaxies. However, discrepancies were mainly obtained for Local Group galaxies. The necessity for a SAGA sample of galaxies comes from the concern that Local Group satellites are representative of typical galaxies in this mass range. The correlation between the bulge-to-total ratio (Javanmardi & Kroupa 2020) and the number of satellites holds for the galaxies with complete satellite census, which helped us identify six such galaxies in the SAGA sample. We confirmed a very strong correlation between the mass of the bulge and the number of satellites (Javanmardi et al. 2019) for these galaxies, including the sample of nearby galaxies for which the total number of satellites is known. We used the former correlation to predict the total number of satellites for SAGA galaxies based on their bulge-to-total ratio. This expected number of satellites makes the correlation between the mass of the bulge and the number of satellites even more significant. This is not unexpected, because the total number of satellites is predicted using another correlation including the mass of the bulge through the bulge-to-total ratio. What came as a surprise is the sign of the correlation, because bulges are assumed to form through accretion of or merger with the surrounding material, including satellites. It would be natural to assume that the larger the bulge, the more satellites have been accreted, but we see quite the opposite: the larger the bulge, the more satellites have survived. This reasoning is valid only if we assume that galaxies with similar baryonic mass evolving in the similar sparse environment have had initially comparable satellite populations.

A correlation that has been found in this work is the one between the total (expected) number of satellites and the total specific angular momentum. Angular momentum is the fundamental galactic property that is preserved during galaxy evolution. It can be redistributed via mergers, but preserved in general. SAGA galaxies with lower angular momentum have more satellites. Dynamical friction that can explain minor mergers would show the opposite behavior. The satellites should transfer their angular momentum to the host galaxy, slowing down and eventually spiraling in toward the central galaxy. The angular momentum of the host would consequently increase, because energy and momentum must be conserved. This scenario would imply that galaxies with more massive bulges should harbor less satellites and have higher angular momentum. Results from this and previous work (Javanmardi et al. 2019; Javanmardi & Kroupa 2020) suggest that the more abundant satellite populations have more massive bulges that are not built up by minor mergers. If the satellites are formed within the disk of baryonic and dark matter around galaxies, carrying away part of its angular momentum, then this would provide a natural explanation of the plane of satellites observed in several nearby galaxies (Pawlowski & Kroupa 2013; Müller et al. 2018). Deeper spectroscopic observations are needed to confirm these correlations because they critically depend on the number of satellites fainter than Leo I of the Milky Way.

Acknowledgments

We acknowledge the financial support of the Ministry of Education, Science and Technological Development of the Republic of Serbia through the contract No. 451-03-68/2022-14/200002. The authors thank the referee for the prompt report and constructive comments. This research made use of data from the SAGA Survey (sagasurvey.org). The SAGA Survey was supported by NSF collaborative grants AST-1517148 and AST-1517422 and by Heising-Simons Foundation grant 2019-1402. This research made use of IRAF and R programming languages. We acknowledge the use of the Extragalactic Distance Database (https://edd.ifa.hawaii.edu). We acknowledge the usage of the HyperLeda database (http://leda.univ-lyon1.fr). Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

References

Appendix A: Table of structural and kinematic properties of SAGA galaxies

Multicomponent modeling of SAGA galaxies with GALFIT code, from which the mass of the bulge is estimated, as well as the structural parameters of the bulge: Sersic index n, effective radius R_eff and ellipticity ϵ = 1 − b/a. In addition, velocity dispersions are retrieved from: the Sloan Digital Sky Survey catalog and the Hyperleda database. Radius of HI disk is estimated using Eq. 5 from gas mass M_HI that was retrieved from the Extragalactic Distance Database.

Appendix B: Surface brightness profiles of SAGA galaxies – decomposition into bulge and disk

Near-infrared images from the Spitzer Heritage Archive in 3.6 and 4.5 μm were retrieved and modeled with concentric ellipses using the IRAF ellipse task. These radial profiles are given in Fig. B.1 for both bands. Isophotal points from the ellipse task that deviate from the radial trend are artifacts from modeling without constraints on magnitude error and the number of valid points (unmasked pixels) inside ellipses. Constraints have not been imposed to recover the faintest parts of galaxies (the furthest points away from the center). They are not used for any analysis, and are given only as a representation of GALFIT modeling results in two dimensions. Images are also modeled with GALFIT code to extract individual components: bulge and disk. Apart from these two components, additional components were often required to obtain successful modeling: point spread function (PSF) for the nucleus and Sersic with low index n for bars and rings. Individual components are given for the 3.6 μm band only, for clarity. Composite models are generated as the sum of the fluxes of individual components: black dash dotted line in 3.6 μm, and gray dash dotted line in 4.5 μm. For comparison, isophotal contours from Muñoz-Mateos et al. (2013) are overplotted in the Fig. B.1. The results from this work agree well with the results from Muñoz-Mateos et al. (2013).

Fig. B.1. Surface brightness profiles of 32 SAGA galaxies in the NIR (3.6 and 4.5 μm). All galaxies are modeled with multiple components including bulge and disk, but also other components such as the nucleus, ring, and bar. Individual components are colored in the same manner in all the subfigures: PSF in orange, bulge in red, disk in green, and bar/ring in violet. Blue circles correspond to radial profile in 3.6 μm, while red circles correspond to 4.5 μm. Individual modeled components are presented for clarity only for 3.6 μm. The composite model is given for both 3.6 μm (black dash dotted line) and 4.5 μm (gray dash dotted line). The figure continues on the next page.

Fig. B.1. continued.

Fig. B.1. continued.

Fig. B.1. continued.

Three galaxies are edge-on and thus excluded from further study: NGC 4348, PGC 013646 and UGC 00903. Their radial profiles could not be modeled well with elliptical contours due to an extreme inclination. They are listed here only for completeness.

All Tables

All Figures

Fig. 1. Radial profile of the NGC 1015 galaxy in 3.6 and 4.5 μm-bands modeled with the IRAF ellipse task, denoted as blue and red points, respectively. Overplotted are orange crosses from Muñoz-Mateos et al. (2015) for the 3.6 μm-band. For clarity, only composite models are given for 3.6 μm (black) and 4.5 μm bands (gray), that correspond to the sum of all individual subcomponents. Individual subcomponents are presented only for the 3.6 μm band (point spread function, bulge and disk).

In the text

Fig. 2. Relation between the B/T ratio and the number of satellites brighter than Mr, 0 < −12.3 (SAGA limit). Red points are taken from Javanmardi & Kroupa (2020) and references therein. Their values on the y-axis are corrected to represent the number of satellites brighter than the SAGA limit. Gray points are measured in this work. The green shaded area is the correlation between bulge-to-total ratio and the number of satellites from Javanmardi & Kroupa (2020), which is found for the total number of satellites for galaxies MW, M 101, M 33, M 31, Cen A, M 81, and M 94.

In the text

Fig. 3. Bulge mass vs. expected number of satellites for SAGA galaxies (dark red crosses and blue asterisks) and the number of satellites brighter than the SAGA limit (gray circles). Overplotted are the galaxies from Javanmardi & Kroupa (2020) with their total number of satellites as red circles. Blue asterisks indicate SAGA galaxies that according to the relation between the B/T ratio and the number of satellites found in Javanmardi & Kroupa (2020) have a complete (or nearly) complete census of satellites.

In the text

	Fig. 4. Stellar mass–size relation for SAGA galaxies: gray circles denote stellar mass with their corresponding errors. Stellar mass is correlated with the radius at which surface brightness at 3.6 μm drops to 25.5 AB mag arcsec−2 (R3.6). Galaxies labeled with their NGC numbers are 2-sigma outliers from the fitted relation (black line).
In the text

	Fig. 5. Baryonic Tully–Fisher relation from McGaugh (2005): the solid line is simply Eq. (5) from McGaugh (2005), not the fit; 1-sigma outliers are labeled with their NGC names, omitting the NGC prefix, with the exception of five-digit number 68743, which corresponds to the PGC prefix. Other points are not labeled for clarity.
In the text

Fig. 6. Relation between dynamical mass and the number of satellites for SAGA galaxies, for which HI line measurements exist in the EDD database: gray circles with error bars denote the number of satellites detected up to the limit of Mr, 0 < −12.3; red crosses give the expected number of satellites if the relation from Javanmardi et al. (2019) holds. Error bars on the y-axis are simply the square root of the number of satellites. Dynamical mass is calculated from Eq. (4), with error bars corresponding to the errors on the rotation velocity only, meaning they are underestimated.

In the text

	Fig. 7. Relation between the total specific angular momentum (Eq. (9)) and the number of satellites for SAGA galaxies, for which HI line measurements exist in the EDD database: gray circles with error bars denote the number of satellites detected up to the SAGA limit (Mr, 0 < −12.3); red crosses give the expected number of satellites if the relation from Javanmardi et al. (2019) holds.
In the text

Fig. B.1. Surface brightness profiles of 32 SAGA galaxies in the NIR (3.6 and 4.5 μm). All galaxies are modeled with multiple components including bulge and disk, but also other components such as the nucleus, ring, and bar. Individual components are colored in the same manner in all the subfigures: PSF in orange, bulge in red, disk in green, and bar/ring in violet. Blue circles correspond to radial profile in 3.6 μm, while red circles correspond to 4.5 μm. Individual modeled components are presented for clarity only for 3.6 μm. The composite model is given for both 3.6 μm (black dash dotted line) and 4.5 μm (gray dash dotted line). The figure continues on the next page.

In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text