© The Authors 2026

1. Introduction

The stochastic nature of galaxy growth within the Λ cold dark matter (ΛCDM) framework favours the formation of dynamically hot central structures, such as classical bulges. These components are typically thought to arise either from early gravitational collapse or through the hierarchical merging of smaller progenitors at high redshift (De Lucia et al. 2011; Brooks & Christensen 2016; Costantin et al. 2021, 2022).

Within this context, the existence of massive disc galaxies (M_* > 10¹⁰ M_⊙) with little or no bulge remains a puzzle (Kautsch 2009; Barazza et al. 2009; Kormendy et al. 2010; Di Teodoro et al. 2023). The specific formation pathways of these bulgeless disc galaxies shape not only their internal structure and morphology (e.g. Sales et al. 2012; Du et al. 2021; Rosas-Guevara et al. 2026) but also influence the properties of their satellite systems.

In the dwarf regime, galaxy formation models and subgrid physics prescriptions are often calibrated using the observed properties of Local Group dwarfs, particularly the satellite galaxies of the Milky Way (MW). Thanks to its proximity, the MW offers a uniquely detailed view of satellite populations (e.g. McConnachie 2012; Simon 2019). Historically, the abundance of satellites on the MW has been at the heart of one of the most well-known challenges to the ΛCDM model, the missing satellites problem (Klypin et al. 1999; Moore et al. 1999). Over time, however, this issue has become better understood. Factors such as subhalo occupation statistics, galaxy formation efficiency, and inclusion of baryonic processes in simulations have all contributed to alleviating, or even solving, the tension originally raised by the discrepancy between DM only simulations and early observations of the MW satellite population (see Sales et al. 2022; Jung et al. 2024 and references there in).

Since satellites and their hosts co-evolve over extended periods (Wetzel et al. 2015; Simpson et al. 2018), satellite populations may carry imprints of their host’s formation history, morphology, and evolution. Their global properties, such as luminosity functions and radial distributions, are sensitive to the physics of reionisation (e.g. Bullock et al. 2000; Okamoto et al. 2008), environmental effects such as tidal and ram-pressure stripping (e.g. Mayer et al. 2006; Emerick et al. 2016), and the structural properties of the host itself (e.g. Garrison-Kimmel et al. 2017).

It also remains unclear whether the MW satellite system is representative of galaxies of similar mass in the Universe or whether it is an outlier. Cosmological simulations suggest that the MW satellite population may be atypical in several respects: The presence of two massive satellites is uncommon (Haslbauer et al. 2024; Buch et al. 2024), and the radial distribution of its satellites is more centrally concentrated than expected (e.g. Carlsten et al. 2020). The MW also exhibits an unusually high quenched fraction, nearly all of its satellites are quenched, exceeding that observed around comparable hosts in surveys such as SAGA and ELVES (Karunakaran et al. 2023; Geha et al. 2024), possibly reflecting the MW’s relatively quiet merger history (Kruijssen et al. 2019). High-resolution simulations similarly show that MW–like halos can display elevated quenched fractions (e.g. Engler et al. 2023; Rodríguez-Cardoso et al. 2025), although such results fall within a broader halo-to-halo scatter. Both simulations and observations consistently indicate substantial system-to-system variation in satellite population properties, even at a fixed host mass (Font et al. 2021; Engler et al. 2023; Mao et al. 2024), underscoring the need for large statistical samples when assessing how typical the MW truly is.

Over the past years, large observational programs such as ELVES (Carlsten et al. 2022), SAGA (Geha et al. 2017), xSAGA (Wu et al. 2022), and others (e.g. Zaritsky et al. 2024) have significantly expanded the sample of galaxies with known satellite populations. These efforts have collectively enabled statistical analyses over samples of tens to hundreds of nearby galaxies, improving our understanding of satellite demographics beyond the Local Group. However, most of these works have focused on general MW-mass hosts, often without isolating the impact of bulge morphology.

The Bulgeless Evolution And the Rise of Discs¹ (BEARD) survey (Méndez-Abreu in prep.), is an observational program targeting a volume-limited (< 40 Mpc) sample of 75 galaxies, of which 54 are massive bulgeless spirals (M_* ≥ 10¹⁰ M_⊙). Its main scientific aim is to recover the past merger history and mass growth of bulgeless galaxies. Bulgeless galaxies are defined here as systems with a bulge-to-disc mass ratio < 0.08 and a bulge-to-total mass ratio < 0.1, following the photometric decomposition of Zarattini (in prep.).

Understanding whether such galaxies pose a challenge to ΛCDM theories is therefore of central importance. The BEARD project is a multi-facility survey, and among its different observations, we carried out a deep photometric survey covering 100 kpc radii around our bulgeless galaxies. The BEARD large field of view and photometric depth enable simultaneous characterisation of both the central galaxies (Marrero-de la Rosa et al. 2026) and their satellite populations (Choque-Challapa in prep.)².

In this work, we use the cosmological hydrodynamical simulation TNG50-1 (Pillepich et al. 2018a; Nelson et al. 2018) to study the satellite populations of massive bulgeless galaxies, and to investigate how their properties correlate with host morphology. This paper is structured as follows: In Sect. 2 we describe the simulations and galaxy identification. In Sect. 3 we present the sample of host bulgeless galaxies, we describe the construction of a comparison control sample and the satellite selection criteria. In Sect. 4 we explore the properties of the satellite galaxies at a population level, comparing across the different samples. Finally, in Sect. 5 we provide a summary of our findings together with our conclusions.

2. Methodology

2.1. The TNG50 simulation

For this study, we used the TNG50-1 simulation (Pillepich et al. 2019; Nelson et al. 2019a) from the IllustrisTNG project (Pillepich et al. 2018a; Marinacci et al. 2018; Naiman et al. 2018; Springel et al. 2018; Nelson et al. 2018, 2019b). IllustrisTNG³ is a suite of cosmological magneto-hydrodynamic simulations that follow the evolution of structure formation in a ΛCDM Universe using the moving mesh code AREPO (Springel 2010). The simulations assume a cosmology consistent with results from the Planck Collaboration (Planck Collaboration XXVII 2016), specifically: h = 0.6774, Ω_m = 0.3089, Ω_Λ = 0.6911, and σ₈ = 0.8159. Initial conditions were generated at redshift z = 127 and evolved forward to z = 0.

Among the different realisations of the IllustrisTNG project, we focused on the highest-resolution run, TNG50-1 (Pillepich et al. 2019). This simulation follows the co-evolution of dark matter, gas, stars, and supermassive black holes within a (51.7 Mpc)³ comoving volume, using 2160³ dark matter particles and an equal number of initial gas resolution elements.

The dark matter particle mass is m_DM = 4.5 × 10⁵ M_⊙, while the mean baryonic mass resolution is m_baryon = 8.5 × 10⁴ M_⊙. Stellar particles inherit similar initial masses and lose a fraction of their mass over time due to stellar evolution.

The gravitational softening length for collisionless particles (dark matter and stars) is ϵ = 575 comoving pc until z = 1, after which it is fixed at ϵ = 288 pc in physical units. Gas cells have adaptive resolution, reaching sizes as small as ∼72 pc in dense, star-forming regions.

TNG50 incorporates detailed subgrid physics to model unresolved processes, including radiative cooling, chemical enrichment, and a multiphase interstellar medium with stochastic star formation (Pillepich et al. 2018b). Stellar feedback is implemented via kinetic galactic winds, while black hole feedback follows a dual-mode model: thermal at high accretion rates and kinetic at low accretion rates (Weinberger et al. 2017; Pillepich et al. 2018b). Thanks to its combination of high resolution and cosmological framework, TNG50-1 is particularly well suited for investigating the demographics and structural properties of satellite galaxies within realistic galactic environments.

2.2. Halo and subhalo identification

All halo and merger tree data used in this work were obtained from the publicly released TNG50-1 catalogues (Nelson et al. 2019b), where a complete description can be found. Here, we provide a brief summary.

Dark matter halos and subhalos in the TNG simulations were identified using the SUBFIND algorithm (Springel et al. 2001), applied to friends-of-friends (FoF) groups constructed from the dark matter distribution with a standard linking length of b = 0.2 times the mean inter-particle separation. The FOF algorithm identifies virialised structures (halos), while SUBFIND subsequently decomposes each halo into gravitationally self-bound substructures (subhalos), including both central and satellite galaxies.

Each subhalo is assigned a host FoF group, and baryonic components (gas, stars, black holes) are associated with the closest bound dark matter particle. The most massive subhalo in a given FoF group is considered the central, while the others are labelled as satellites. This inclusive halo-finding approach makes it possible to track galaxies within the potential wells of larger systems; however, we note that subhalo identification remains a complex challenge. Recent studies (e.g. Mansfield et al. 2024; Diemer et al. 2024; Forouhar Moreno et al. 2025) have shown that many halo finders can struggle to recover subhalos after pericentric passages, when there is substantial mass loss and blending with the central host density, making their detection more difficult.

To study the evolutionary history of halos and galaxies, we used merger trees constructed via the SUBLINK algorithm (Rodriguez-Gomez et al. 2015). SUBLINK links subhalos across snapshots by tracing the most bound particles, enabling the reconstruction of the formation and accretion history of both centrals and satellites.

3. Sample selection

3.1. Host galaxies

We are interested in exploring the properties of the satellite population of BEARD-like galaxies. Thus, from the TNG50-1 simulation we identified all the central subhalos whose stellar mass is between 1 × 10¹⁰ M_⊙ and 5 × 10¹¹ M_⊙. Close galaxy pairs would indeed affect our measurements of the satellite populations, thus we decided to include an isolation criterion: We excluded from the sample any host galaxy that has a companion within 500 kpc whose stellar mass is greater than 80% of the candidate’s. This selection resulted in a parent sample of 538 galaxies.

In addition, we imposed a morphological constraint to select only bulgeless galaxies. To this aim, we split the parent sample according to its bulge-to-disc mass ratio (B/D). We obtained the bulge-to-disc mass ratio from the kinematic decomposition of Zana et al. (2022). Here, we provide a brief description of their method, but we refer the interested reader to the original work for further details. The decomposition is based on analysing the binding energy and circularity (η ≡ j_z/j_circ(E) with j_z, the angular momentum perpendicular to the plane, and j_circ, the angular momentum of a circular orbit with energy E) distribution of stellar particles to classify them into five main structures: thin disc, thick disc, pseudo-bulge, bulge, and stellar halo. The classification begins by separating stellar particles into more and less bound components based on their binding energy. A Monte Carlo sampling technique is then applied to the most bound sample, assigning particles to the bulge and selecting those with η < 0 along with their symmetric positive counterparts⁴. The remaining particles are classified as part of the thin disc if their circularity is sufficiently high (η > 0.7) or as part of the pseudo-bulge. A similar approach is used for the less bound component to determine whether particles belong to the stellar halo, thin disc, or thick disc.

This procedure does not take into account the possible existence of a bar. The star particles forming a bar (if any) are likely assigned to the pseudo-bulge component, with some contamination in all the other ones.

This kinematic decomposition is not likely to be a one-to-one comparison with photometric decompositions (see the discussion in Cristiani et al. 2024) such as the one used in the BEARD survey (Zarattini in prep.). Thus, a hard cut in bulge-to-disc mass ratio equal to the one used in observational studies may not be adequate. In contrast we decide to get as our sample of BulgeLess (BEARD-Like, BL) galaxies those that lie in the first quartile of the bulge-to-disc mass ratio distribution (Fig. 1). Note that in this work the disc component includes thin and thick disc, while the bulge component only includes the classical bulge. With this approach, we effectively select the most extreme BL galaxies in a volume limited sample. Moreover, we select as BL galaxies those whose bulge mass is only as much as ≲6% per-cent of the disc mass, close to the 8% observed for our MW (Shen et al. 2010).

Fig. 1. Top panel: Distribution of host galaxy stellar mass. Bottom panel: Bulge-to-disc (B/D) mass ratio of host galaxies as a function of stellar mass. Indigo diamonds show the relation from Méndez-Abreu et al. (2021), obtained under the assumption B + D = T (see Section 3.1). Both panels display the full parent sample (black), BL sample (blue), the BD sample (orange), and, for clarity, one realisation of the CS (green). In the bottom panel, black horizontal dashed lines mark the first and fourth quantiles of the bulge-to-disc distribution, which define the division between BL and BD galaxies.

To understand the properties and particularities of the satellite population of BL galaxies, we need a comparison sample. For doing so, we decide to compare with bulge dominated (BD) galaxies (the fourth quartile in the bulge-to-disc mass ratio distribution, B/D ≳ 0.38).

However, as expected from previous works (e.g. Weinzirl et al. 2009; Méndez-Abreu et al. 2021), it is clear from Fig. 1 that there is a mass dependence of the bulge-to-disc mass ratio. In the original study from Méndez-Abreu et al. (2021), the authors report a positive correlation between the bulge-to-total mass ratio (B/T) and stellar mass using data from the CALIFA IFS survey (Sánchez et al. 2016). Here, in Fig. 1, we have transformed their data by assuming a two component galaxy (only bulge and disc, B + D = T), and thus the bulge-to-disc mass ratio can be expressed as B/D = (B/T)/(1 − (B/T)). This transformation is fully consistent with their formalism, as they assume a two-component galaxy model.

As the number of satellite galaxies is expected to scale with the host halo mass (e.g. Gao et al. 2004; Klypin et al. 2011), and therefore with stellar mass, the mass dependence of the bulge-to-disc mass ratio introduces a bias in the comparison of the satellite population properties of the BL and BD samples. In order to overcome this issue, we produced a third sample whose stellar mass distribution is statistically similar to the stellar mass distribution of the BL sample. We refer to this new sample as control sample (CS). In Fig. 2 we show a RGB composition of the stellar component of a random selection of host galaxies belonging to each of the three samples. In this mosaic the morphological difference between the samples and massive discs of the BL sample is clearly visualised.

Fig. 2. Face-on color-composite images of randomly selected galaxies from the samples used in this work. Images have been generated with the python module pynbody (Pontzen et al. 2013). The three channels are the stellar surface brightness in the i, v, and u Johnston filters. From top to bottom: BL, BD, and CS galaxies. From left to right, galaxies are sorted by stellar mass in increasing order. Each image has a width of 50 kpc. Stellar mass and B/D is given for each galaxy.

The CS has been defined as follows. First, we estimate the probability distribution of the (log₁₀) stellar mass of the BL (P_BL(log₁₀M_*)) and BD (P_BD(log₁₀M_*)) samples via kernel density estimation (solid lines in top panel of Fig. 1). Then we sample, with replacement, N galaxies (form the BD sample) with a sampling probability proportional to P_BL(log₁₀M_*)/P_BD(log₁₀M_*). In this case N = 135, in order to have the same number of host galaxies in all the samples of the analysis. This has the drawback that the same galaxy may be selected more than once. However, the method ensures that too massive galaxies, which are overrepresented in the BD (compared to the BL) sample, have low sampling probabilities, while those galaxies that are under-represented (i.e. low mass galaxies) have high sampling probabilities.

The resulting CS has a (log₁₀) stellar mass distribution indistinguishable from the BL sample. In Fig. 1 we show one realisation of the CS compared to the BL. In this work, except when stated otherwise, confidence regions regarding the CS were obtained by repeated measurements over 1000 different realisations of the CS. This prevents us from being biased towards any particular realisation of the CS. A Kolmogorov–Smirnov (KS) test between the two samples returns a median maximum distance of 0.12 with a median p-value of 0.3, not being able to distinguish between the two populations, as expected as the two samples are equal by construction. Despite the fact that our selection does not, a priori, enforce matching total halo masses, we have verified that the halo mass distributions of the BL and CS samples are statistically consistent (see Appendix A), and that their mass accretion histories are also nearly identical (see Appendix B). This confirms that the results presented here are not influenced by differences in the underlying host properties.

3.2. Satellite galaxies

We identified as satellite every subhalo located within a 3D aperture of 300 kpc from the host. This aperture is approximately the virial radius of a MW-like halo⁵ and has been extensively used in the literature for satellite identification in both simulations (e.g. Samuel et al. 2020; Engler et al. 2021; Font et al. 2022) and observations (e.g. Ruiz et al. 2015; Mao et al. 2021, 2024; Carlsten et al. 2022). Thus, for consistency with previous works, we decided to fix this aperture in the search of satellite galaxies.

We also imposed some refinement criteria by removing all the subhalos that do not have cosmological origin, i.e. those when traced back in time were formed inside another halo (this has been done according to the SubhaloFlag in the public data release). Furthermore, we removed all the subhalos with fewer than 50 star particles, to ensure we are not affected by resolution issues. Finally, we mimicked BEARD’s expected magnitude limit by imposing a further cut in the r band absolute magnitude removing all satellites dimmer than −11 mag in the r band (Choque-Challapa in prep.). This selection implies an effective lower stellar mass cut of the satellite galaxies of ∼2 × 10⁶ M_⊙⁶. As we are only focused in massive well resolved satellites we do not include any correction from orphan subhalos.

This procedure lead to a total of 515 satellites in the BL samples, 1051 in the BD sample, and an average of 506 ± 50 in the CS (Table 1). The strong difference in the satellite number counts between the BL and BD samples highlights the need of the mass resampling described above (Sect. 3.1).

For comparison throughout this work, we use the catalogue of Local Group dwarfs from Pace (2025)⁷. We applied selection criteria consistent with those used for the TNG50 satellite galaxies: We included only galaxies located within 300 kpc of their host and with stellar masses greater than 2 × 10⁶ M_⊙.

4. Results and discussion

4.1. Satellite number counts

The abundance of satellite galaxies remains sensitive to the physics implemented in simulations. While modern cosmological simulations reproduce the observed number of MW satellites (e.g. Engler et al. 2021) and the missing satellite problem has largely diminished (Sales et al. 2022); some authors even consider it solved (Jung et al. 2024). Nevertheless, it remains important to examine how satellite number counts scale with host galaxy properties. Here, we focus on the relation between the number of satellites and host stellar mass (Fig. 3).

Fig. 3. Number of satellite galaxies as function of the host’s galaxy stellar mass colour-coded by the host galaxy bulge-to-total mass ratio (B/T). Note that the y-axis switches from linear to logarithmic scaling at Nsat = 1. The black line shows a power law fit, while the shaded region indicates the 95% confidence interval, estimated from the NB regression. For comparison we include the values of the MW and Andromeda galaxies. The stellar masses have been obtained from Licquia & Newman (2015) and Sick et al. (2015) respectively. The satellite number counts were obtained from the Pace (2025) database, mimicking our satellite selection criteria (see Section 3.2).

We modelled the number of satellites, N_sat, as a power-law function of the host stellar mass, M_*, using the following form:

$$\begin{array}{c}\hfill ln\left(\mathbb{E}\left[N_{\mathrm{sat}}\right]\right)={\beta }_0+{\beta }_{1}ln\left(\frac{M_{\ast }}{{10}^{10}{\mathrm{M}}_{\odot }}\right),\end{array}$$(1)

where, β₀ and β₁ are regression coefficients estimated through a generalised linear model (GLM) using a log link function (e.g. McCullagh & Nelder 1989; Dobson & Barnett 2018).

Motivated by previous studies (e.g. Ruiz et al. 2015; Müller & Crosby 2023) suggesting that satellite abundance may also correlate with bulge mass fraction ($B/T=M_{*}^{\mathrm{Bulge}}/M_{*}$), we tested a second regression model in which the slope of the power-law depends linearly on B/T:

$$\begin{array}{c}\hfill ln\left(\mathbb{E}\left[N_{\mathrm{sat}}\right]\right)={\beta }_0+({\beta }_1+{\beta }_{2}B/T)ln\left(\frac{M_{\ast }}{{10}^{10}{\mathrm{M}}_{\odot }}\right).\end{array}$$(2)

To perform the regression, we used both Poisson and Negative Binomial (NB) likelihoods. Although Poisson regression is the standard approach for count data, our residuals showed evidence of overdispersion, where the variance exceeds the mean. In such cases, the Poisson model may underestimate the parameter uncertainties, leading to biased inferences. To address this, we adopted the NB likelihood, which includes an extra dispersion parameter, α⁸, which allows for a more accurate estimation of uncertainties.

We fit a series of Bayesian regression models using the bambi (Capretto et al. 2022) Python package. The priors and regression results are summarised in Table 2. For model comparison, we employ the expected log pointwise predictive density (ELPD), estimated via Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO, Vehtari et al. 2017), as implemented in ArviZ (Kumar et al. 2019). The NB models clearly outperform their Poisson counterparts. Among the NB models, the simpler version (i.e. without dependence on B/T) is slightly preferred.

The presence of over-dispersion is expected as the number of satellite galaxies is more tightly correlated with its host dark matter halo mass (e.g. FoF group mass) than with the central stellar mass alone (see Jiang & van den Bosch 2017, for an analysis on the non-Poissonity of the subhalo number counts), introducing an additional source of scatter in the relation. If we had ignored over-dispersion (i.e. using Poisson likelihood only), we would have inferred a statistically significant dependence on B/T.

Interestingly, our results are in contrast with those of Müller & Crosby (2023), who studied TNG100 and found a significant dependence on B/T. The difference between these two works it is likely related to the resolution effects and simulation volume. In our analysis we use TNG50 because of its higher resolution compared with TNG100. As discussed in Sect. 3.2, this higher resolution, together with our restrictive selection criteria, ensures that our results are not affected by the resolution–driven incompleteness that may influence the sample in Müller & Crosby (2023). However, the simulation volume of TNG50 is eight times smaller, implying that the most massive host galaxies may be under-represented in our sample. This possible incompleteness at the high–mass end is important to bear in mind, since the strongest correlation between bulge mass fraction and satellite number counts reported by Müller & Crosby (2023) happens to be stronger in the most massive hosts.

In summary, we do not find a significant correlation between host morphology and satellite abundance. The apparent difference between the BL and BD samples is largely explained by their different stellar masses; once we account for mass effects (i.e. when comparing the BL sample with the CS), the satellite number counts become comparable (Table 1).

4.2. Luminosities and colours

Galaxy magnitudes were computed by adding the contribution of all star particles within each subhalo. In IllustrisTNG, each star particle represents a single stellar population with a Chabrier IMF (Chabrier 2003); its initial metallicity is inherited from the gas from which it formed. Stellar photometry is based on BC03 SSP models (Bruzual & Charlot 2003) with Padova isochrones. The default stellar particle magnitudes are intrinsic (dust-free) and reported in standard photometric bands (e.g. SDSS r band).

In Fig. 4 we show a colour–magnitude diagram of the satellite galaxies of this study. The majority of our satellite galaxies present redder colours, as expected from the large quenched fractions of satellite galaxies in this mass range, which is of the order of ∼80 percent (e.g. Donnari et al. 2021; Mao et al. 2024). No clear differences are observed in the overall colour distribution between the different satellite populations; however, an excess of faint satellites is noticeable in the BL population.

Fig. 4. Colour–magnitude diagram of satellite galaxies. Satellites from BL galaxies are in blue, from BD galaxies in orange, and from the CS in green. Each point represents a satellite galaxy in the sample. While the colour distributions are similar across populations, BL satellites show a noticeable excess at faint magnitudes. For visualisation purposes only one realisation of the CS is shown. Top and right panels show the 1D distribution of r band magnitude and colour respectively. Distributions are shown via KDE (thick lines) and normalised histograms (thin lines).

We decided to characterise the differences between satellite populations by estimating their luminosity functions (LFs). For this purpose, we adopted a Schechter (1976) function parametrization:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill n(M;\varphi ,M^{\mathrm{knee}},\alpha )& =0.4ln(10)\varphi {10}^{0.4(M^{\mathrm{knee}}-M)(\alpha +1)}\hfill \\ \hfill & \times exp[-{10}^{0.4(M^{\mathrm{knee}}-M)}],\hfill \end{array}\end{array}$$(3)

where M is the satellite absolute magnitude in the r band, M^knee is the characteristic magnitude marking the knee of the luminosity function, ϕ is the normalisation, and α is the faint-end slope.

To fit the model to our data, we employed the extended maximum likelihood method introduced in Andreon et al. (2005), which accounts for variable completeness limits across host galaxies. The corresponding log-likelihood is given by

$$\begin{array}{c}\hfill ln\mathcal{L}(\{M_{\mathit{ij}}\}\mid \varphi ,M^{\mathrm{knee}},\alpha )={\displaystyle \sum _{j=1}^{N_{\text{hosts}}}}({\displaystyle \sum _{i=1}^{N_{\text{sat}}^j}}lnn(M_{\mathit{ij}};\varphi ,M^{\mathrm{knee}},\alpha )-s_j),\end{array}$$(4)

with the correction term s_j defined as:

$$\begin{array}{c}\hfill s_j={\displaystyle {\int }_{M_{\mathrm{low}}^j}^{M_{\mathrm{high}}^j}}n(M;\varphi ,M^{\mathrm{knee}},\alpha )\mathrm{d}M.\end{array}$$(5)

Here, we fixed the lower magnitude limit to M_low = −11 mag, consistent with our satellite selection threshold. The upper limit is fixed at M_high = −100 mag, following the prescription in Negri et al. (2022). We verified that this choice does not significantly impact the derived parameters, since bright satellites are rare⁹.

We adopted wide uninformative priors for the Schechter function parameters: the faint-end slope α, characteristic magnitude M^*, and normalisation log₁₀ϕ. Specifically, the priors are uniform over the ranges −5 < α < 5, −30 < M^knee/mag < 0, and −25 < log₁₀ϕ/dex < 15. The posterior distribution was sampled using Markov chain Monte Carlo (MCMC) method, as implemented in the emcee Python package (Foreman-Mackey et al. 2013), which employs an affine-invariant ensemble sampler.

The resulting luminosity functions are shown in Fig. 5. Posterior estimates of the parameters are summarised in Table 3. The BD sample clearly shows a higher satellite abundance, as seen in the offset of its luminosity function relative to the BL and CS samples. This emphasises the importance of the CS in distinguishing mass effects from morphological ones.

Fig. 5. Stacked satellite LFs of the BD (orange), BL (blue) and one realisation of the CS. Shaded regions indicate the 68% and 95% confidence intervals. Error bars in the data points indicate the 1σ Poisson uncertainty (y-axis) and the bin width (x-axis). Luminosity functions have not been fitted to these data-points but to the full distribution (see Sect. 4.2 for details). For visualisation purposes for the CS we show only one random realisation. Inset plot: Mean and scatter (1σ and 2σ) of the marginalised distribution of the faint end slope for the BL and BD samples. For the CS sample we show the distribution of the maximum a posteriori estimates of the faint end slope α over the 1000 realisations.

A visual inspection of the luminosity functions suggests that the CS exhibits an overabundance of bright satellites and a deficit of faint satellites compared to the BL sample. This trend is more clearly reflected in the marginalised posterior distributions of the faint-end slope α (inset panel of Fig. 5), where we find a ∼2.8σ discrepancy: the BL sample shows a steeper slope, indicating an excess of faint satellites. Of the three samples, the BL has a faint-end slope consistent within 1σ with the MW best-fitting value of ${\alpha }_{\mathrm{MW}}=-1.{34}_{-0.03}^{+0.04}$ (Nadler et al. 2019), while the other two (BD and CS) are slightly offset at the ∼1.5σ level (see Table 3).

The remaining Schechter parameters (M^* and ϕ) are consistent between the BL and CS samples. This implies that both samples host a similar total number of satellites, in agreement with the results presented in Sect. 4.1, where no significant dependence of satellite abundance on B/T was found. However, while the normalisation and characteristic magnitude remain similar, the satellites in BL systems are systematically fainter.

Moreover, when examining the stellar mass ratio of the most massive satellite to its host (Fig. 6), we see a clear dichotomy that complements with the luminosity function differences. There is a ∼1 dex difference between the most likely max-mass ratio between the BL and CS samples. Not only that but BL systems almost never host a satellite more massive than ∼10% of its host stellar mass, while both the BD and CS samples exhibit a significant heavy tail in the most massive satellite mass fraction distribution.

Fig. 6. Distribution of the ratio between the stellar mass of the most massive satellite and its host galaxy for the BL sample (blue), BD sample (orange), and CS (green). Error bars indicate 1σ Poisson uncertainty. Solid lines are KDE estimations. The green bands encompass the 68% and 95% confidence regions of the estimated PDFs over the 1000 realisations of the CS.

4.3. Satellite spatial distribution

Another parameter of interest when characterising satellite populations is their spatial distribution, both radial and angular. In Fig. 7, we show the radial distribution of the three samples analysed in this work.

Fig. 7. Stacked empirical cumulative distribution function of the radial position of the satellites galaxies of the BL (blue), BD (orange), and CS (green) samples. The CS shows the 1σ and 2σ confidence intervals estimated over 1000 re-samplings as green bands. The dotted lines indicate the 1σ scatter due to halo-to-halo variance. Black line shows the distribution of MW satellites whose positions and errors have been obtained from Pace (2025). The shadow black region indicates the 1σ confidence region due to observational uncertainties in the distance measurements. Dashed lines indicate the best fit to an NFW radial distribution. The resulting concentrations c300 kpc from the fit are shown in the inserted panel in the bottom right of the figure. The MLE of c300 kpc for the BL (blue) and BD (orange) samples are shown as a vertical line, the shaded bands indicate 68% and 95% confidence region. Regarding the CS we show the distribution of the 1000 MLE estimations of c300kpc as a green PDF.

The BL sample shows a more centrally concentrated satellite distribution compared to both the BD sample and CS. On average, half of its satellites lie within ∼120 kpc, indicating a more compact spatial configuration relative to the other two samples (see Table 4).

For comparison, we also include in Fig. 7 the radial distribution of MW satellites from the Local Volume Database (Pace 2025). The shaded regions represent the 68% confidence interval, obtained by propagating the uncertainties in the reported distances of the MW satellites.

It is clear that the MW satellites are significantly more centrally concentrated than any of the simulated samples considered in this work. This discrepancy is a well-known issue (e.g. Yniguez et al. 2014; Carlsten et al. 2020) and it seems to be related to the recent infall of the large magellanic cloud (LMC, Patel et al. 2024) due to both: its gravitational influence, and the group infall of satellite galaxies associated with the LMC.

We do find, for each sample, a large halo-to-halo variance (i.e. variations among individual host galaxies within the same sample) which can account for variations of ∼40% in the estimation of the half distance (see Table 4). We indicate such halo-to-halo variations in Fig. 7 with dotted lines. Even with this large halo-to-halo scatter, within the inner 100 kpc, the MW satellite distribution remains mildly (≳1σ) elevated relative to the expected halo-to-halo scatter, but aligns more closely with the BL sample.

Alternatively, we decided to describe the satellite radial distribution with an NFW profile (Navarro et al. 1997), as it allows for a better comparison with the dark matter halo concentration. The fit to the satellite distribution was performed by maximising the log-likelihood defined through the normalised differential enclosed mass profile:

$$\begin{array}{c}\hfill \mathcal{L}(\{r_{\mathit{ij}}\}\mid c)={\displaystyle \prod _{j=1}^{N_{\mathrm{hosts}}}}{\displaystyle \prod _{i=1}^{N_{\mathrm{sat}}^j}}\frac{1}{M_{200}(c)}\frac{\mathrm{d}M(r_{\mathit{ij}},c)}{\mathrm{d}r},\end{array}$$(6)

where M₂₀₀(c) is the total mass enclosed within R₂₀₀ for an NFW profile with concentration c, and

$$\begin{array}{c}\hfill \frac{\mathrm{d}M(r,c)}{\mathrm{d}r}=4\pi r^2{\rho }_{\mathrm{NFW}}(r,c),\end{array}$$(7)

with ρ_NFW(r, c) being the standard NFW density profile. For the radial fits, we fixed R₂₀₀ ≡ R_300 kpc = 300 kpc which corresponds to the maximal aperture within which we identify satellite galaxies. The concentration is defined in terms of the radius where the logarithmic slope of the density profile equals −2 (r₋₂), such that: c ≡ c_300 kpc = R₂₀₀/r₋₂ = 300 kpc/r₋₂.

It is worth noting that, while the NFW fit broadly captures the general shape of the radial distribution of satellites, in the BL sample it systematically deviates at both small (r ≲ 20 kpc) and large radii (r ≳ 150 kpc). Specifically, the NFW model overpredicts the number of satellites in the central regions and underpredicts it in the outskirts. This indicates that satellites in the BL population are slightly less centrally concentrated than a standard NFW profile would suggest. However, these deviations are minimal, amounting to fewer than 5%.

A direct fit to the satellite radial distributions (assuming an NFW distribution) yields a concentration of log₁₀c_300kpc = 0.79 ± 0.05 for the BL sample, and log₁₀c_300kpc = 0.48 ± 0.06 for the CS sample. When compared to the concentration of the host galaxy’s dark matter halo, these values are lower by about ∼0.5 dex, a difference similar to the ∼0.3 dex offset reported by McDonough & Brainerd (2022) using the TNG100 simulation.

This difference in the spatial concentration of satellites does not appear to mirror the concentration of their host dark matter halos. The average halo concentrations¹⁰ are nearly identical for the BL and CS samples, with values of ${log}_{10}{c}_{200,c}^{\mathrm{DM}}=1.11\pm 0.01$ and ${log}_{10}{c}_{200,c}^{\mathrm{DM}}=1.10\pm 0.02$, respectively. Instead, the difference seems to reflect variations in the angular momentum of the host systems. This applies to both the baryonic component, manifested in their different morphologies, and the dark matter component, with the average spin parameter¹¹λ_DM being roughly 0.2 dex lower in the CS sample than in the BL sample (see Appendix C). This is consistent with the results of Rosas-Guevara et al. (2026), who, in a study of BD galaxies, found that such systems preferentially undergo mergers that coherently contribute to their angular momentum. We explore the angular distribution of the satellite systems in the following sections.

4.4. Orbital alignment

In Fig. 8 we explore the angular distribution of satellite orbits. To quantify the degree of alignment, we computed the angle between each satellite orbital angular momentum (L_Sat) and the host galaxy stellar angular momentum vector (L), measured within one stellar half-mass radius:

$$\begin{array}{c}\hfill \theta =arccos\frac{\mathit{L}·{\mathit{L}}_{\mathrm{Sat}}}{|\mathit{L}||{\mathit{L}}_{\mathrm{Sat}}|}.\end{array}$$(8)

Fig. 8. Distribution of satellite orbital alignment with respect to host stellar angular momentum computed within 1 half mass radius at z = 0. BL in blue, BD in orange, and CS in green. Error bars in the histogram indicate 1σ Poisson scatter. Solid thick lines show a best fit to a beta distribution. Green band indicates the 1σ and 2σ interval of the best fitted beta distribution to the 1000 realisations of the CS. The black dashed line indicates what we would expect from an isotropic distribution of satellite orbits We include for comparison the orbital orientation of MW satellites (indigo dots) obtained from the Pace (2025) dataset. Note that for visual purposes we have added random vertical offset to each satellite.

The dashed black line in Fig. 8 represents the expectation for a completely isotropic angular distribution. We find that satellites in the BL sample exhibit a clear excess at small angles, indicating preferential alignment with the host galaxy angular momentum. This suggests that satellites in BL systems tend to orbit in a flattened (but thick), co-rotating structure, a planar configuration aligned with the disc plane.

As a reference, we include in Fig. 8 the orbital alignments of the Milky Way (MW) satellites. The strong clustering of the MW satellites reflects the well-known plane of satellites identified in the MW system. Interestingly, the orientation of the MW satellite plane (61.0 ± 1.5 deg in our selected sub-sample of MW satellites) approximately coincides with the most likely alignment in the BL sample¹². However, we emphasise that connecting the distribution of orbital alignments to individual planes of satellites is not straightforward and requires a dedicated, case-by-case analysis of plane identification, which is beyond the scope of this work and deferred to future study.

This angular alignment complements the findings from the radial distribution (Sect. 4.3), where BL galaxies showed more concentrated satellite distributions. Together, these results suggest that the satellite systems of BL hosts have dynamically colder and more coherent orbital structures than those of the CS or BD samples.

4.5. Satellite infall times and alignment evolution

We estimate the satellite infall time as the time when the distance between the satellite and the host galaxy is smaller than the R₂₀₀(z) of its host. Since we use a fixed physical aperture at z = 0 to identify satellites, in some cases the satellite candidate has never crossed R₂₀₀(z). We call these pre-infall satellites.

The fraction of pre-infall satellites is: 1.6 ± 0.6%, 0.9 ± 0.2 and 2.3 ± 0.8% in the BL, BD, and CS, respectively. The values for the BL and CS populations are consistent within 1σ, but the BD sample shows a noticeably smaller fraction. This reduction in the BD sample is primarily due to an aperture selection bias driven by the host-mass distribution: hosts in BD are more massive and have larger R₂₀₀, so our fixed physical aperture at z = 0 contains fewer galaxies that lie outside the host R₂₀₀.

The three populations under study show qualitatively similar distributions of infall times (Fig. 9). All exhibit a bimodal distribution, with one population of long-lived satellites that had its infall approximately at z ∼ 1, and a second population of satellites that have infallen recently (at cosmic times > 9 − 10 Gyr) or are on their first approach. This bimodal distribution is expected (see Simpson et al. 2018) and is due to a selection effect. If we took into account splash-back and already merged galaxies, this bimodality would be expected to vanish.

Fig. 9. Satellite infall time distribution. Infall times are estimated as the first moment the satellite lies within R200 of its host. BL in blue, BD in orange, and CS in green. Histogram error bars indicate 1σ Poisson scatter. Solid lines are KDE estimates of their respective distributions. Shaded bands in the CS indicate to the 1σ and 2σ regions over the 1000 realisations.

The long-lived population of satellites infalls slightly later in the BL sample. There is also an enhancement in the recently infallen satellites in the CS sample.

We repeated the same experiment as in Sect. 4.4 and we explored the orbital alignment of the satellites but at infall time (Fig. 10 top panel). Interestingly, despite the observed angular alignment of satellites in the BL population at z = 0, we do not find strong evidence of an anisotropic accretion of the satellite galaxies. This is in line with previous studies that found that anisotropic accretion is less prominent in less massive (M_halo ≲ 10¹³ M_⊙) halos (e.g. Tenneti et al. 2021).

Fig. 10. Top panel: Distribution of satellite orbital alignment, at infall time, with respect to host stellar angular momentum computed within 1 half mass radius. BL in blue, BD in orange and CS in green. Error bars in the histogram indicate 1σ Poisson scatter. Solid thick lines show the best fitting beta distribution. Green band indicates the 1σ and 2σ interval of the best fitted beta distribution to the 1000 realisations of the CS. Black dashed line indicates what we would expect from an isotropic distribution of satellite orbits. Bottom panel: Evolution of the mean (median) orbital alignment relative to the satellite infall time as a solid (dashed) line. We only show bins with more than 100 satellites. The shaded band indicates the 1σ uncertainty estimated via bootstrap resampling. Black horizontal thin dashed line at 90 deg indicates the expected value of an isotropic distribution.

Our results suggest that the coherent satellite configuration observed in BL hosts may instead arise from long-term orbital evolution within the host potential. We show this in the bottom panel of Fig. 10, where we follow the average orbital orientation of the satellites since their infall time, estimated as:

$$\begin{array}{c}\hfill \langle \theta \rangle (\mathrm{\Delta }t)=\arg (\frac{1}{N_{\mathrm{sat}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{sat}}}}exp\left(i{\theta }_j(\mathrm{\Delta }t)\right)),\end{array}$$(9)

where the summation is done over all the satellites of the sample, i is the imaginary unit, θ_j is the orbital alignment of the j-th satellite at a time Δt after its infall. The arg function extracts the angle of the resulting complex number. The uncertainties were estimated via bootstrap resampling.

As shown in Fig. 10, before infall the orbits are, on average, at 90 deg with respect to the host angular momentum, reflecting the isotropy of the distribution. There is a small systematic deviation (before infall) in the BL sample of the order of ∼5 deg which likely reflects the larger angular momentum environment this massive disc galaxies live in (see Di Cintio et al. 2019; Rosas-Guevara et al. 2026).

After infall, the satellite orbits of all three samples tend to steadily align with their host angular momentum at a similar rate: −2.5 ± 0.2 deg Gyr⁻¹ (BL), −2.1 ± 0.13deg Gyr⁻¹ (BD) and −2.36 ± 0.17 deg Gyr⁻¹ (CS). This persistent deviation from isotropy is primarily driven by satellites gradually aligning their orbits with the principal axes of their hosts. Such behaviour is consistent with classical studies of satellite–host interactions, where tidal torques, dynamical friction, and resonant angular-momentum exchange drive orbital decay, and inclination changes, progressively coupling satellite orbits to the symmetry plane of the central potential (e.g. Quinn & Goodman 1986; Peñarrubia et al. 2004; Read et al. 2008; Ogiya & Burkert 2016).

At later times (several gigayears after infall), the BD and CS samples depart from the initial steady alignment: in the BD sample the alignment rate weakens down to −0.64 ± 0.12 deg Gyr⁻¹, whereas in the CS sample the process effectively stalls or even reverses, reaching 0.9 ± 0.2 deg Gyr⁻¹. In contrast, the BL satellites maintain their initial alignment rate over the full interval with an alignment rate of −2.5 ± 0.3 deg Gyr⁻¹.

We interpret this divergence as evidence that, although all systems undergo a common secular alignment after infall, this evolution is still susceptible to dynamical perturbations. The infall of massive companions, more frequent in BD and CS hosts, can disturb the coherence of the satellite system and temporarily weaken or stall the alignment, whereas the more weakly perturbed BL hosts preserve the steady trend over several gigayears. Massive accretion events therefore act to modulate an intrinsically morphology-dependent secular process rather than fully determine the observed differences. We assess this interpretation explicitly in Appendix D, where we separate hosts according to the presence or absence of massive infall and compare their alignment histories.

5. Summary and conclusions

In this study we investigated the satellite populations of BL galaxies using the high-resolution TNG50-1 cosmological simulation. Our analysis focused on BL systems that are morphological analogues of the Milky Way and the primary targets of the BEARD multi-facility survey. Within TNG50-1 simulations we identified MW and BEARD-like analogues by selecting central galaxies with stellar mass between 1 × 10¹⁰ M_⊙ and 5 × 10¹¹ M_⊙ that fall in the first quartile of the bulge-to-disc mass-ratio distribution. For comparison, we also constructed two reference samples: a BD sample drawn from the upper quartile of the bulge-to-disc distribution, and a stellar-mass–matched bulge dominated CS (Fig. 1). By systematically comparing the satellite populations of BL with both BD and CS, we uncovered a set of robust trends:

• Satellite abundance: We do not find a statistically significant dependence of the number of satellite galaxies with bulge mass fraction. The difference of satellite number counts between the BD and BL is due to differences in the mass distribution (Fig. 3).

• Luminosity functions and stellar mass: The satellite luminosity function of BL hosts displays a significantly (2.8σ) steeper faint-end slope compared to the CS, indicating an overabundance of faint satellites in BL galaxies (Fig. 5) and an excess of more massive satellites in the CS (Fig. 6). Satellite luminosity functions of BL hosts are in better agreement (< 1σ) with current MW estimations, compared to the CS satellite galaxies (∼1.5σ).

• Spatial distribution: Satellites of BL galaxies are more centrally concentrated than those around BD and CS hosts, with d_half ∼ 120 kpc and log₁₀c_300 kpc ≈ 0.8. Despite this, there is substantial halo-to-halo variance, which can change d_half by roughly forty per cent. When compared to the Milky Way, whose satellites have d_half = 49–84 kpc and log₁₀c_300 kpc ≈ 1.2, all simulated samples remain less centrally concentrated, although the BL population provides the closest match.

• Orbital alignment: A preferential alignment between the orbital angular momentum of satellites and the stellar angular momentum of the host galaxy is found in BL systems, indicating a preference for co-rotating satellite configurations, in contrast to the more isotropic orbital configurations seen in BD and CS hosts (Fig. 8). Satellite orbital alignment follows a secular post-infall evolution that tends to align satellite orbits with the host disc (Fig. 10). At later times, the alignment weakens or stalls in CS hosts, and to a lesser extent in BD hosts, while remaining steady in BL systems. Our extended analysis (Appendix D) indicates that the infall of massive satellite companions can perturb this secular evolution, contributing to the observed slowdown of the alignment in affected systems. In BL hosts, where massive satellites are rare, the secular alignment proceeds uninterrupted over several gigayears. These results indicate that the efficiency of orbital alignment is intrinsically morphology-dependent, with massive accretion acting to modulate rather than fully determine the trends.

Overall, our findings suggest that the morphological state of a galaxy, particularly the absence of a bulge, is closely linked to the properties of its satellite system. The steeper luminosity functions, enhanced angular alignment, and spatial concentration of satellites in BL hosts point to a more dynamically coherent and possibly quiescent evolutionary history (in agreement with the recent work of Rosas-Guevara et al. 2026).

These results provide a useful theoretical benchmark for future and current observational campaigns (e.g. BEARD, ELVES, SAGA), and highlight the potential of satellite populations as tracers of the formation history and internal structure of their host galaxies.

Acknowledgments

We thank the anonymous referee for their constructive comments, which improved the clarity of the paper. SCB and JMA acknowledges the support of the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (MCIN/AEI/10.13039/501100011033) under grant nos. PID2021-128131NB-I00 and CNS2022-135482 and the European Regional Development Fund (ERDF) ‘A way of making Europe’ and the ‘NextGenerationEU/PRTR’. AdLC acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through RYC2022-035838-I and PID2021-128131NB-I00 (CoBEARD project). EAG acknowledges support from the Agencia Espacial de Investigación del Ministerio de Ciencia e Innovación (AEI-MICIN) and the European Social Fund (ESF+) through a FPI grant PRE2020-096361. MCC acknowledges the support of AC3, a project funded by the European Union’s Horizon Europe Research and Innovation programme under grant agreement No 101093129. MCC acknowledges financial support from the Spanish Ministerio de Ciencia, Innovación y Universidades (MCIU) under the grant PID2021-123417OB-I00 and PID2022-138621NB-I00. SZ acknowledges the financial support provided by the Governments of Spain and Aragón through their general budgets and the Fondo de Inversiones de Teruel, the Aragonese Government through the Research Group E16_23R, and the Spanish Ministry of Science and Innovation and the European Union – NextGenerationEU through the Recovery and Resilience Facility project ICTS-MRR-2021-03-CEFCA. NCC acknowledges a support grant from the Joint Committee ESO-Government of Chile (ORP 028/2020) and the support from FONDECYT Postdoctorado 2024, project number 3240528. EMC and AP acknowledge the support by the Italian Ministry for Education University and Research (MUR) grant PRIN 2022 2022383WFT “SUNRISE” (CUP C53D23000850006) and Padua University grants DOR 2022-2024. FP acknowledges support from the Horizon Europe research and innovation programme under the Maria Skłodowska-Curie grant “TraNSLate” No 101108180, and from the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (MCIN/AEI/10.13039/501100011033) under grant (PID2021-128131NB-I00) and the European Regional Development Fund (ERDF) “A way of making Europe”. ADC kindly thanks the Spanish Ministerio de Ciencia, Innovacion y Univerasidades through grant CNS2023-144669, programa Consolidacion Investigadora. J.R. acknowledges financial support from the Spanish Ministry of Science and Innovation through the project PID2022-138896NB-C55 and from the Plan Propio de Investigación 2025 submodalidad 2.3 of the University of Córdoba. The author(s) wish to acknowledge the contribution of the IAC High-Performance Computing support team and hardware facilities to the results of this research. This work has extensively used the following software: numpy (Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), emcee (Foreman-Mackey et al. 2013), ArviZ (Kumar et al. 2019), bambi (Capretto et al. 2022), pynbody (Pontzen et al. 2013) and numba (Lam et al. 2015).

References

Appendix A: Halo mass distribution

As discussed in Section 4.1, halo mass is the primary factor governing both the number of satellites and their properties. In this section, we therefore present the halo mass distributions of the samples explored in this work (Fig. A.1). BD galaxies tend to reside in more massive halos, by roughly ∼0.5 dex, compared to the BL sample. However, the halo mass distributions of the galaxies of the CS and BL samples are statistically indistinguishable. A Kolmogorov–Smirnov test yields a median p-value of 1.380, confirming that the two distributions are fully consistent with being drawn from the same parent population.

Fig. A.1. Top panel: Distribution of host galaxy stellar mass. Bottom panel: Bulge-to-disc (B/D) mass ratio of host galaxies as a function of total halo mass. Both panels display the full parent sample (black), the BL sample (blue), the BD sample (orange), and, for clarity, one realisation of the CS sample (green). In the bottom panel, black horizontal dashed lines mark the first and fourth quantiles of the bulge-to-disc distribution, which define the division between BL and BD galaxies.

Appendix B: Mass Accretion Histories

In this appendix, we present the mass accretion histories (MAHs) of the galaxy samples analysed in this work. The MAHs include both baryons and dark matter and are normalised to the final mass of each galaxy (Fig. B.1).

Fig. B.1. Mass accretion histories (baryons plus dark matter), normalised to the final mass. Shaded bands indicate the halo-to-halo scatter, enclosing the 16th and 84th percentiles of the MAHs for each galaxy. Dashed lines are included to delineate the same region and aid visualisation in cases where the shaded bands overlap.

BD hosts, on average, form slightly more slowly than galaxies from the BL and CS samples, reaching half of their stellar mass approximately 0.4 Gyr later (see Table B.1). Nevertheless, there is substantial variability due to the stochastic nature of individual MAHs, with the 1σ scatter in t₅₀ being ∼2.3, Gyr for BD, ∼2.5 Gyr for CS, and ∼1.9 Gyr for BL galaxies.

Although the scatter in MAHs is larger in the CS sample than in the BL sample, the average evolution of galaxies in these two samples is nearly identical. Despite this, the growth of the stellar component (not shown) differs between the samples, with the CS and BD galaxies exhibiting faster stellar mass growth, in agreement with the findings of Rosas-Guevara et al. (2026). This indicates that differences in the assembly histories of host galaxies cannot account for the observed variations in the properties of their satellite populations.

Appendix C: Host Halo DM concentrations and spin

Dark matter halo concentrations, $c_{200,c}^{\mathrm{DM}}$, were obtained from Anbajagane et al. (2022), where they were estimated by fitting an NFW profile to the dark matter density profile of each subhalo out to R_200, c, i.e. the radius enclosing a density 200 times the critical density of the Universe.

For the spin parameter, we adopt the definition of Bullock et al. (2001):

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{dm}}=\frac{j_{200,c}}{\sqrt{2}{V}_{200,c}{R}_{200,c}},\end{array}$$(C.1)

where j_200, c is the specific dark matter halo angular momentum, $V_{200,c}=\sqrt{\frac{G{M}_{200,c}}{R_{200,c}}}$, G the gravitational constant, R_200, c is the virial radius, and M_200, c is the dark matter halo mass. For consistency with the concentration definition, all these quantities were measured within an aperture of R_200, c.

Fig. C.1 shows a corner plot of the distributions and correlations between the dark matter halo concentration, $c_{200,c}^{\mathrm{DM}}$, and spin parameter, λ_dm, for the BL, BD, and CS. The expected distribution of dark matter halo spin from Bullock et al. (2001) is also indicated for reference.

Fig. C.1. Corner plot showing the joint and marginal distributions of dark matter halo concentration, $c_{200,c}^{\mathrm{DM}}$, and spin parameter, λ, for the BL (blue), BD (orange) and CS (green). Diagonal panels show the normalised marginal distribution. Shaded bands indicate 1σ and 2σ confidence regions over the 1000 realisations of the CS. Concentrations are obtained via an NFW fit to the DM density profile from Anbajagane et al. (2022) The 2D contours indicate the 0.5σ, 1σ, 1.5σ and 2σ confidence regions. Black line dashed in bottom right panel indicates the expected spin parameter λdm distribution from Bullock et al. (2001).

As expected from the concentration-mass relation (Dutton & Macciò 2014) the BD sample shows smaller concentrations than both BL and CS.

CS and BL samples have approximately equal average concentrations, with $log{c}_{200,c}^{\mathrm{DM}}$ −1.10 ± 0.02 and −1.114 ± 0.008 respectively. While the average are similar, the scatter is slightly larger in the CS, 0.22 ± 0.03 compared to the ∼0.097 ± 0.009 found in the BL. The uncertainties in the average and scatter are bootstrap estimates.

Regarding the DM halo spins, there is clear systematic difference between the BL galaxies and the CS. Having the CS (and the BD galaxies), extremely low values, even below the expectations from Bullock et al. (2001).

In summary, DM halos of BL and CS galaxies are radially similar, but show systematic differences in their angular momentum.

Appendix D: The effect of the massive infall in orbital alignment

In this section we examine the impact of massive infalling satellites on the evolution of the orbital alignment of satellite systems. To this end, we analyse three complementary selections:

• satellites in hosts that experience a massive satellite infall (Massive Accretion);

• satellites in hosts without such an event (No Massive Accretion);

• a cleaned sample in which the orbital evolution of each satellite is followed only up to the moment when a massive infall occurs in its host, if applicable.

If perturbations induced by massive satellites are the dominant driver of alignment changes, we expect the No Massive Accretion sample to follow the intrinsic secular alignment reported in Section 4.4, whereas in the Massive Accretion sample a slowdown or stalling of the alignment should appear a few gigayears after the infall of the massive companion. For the cleaned sample, we expect the evolution to track that of the full population until the time when the effects of the massive infall begins to be present in the Massive Accretion systems.

Figure D.1 shows the alignment evolution of these three selections: Massive Accretion (green), No Massive Accretion (red), and cleaned (orange); compared to the full sample (blue). Here, a massive infall is defined as a satellite–to–host mass ratio greater than > 0.1 at the time when the satellite reaches its maximum total mass. We explored higher thresholds, which yield qualitatively similar trends, although the statistics of the Massive Accretion sample naturally decrease as the threshold increases.

Fig. D.1. Mean orbital–alignment evolution for satellites in the BL (top), BD (middle), and CS (bottom) samples. Curves correspond to the full sample with all the satellites (blue), satellites in hosts that underwent massive satellite infalls (Massive Accretion; red), satellites in hosts without a massive infall (No Massive Accretion; green), and the cleaned sample tracking satellite orbits up to the moment of the massive infall if any (orange). Any satellite with a maximum mass ratio with the host galaxy > 0.1 is considered a massive infall. The comparison highlights systematic differences between the subsamples and allows for assessment of the impact of massive satellites on the alignment trends.

Clear systematic differences emerge between the Massive Accretion and No Massive Accretion populations: satellites in hosts with Massive Accretion tend to be, on average, less aligned than those in No Massive Accretion hosts. Interestingly, the divergence between the two samples begins roughly ∼1 Gyr before the infall time, suggesting that the perturbative influence of the massive satellite is already felt by the satellite system prior to the companion crossing the host’s virial radius.

The dependence of these effects on host morphology also reveals notable differences. In BL galaxies (top panel of Fig. D.1), the expected behaviour is most clearly visible: both the full sample and the Massive Accretion sample show a pronounced slowdown in their alignment evolution, whereas this feature is absent in the No Massive Accretion and cleaned samples.

The BD hosts (middle panel) show a more complex behaviour. Their average alignment is closer to isotropy (Section 4.4), making the slowdown in the Massive Accretion sample less prominent, though still present. Meanwhile, the No Massive Accretion and cleaned samples are almost indistinguishable and lie only a few degrees below the full-sample trend, so the significance of the difference is modest. A similar qualitative picture applies to the CS sample (bottom panel), but in this case the sample variance is too large relative to the expected signal to draw firm conclusions.

In summary, the differential slowdown associated with the infall of a massive satellite is strong and unambiguous in the BL hosts, and present but less statistically significant in the BD hosts due to their lower intrinsic alignment levels. For the CS sample, the effect remains inconclusive because of the large sampling variance. These findings also suggest that the underlying mechanism driving the secular alignment differs in strength (or potentially in nature) between BL and BD galaxies.

All Tables

All Figures

Fig. 1. Top panel: Distribution of host galaxy stellar mass. Bottom panel: Bulge-to-disc (B/D) mass ratio of host galaxies as a function of stellar mass. Indigo diamonds show the relation from Méndez-Abreu et al. (2021), obtained under the assumption B + D = T (see Section 3.1). Both panels display the full parent sample (black), BL sample (blue), the BD sample (orange), and, for clarity, one realisation of the CS (green). In the bottom panel, black horizontal dashed lines mark the first and fourth quantiles of the bulge-to-disc distribution, which define the division between BL and BD galaxies.

In the text

Fig. 2. Face-on color-composite images of randomly selected galaxies from the samples used in this work. Images have been generated with the python module pynbody (Pontzen et al. 2013). The three channels are the stellar surface brightness in the i, v, and u Johnston filters. From top to bottom: BL, BD, and CS galaxies. From left to right, galaxies are sorted by stellar mass in increasing order. Each image has a width of 50 kpc. Stellar mass and B/D is given for each galaxy.

In the text

Fig. 3. Number of satellite galaxies as function of the host’s galaxy stellar mass colour-coded by the host galaxy bulge-to-total mass ratio (B/T). Note that the y-axis switches from linear to logarithmic scaling at Nsat = 1. The black line shows a power law fit, while the shaded region indicates the 95% confidence interval, estimated from the NB regression. For comparison we include the values of the MW and Andromeda galaxies. The stellar masses have been obtained from Licquia & Newman (2015) and Sick et al. (2015) respectively. The satellite number counts were obtained from the Pace (2025) database, mimicking our satellite selection criteria (see Section 3.2).

In the text

Fig. 4. Colour–magnitude diagram of satellite galaxies. Satellites from BL galaxies are in blue, from BD galaxies in orange, and from the CS in green. Each point represents a satellite galaxy in the sample. While the colour distributions are similar across populations, BL satellites show a noticeable excess at faint magnitudes. For visualisation purposes only one realisation of the CS is shown. Top and right panels show the 1D distribution of r band magnitude and colour respectively. Distributions are shown via KDE (thick lines) and normalised histograms (thin lines).

In the text

Fig. 5. Stacked satellite LFs of the BD (orange), BL (blue) and one realisation of the CS. Shaded regions indicate the 68% and 95% confidence intervals. Error bars in the data points indicate the 1σ Poisson uncertainty (y-axis) and the bin width (x-axis). Luminosity functions have not been fitted to these data-points but to the full distribution (see Sect. 4.2 for details). For visualisation purposes for the CS we show only one random realisation. Inset plot: Mean and scatter (1σ and 2σ) of the marginalised distribution of the faint end slope for the BL and BD samples. For the CS sample we show the distribution of the maximum a posteriori estimates of the faint end slope α over the 1000 realisations.

In the text

	Fig. 6. Distribution of the ratio between the stellar mass of the most massive satellite and its host galaxy for the BL sample (blue), BD sample (orange), and CS (green). Error bars indicate 1σ Poisson uncertainty. Solid lines are KDE estimations. The green bands encompass the 68% and 95% confidence regions of the estimated PDFs over the 1000 realisations of the CS.
In the text

Fig. 7. Stacked empirical cumulative distribution function of the radial position of the satellites galaxies of the BL (blue), BD (orange), and CS (green) samples. The CS shows the 1σ and 2σ confidence intervals estimated over 1000 re-samplings as green bands. The dotted lines indicate the 1σ scatter due to halo-to-halo variance. Black line shows the distribution of MW satellites whose positions and errors have been obtained from Pace (2025). The shadow black region indicates the 1σ confidence region due to observational uncertainties in the distance measurements. Dashed lines indicate the best fit to an NFW radial distribution. The resulting concentrations c300 kpc from the fit are shown in the inserted panel in the bottom right of the figure. The MLE of c300 kpc for the BL (blue) and BD (orange) samples are shown as a vertical line, the shaded bands indicate 68% and 95% confidence region. Regarding the CS we show the distribution of the 1000 MLE estimations of c300kpc as a green PDF.

In the text

Fig. 8. Distribution of satellite orbital alignment with respect to host stellar angular momentum computed within 1 half mass radius at z = 0. BL in blue, BD in orange, and CS in green. Error bars in the histogram indicate 1σ Poisson scatter. Solid thick lines show a best fit to a beta distribution. Green band indicates the 1σ and 2σ interval of the best fitted beta distribution to the 1000 realisations of the CS. The black dashed line indicates what we would expect from an isotropic distribution of satellite orbits We include for comparison the orbital orientation of MW satellites (indigo dots) obtained from the Pace (2025) dataset. Note that for visual purposes we have added random vertical offset to each satellite.

In the text

	Fig. 9. Satellite infall time distribution. Infall times are estimated as the first moment the satellite lies within R200 of its host. BL in blue, BD in orange, and CS in green. Histogram error bars indicate 1σ Poisson scatter. Solid lines are KDE estimates of their respective distributions. Shaded bands in the CS indicate to the 1σ and 2σ regions over the 1000 realisations.
In the text

Fig. 10. Top panel: Distribution of satellite orbital alignment, at infall time, with respect to host stellar angular momentum computed within 1 half mass radius. BL in blue, BD in orange and CS in green. Error bars in the histogram indicate 1σ Poisson scatter. Solid thick lines show the best fitting beta distribution. Green band indicates the 1σ and 2σ interval of the best fitted beta distribution to the 1000 realisations of the CS. Black dashed line indicates what we would expect from an isotropic distribution of satellite orbits. Bottom panel: Evolution of the mean (median) orbital alignment relative to the satellite infall time as a solid (dashed) line. We only show bins with more than 100 satellites. The shaded band indicates the 1σ uncertainty estimated via bootstrap resampling. Black horizontal thin dashed line at 90 deg indicates the expected value of an isotropic distribution.

In the text

Fig. A.1. Top panel: Distribution of host galaxy stellar mass. Bottom panel: Bulge-to-disc (B/D) mass ratio of host galaxies as a function of total halo mass. Both panels display the full parent sample (black), the BL sample (blue), the BD sample (orange), and, for clarity, one realisation of the CS sample (green). In the bottom panel, black horizontal dashed lines mark the first and fourth quantiles of the bulge-to-disc distribution, which define the division between BL and BD galaxies.

In the text

	Fig. B.1. Mass accretion histories (baryons plus dark matter), normalised to the final mass. Shaded bands indicate the halo-to-halo scatter, enclosing the 16th and 84th percentiles of the MAHs for each galaxy. Dashed lines are included to delineate the same region and aid visualisation in cases where the shaded bands overlap.
In the text

Fig. C.1. Corner plot showing the joint and marginal distributions of dark matter halo concentration, $c_{200,c}^{\mathrm{DM}}$, and spin parameter, λ, for the BL (blue), BD (orange) and CS (green). Diagonal panels show the normalised marginal distribution. Shaded bands indicate 1σ and 2σ confidence regions over the 1000 realisations of the CS. Concentrations are obtained via an NFW fit to the DM density profile from Anbajagane et al. (2022) The 2D contours indicate the 0.5σ, 1σ, 1.5σ and 2σ confidence regions. Black line dashed in bottom right panel indicates the expected spin parameter λdm distribution from Bullock et al. (2001).

In the text

Fig. D.1. Mean orbital–alignment evolution for satellites in the BL (top), BD (middle), and CS (bottom) samples. Curves correspond to the full sample with all the satellites (blue), satellites in hosts that underwent massive satellite infalls (Massive Accretion; red), satellites in hosts without a massive infall (No Massive Accretion; green), and the cleaned sample tracking satellite orbits up to the moment of the massive infall if any (orange). Any satellite with a maximum mass ratio with the host galaxy > 0.1 is considered a massive infall. The comparison highlights systematic differences between the subsamples and allows for assessment of the impact of massive satellites on the alignment trends.

In the text