© The Authors 2026

1. Introduction

In the past decades, a significant part of the cosmological investigation has been dedicated to shedding light on the central topic of the origin and the growth of cosmic structures. These efforts led to the formulation of the so-called hierarchical model of structure formation (White & Rees 1978). According to this, structure formation is the result of the growth of primordial overdensities in the matter distribution that arise at the end of the inflationary era from initial quantum fluctuations in the early Universe. In this framework, gravitationally bound systems, such as galaxies and galaxy clusters, are the product of gravitational infall, the progressive merging of dark matter halos, and the catching up of baryons (White & Rees 1978; Bardeen et al. 1986; Tormen 1998).

Several physical processes contribute to the formation and assembly history of these complex systems. Among them, the long-range nonlinear gravitational interaction sets the dynamic of the overall matter distribution, while different radiative processes, including stellar emission and gas cooling, affect the baryonic matter in different phases. Finally, a series of feedback phenomena, such as supernova explosions (Efstathiou 2000; Smith et al. 2018), stellar winds (Springel & Hernquist 2003; Davé et al. 2011; Hopkins et al. 2012), and active galactic nucleus emission (Bower et al. 2006; Weinberger et al. 2018), reshape the matter distribution as well as the star formation rate (Scharré et al. 2024). Due to the large number of interactions and physical processes involved, an analytical treatment of structure formation is currently unfeasible.

For this reason, numerical simulations have been developed in a wide range of sizes and resolutions that target well-defined scientific cases. They offer important tools for investigating the statistical properties of cosmological systems and testing models that will be applied to observations and real datasets. In particular, for a given set of cosmological parameters, N-body simulations trace the evolution of the dark component alone, following the gravitational interaction of dark matter particles over cosmic time (Springel et al. 2001b; Springel 2005), while hydrodynamical simulations also include baryons, allowing us to understand the complex interplay between their different components, such as gas, stars, and black holes (Duffy et al. 2010; Vogelsberger et al. 2014a; Chisari et al. 2018; Springel et al. 2018; Vogelsberger et al. 2020).

The clustering of large-scale structure (LSS) is recognized today as a fundamental cosmological probe because during different periods of the cosmic expansion, it keeps the memory of the footprint left by the initial conditions of the Universe. In practice, this offers us the possibility to provide constraints on the main cosmological parameters, such as the matter density content of the Universe and the normalization of the primordial fluctuations (Hildebrandt et al. 2020; Asgari et al. 2021; Abbott et al. 2022; Amon et al. 2022; Lesci et al. 2022; Secco et al. 2022; Romanello et al. 2024). In addition, the clustering of LSS is typically combined with other cosmological probes, such as number counts of cosmic structures (To et al. 2021), or a simultaneous study of different two-point statistics (Heymans et al. 2021), which includes cosmic shear (Kilbinger 2015; Hildebrandt et al. 2020; Asgari et al. 2021) and galaxy–galaxy lensing, and can potentially highlight differences in the concordance cosmological paradigm.

Currently, the most widely acknowledged cosmological scenario is the Λ cold dark matter (ΛCDM) model, which assumes that dark matter particles are in a ‘cold’ version, namely in the form of very massive nonrelativistic candidates, such as weakly interacting massive particles, with a mass m_X > 1 GeV, or condensates of light axions, with m_X ≲ 10⁻³ eV (Strigari 2013). The ΛCDM model is consistent with observations at various scales that range from the typical intergalactic distances to the size of the cosmological horizon. It also provides robust explanations for observations of the cosmic microwave background, for the abundances of light elements as products of the primordial nucleosynthesis, and for the accelerating expansion of the Universe. Moreover, it reasonably describes the distribution and statistical properties of the LSS.

However, some possible tensions have been suggested by observations on the galactic and sub-galactic scales, of the order of kiloparsec (see Weinberg et al. 2015; Bullock & Boylan-Kolchin 2017; Del Popolo & Le Delliou 2017 for reviews). Among them, the ‘missing satellites’ (Klypin et al. 1999) and the ‘too-big-to-fail’ problems (Boylan-Kolchin et al. 2011), which are related to the fact that the abundance of satellites around the Milky Way does not match the number of galactic subhalos measured in N-body simulations, and the ‘cusp-core’ problem, which arises because the observed halo density profiles are less concentrated than the predicted cusped Navarro-Frenk-White profiles (NFW, Navarro et al. 1996). Although a recent galactic census has significantly alleviated the missing satellite problem (Kim et al. 2018), the issue remains open scientifically. Usually, baryonic phenomena are invoked as possible solutions (Brooks et al. 2013); for example, photoionization feedback to quench star formation in low-mass halos (Busha et al. 2010; Castellano et al. 2016; Yue et al. 2018), and supernova feedback and tidal stripping to expel gas from the galactic environment, resulting in a redistribution of matter that erases the internal cusps of density profile and in a reduction of the number of faint satellites (Brooks et al. 2013; Bullock & Boylan-Kolchin 2017).

Another perhaps more intriguing possibility is the exploration of the dark sector by investigating the macroscopic consequences of alternative cosmological scenarios that postulate the existence of warm dark matter (WDM, Bode et al. 2001; Viel et al. 2005; Schneider et al. 2012; Lovell et al. 2014) or self-interacting dark matter particles (SIDM, Spergel & Steinhardt 2000; Vogelsberger et al. 2012). In the first case, the growth of cosmic perturbations smaller than a characteristic free streaming length is suppressed by the thermal motion of collisionless dark matter particles. This naturally produces shallower density profiles and fewer low-mass halos (Menci et al. 2016). This observationally translates into a turnover at the faint end of the rest-frame galaxy UV luminosity function (Corasaniti et al. 2017) and in a general delay of galaxy formation (Dayal et al. 2015; Menci et al. 2018), with a possible effect on the timeline of the reionization process (Dayal et al. 2017; Carucci & Corasaniti 2019; Romanello et al. 2021). In the latter case, dark matter particles have a nonzero cross section of self-interaction, σ, which yields a non-negligible scattering probability that allows for a transfer and a redistribution of energy and momentum between different parts of dark matter halos (see Tulin & Yu 2018; Adhikari et al. 2025, for reviews). As a result, the inner region of the density profile tends to become cored and isothermal, whereas the outer part, in which the density and the interaction rate are lower, is well described by an NFW profile. In more realistic SIDM models, self-interaction is determined by a massive mediator and σ decreases with the relative velocities of the interacting particles (Vogelsberger et al. 2012; Correa 2021; Adhikari et al. 2025). Systems such as galaxy clusters, which have a large velocity dispersion, are therefore less affected than faint galaxies. As a consequence, SIDM affects on kiloparsec rather than on megaparsec scales, which agrees with what is expected from the observational evidence in the local Universe.

This publication is part of a series of works that make use of a suite of cosmological simulations (Despali et al. 2025) to study and characterize the properties of dark matter halos, such as their shapes (Giocoli et al. 2026), gas distribution (Zhang et al. 2026), concentrations and density profiles (Despali et al. 2026), in alternative dark matter models. The management of halo catalogs, the measurements of the relevant two-point statistics, the cosmological modeling, and the Bayesian inference presented in this study were performed with COSMOBOLOGNALIB¹ (Marulli et al. 2016), a set of free software libraries for C++ and Python.

2. The AIDA-TNG simulations

The AIDA-TNG sample (Alternative Dark Matter in the TNG universe, Despali et al. 2025) is a set of cosmological simulations designed to study the effect of alternative dark matter models (WDM and SIDM) on multiple scales. Boxes of 110.7 and 51.7 Mpc on a side, corresponding to 75 and 35 Mpc h⁻¹, are presented in Despali et al. (2025); we refer to them as “100/” and “50/”. Each box is simulated in up to six dark matter scenarios, as well as with and without baryons, and labeled according to its decreasing mass resolution: m_DM = 4.8 × 10⁷ M_⊙ h⁻¹ for 100/A, m_DM = 2.9 × 10⁶ M_⊙ h⁻¹ for 50/A, and 2.3 × 10⁷ M_⊙ h⁻¹ for 50/B, respectively. The simulations are listed in Table 1, with respective details and available cosmologies. The baryonic physics is included through the TNG galaxy formation model (Weinberger et al. 2017; Pillepich et al. 2018b).

The starting redshift of the simulations is z = 127, with initial conditions derived with the Zel’dovich approximation, as also in the IllustrisTNG simulations (Vogelsberger et al. 2014b; Pillepich et al. 2018a). The friends-of-friends (FOF) and SUBFIND algorithms (Springel et al. 2001a) have been used to identify dark matter halos and subhalos, respectively. First, the FOF group finder is run. We then refer to subhalos as sets of locally overdense, gravitationally bound particles that are hosted in each FOF detection. The position of central subhalos coincides with the minimum of the gravitational potential, that is, with the centers of the FOF halo (Pillepich et al. 2018a). Therefore, throughout this paper, we consider satellites and subhalos to be synonyms.

For the purpose of our work, we focus on dark matter-only 50/A and 50/B boxes at z = 2, z = 1, z = 0.5 and z = 0, in order to take advantage of the resolution and of the largest variety of dark matter models, while leaving further study on the baryonic influence to the future. The cosmological parameters of the simulations are taken from Planck Collaboration XXIV (2016): Ω_m = 0.3089, Ω_Λ = 0.6911, Ω_b = 0.0486, H₀ = 0.6774, and σ₈ = 0.8159.

The AIDA-TNG project explores the same cosmological volumes in different alternative dark matter scenarios. In particular, for WDM cosmologies, dark matter is considered a ‘thermal relic’, resulting from the freeze-out of a dark component in thermal equilibrium with the early Universe. From a physical point of view, the free streaming of dark matter particles determines a suppression in the power spectrum, P(k), which can be modeled through the transfer function, T(k) (Bode et al. 2001; Viel et al. 2005; Schneider et al. 2012),

$$\begin{array}{c}\hfill T(k)\equiv \sqrt{\frac{P_{\mathrm{WDM}}(k)}{P_{\mathrm{CDM}}(k)}}={[1+{(ϵk)}^{2\nu }]}^{-5/\nu },\end{array}$$(1)

where ν = 1.2, and

$$\begin{array}{c}\hfill ϵ=0.048{\left(\frac{{\mathrm{\Omega }}_X}{0.4}\right)}^{0.15}{\left(\frac{h}{0.65}\right)}^{1.3}{\left(\frac{\mathrm{keV}}{m_X}\right)}^{1.15}.\end{array}$$(2)

Here, Ω_X is the density parameter of the dark matter. We consider three thermal WDM cosmologies, with particle mass corresponding to m_X = 1 − 3 − 5 keV, respectively. Lighter models such as WDM1 and WDM3, with m_X = 1 − 3 keV, have been ruled out by Lyman-α, number counts, UV luminosity function observations (Viel et al. 2013; Menci et al. 2016; Corasaniti et al. 2017; Dekker et al. 2022; Villasenor et al. 2023) and are used only for a broader characterization of WDM abundances and clustering. Furthermore, in WDM cosmologies, the halo mass function experiences a suppression given by (Lovell et al. 2014)

$$\begin{array}{c}\hfill \frac{n_{\mathrm{WDM}}}{n_{\mathrm{CDM}}}={(1+\gamma \frac{M_{\mathrm{hm}}}{M})}^{\beta },\end{array}$$(3)

where γ = 2.7, β = −0.99 and M_hm is the half-mode mass, namely, the mass at which the transfer function in Eq. (1) is equal to 1/2.

We consider SIDM models with a constant cross section (SIDM1, with σ/m_X = 1 cm² g⁻¹) and, more realistically, with the velocity-dependent cross section presented in Correa (2021, vSIDM), which rapidly declines at the typical velocity dispersion of massive bound systems. The model fits to observations of central densities of dwarf spheroidal galaxies, and it is based on a self-interaction described by a particle mass m_X = 53.933 ± 9.815 GeV and a massive mediator in a Yukawa potential, with mass m_φ = 6.605 ± 0.435 MeV (Correa 2021).

In Fig. 1 we show the spatial distribution of dark matter halos in box 50/B at different redshifts and cosmologies. The nodes and filaments become progressively clumpier toward low z, as they collect matter from the background. In particular, WDM models show a clear progressive suppression in the number density of bound structures. SIDM cosmologies, on the other hand, are not represented, since, upon qualitative inspection, they do not exhibit a significant difference with respect to ΛCDM.

Fig. 1. From top to bottom: Spatial distribution of dark matter halos in box 50/B for the CDM, WDM3, and WDM1 cosmologies. The color bar shows the corresponding redshift in the range 0 < z < 2.

3. The halo occupation distribution model

Discrete structures, such as galaxies and dark matter halos, arise at the overdensity peaks of the Universe; therefore, they are biased tracers of the underlying matter distribution. On large scales, this reduces to a linear, scale-independent relation between the halo and the dark matter density contrast fields, δ_h(x) and δ(x), through a bias factor, b_h(z),

$$\begin{array}{c}\hfill {\delta }_{\mathrm{h}}(\mathit{x})=b_{\mathrm{h}}(z)\delta (\mathit{x}).\end{array}$$(4)

On smaller scales, however, it is challenging to accurately model the clustering of halos because nonlinear effects emerge.

In recent years, the halo occupation distribution (HOD) formalism has been widely used in number count and clustering studies because it represents a simple and powerful theoretical tool for approximating galaxy abundances, bias, and small-scale galaxy clustering (Zheng et al. 2005; Giocoli et al. 2010; Asgari et al. 2023). It is based on the halo model, namely the idea that all the matter in the Universe is collected into halos (Cooray & Sheth 2002). Under this hypothesis, the question comes down to understanding how dark matter halos cluster together, based on their classification as central or satellites, and how luminous tracers, such as galaxies, correspondingly populate them.

3.1. Number counts

According to the HOD model, the average number of halos is given by

$$\begin{array}{c}\hfill ⟨N_h⟩=⟨N_c⟩+⟨N_s⟩,\end{array}$$(5)

where ⟨N_c⟩ is the average number of central halos, while ⟨N_s⟩ represents the average number of satellites. Unlike galaxy clustering, where ⟨N_c⟩ is parameterized according to the details of the galaxy formation process and of the uncertainties in the luminosity-mass relation, for dark matter halos it is simply given by a Heaviside step function (Zehavi et al. 2005),

$$\begin{array}{c}\hfill ⟨N_c⟩=\mathrm{\Theta }(M-M_{\min }),\end{array}$$(6)

where M_min is the mass cut of the sample, below which the number of central and satellite halos is equal to zero. There must be one central halo in order to have satellite halos, and thus, for masses higher than M_min, where ⟨N_c⟩ = 1, we can write the number of subhalos as follows:

$$\begin{array}{c}\hfill \langle N_s\rangle =⟨N_c⟩{\left(\frac{M-M_0}{M_1}\right)}^{\alpha },\end{array}$$(7)

where M₁ and α are the normalization and slope of the power law, respectively, while the parameter M₀ is included in more general HOD models, as a satellite truncation mass, which is different from M_min. As we also find in our analysis, however, it is generally poorly constrained (Watson et al. 2009, 2012; Piscionere et al. 2015; Skibba et al. 2015).

In this sense, the halo occupation number, ⟨N_h⟩, gives us a first indication of the HOD parameters. Moreover, it allows us to test whether the distribution of satellites within halos of a given mass follows Poissonian statistics. In particular, we can investigate the average number of pairs associated with the same halo, which brings clustering information and is the base for the one-halo term of the two-point correlation function (Zheng et al. 2005),

$$\begin{array}{c}\hfill ⟨N_h(N_h-1)⟩=2⟨N_c{N}_s⟩+⟨N_s(N_s-1)⟩=⟨N_s^2⟩+⟨N_s⟩.\end{array}$$(8)

For a Poissonian distribution, this reduces to ⟨N_h(N_h − 1)⟩ = ⟨N_s⟩² + 2⟨N_s⟩. Therefore, for massive halos that contain satellites, the ratio $⟨N_h(N_h-1)⟩/{⟨N_h⟩}^2\stackrel{1}{\to }$, as we discuss in Sect. 4.1.

3.2. Radial distribution of subhalos

It is often assumed that the distribution of tracers follows the distribution of dark matter. For this reason, an NFW profile (Navarro et al. 1996) is generally adopted in clustering studies. However, it might not be sufficiently accurate to describe the small-scale distribution of subhalos, introducing systematic uncertainties in the parameter values of α and M₁ and in cosmology. In particular, we expect the concentration of satellite halos, c_h = f_c c_200c, to differ from that of dark matter, with f_c as a free parameter of our model². Therefore, following the approach of More et al. (2009), Watson et al. (2009, 2012) and Piscionere et al. (2015), we consider a generalized NFW (gNFW),

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_s}{{\left(\frac{r}{r_s}\right)}^{\gamma }{(1+\frac{r}{r_s})}^{3-\gamma }},\end{array}$$(9)

where γ corrects the slope of the profile and r_s = R_200c/(f_c c_200c) is the scale radius, which varies depending on the halo mass and concentration. To estimate the latter, we use the c_200c − M_200c relation from Duffy et al. (2008). Finally, ρ_s = M_200c/(4πr_s³∫₀^R_200c/f_cx^2 − γ/(1+x)^3 − γdx) is the characteristic density, where x ≡ r/r_s. For example, for a luminous red galaxy sample, the parameter γ is found to be an increasing function of luminosity, while f_c is a decreasing one (Watson et al. 2012). We retrieve the NFW profile when both γ and f_c are equal to one. A value of f_c smaller than one boosts r_s, eventually to a radius greater than R_200c, which means that the density profile keeps its inner slope −γ even at larger radial separations, behaving like a pure power law (Watson et al. 2009, 2012). Therefore, in order to explore the full radial profile, the upper limit of the integral is adjusted accordingly. Differences in γ and f_c among alternative cosmologies also account for variations of the profile shape and of the relative number of subhalos in the inner part of dark matter groups, respectively. As we show in the next section, this also impacts on the clustering properties of central and satellite halos.

3.3. Clustering

We use the HOD formalism to model the small-scale clustering of dark matter halos in cold and alternative dark matter scenarios. In particular, the halo model distinguishes two contributions to the correlation function,

$$\begin{array}{c}\hfill \xi (r)={\xi }_{1h}(r)+{\xi }_{2h}(r),\end{array}$$(10)

where ξ_1h(r) is the one-halo term, accounting for the clustering of matter or tracers hosted by the same halo, while ξ_2h(r) is the two-halo term, given by the pairs counted between different halos. The same formalism is equivalently written in Fourier space as

$$\begin{array}{c}\hfill P(k)=P_{1h}(k)+P_{2h}(k).\end{array}$$(11)

In Eq. (11), the former term dominates at small scales, in a full nonlinear range, that is, below the typical size of dark matter halos, while the latter describes large-scale clustering, where the impact of nonlinearities is smaller and P(k) = b_h²P_lin(k). The linear power spectrum, P_lin(k), is computed according to the approximation of Eisenstein & Hu (1998). This occurs for separations greater than the size of the largest dark matter halo.

Dealing with two different populations of tracers, in the one-halo term, we should separate the contribution of central-satellite and satellite-satellite pairs to the power spectrum, P_cs and P_ss, respectively,

$$\begin{array}{c}\hfill P_{1h}(k)=P_{\mathit{cs}}(k)+P_{\mathit{ss}}(k).\end{array}$$(12)

In particular,

$$\begin{array}{c}\hfill P_{\mathit{cs}}(k,z)=\frac{2}{n_h^2(z)}{\displaystyle {\int }_{M_{min}}^{M_{max}}}⟨N_c{N}_s⟩n(M,z){\stackrel{\sim }{u}}_s(k,M,z)\mathrm{d}M.\end{array}$$(13)

Here, the central and satellite occupation functions are independent, and therefore, ⟨N_cN_s⟩=⟨N_c⟩⟨N_s⟩. The normalized distribution of satellites within the halo, in Fourier space, is ${\stackrel{\sim }{u}}_s(k,M,z)$, directly derived by numerically transforming Eq. (9), whose parameters are estimated by fitting the stacked satellite radial density profile measured in Sect. 4.2. This measurement is possible only with an adequate number of subhalos, which requires a group mass M_200c > M₁. In addition, halos of mass lower than M₁ have R_200c smaller than the minimum separation at which we measured the two-point correlation function in Sect. 4.3. From a physical point of view this is equivalent to exclude the one-halo contribution for groups with a central halo only.

For the halo mass function, n(M, z), we adopt the functional form provided by Tinker et al. (2008). It can be integrated to compute the total halo number density,

$$\begin{array}{c}\hfill n_h(z)={\displaystyle {\int }_{M_{min}}^{M_{max}}}n(M,z)⟨N_h⟩\mathrm{d}M,\end{array}$$(14)

in which ⟨N_h⟩ is expressed through Eq. (5). The satellite-satellite power spectrum is

$$\begin{array}{c}\hfill P_{\mathit{ss}}(k,z)=\frac{1}{n_h^2(z)}{\displaystyle {\int }_{M_{\min }}^{M_{max}}}⟨N_s(N_s-1)⟩n(M,z){\stackrel{\sim }{u}}_s^2(k,M,z)\mathrm{d}M.\end{array}$$(15)

Finally, the two-halo power spectrum is given by

$$\begin{array}{cc}\hfill P_{2h}(k,z)& =P_{\mathrm{lin}}(k,z)[\frac{1}{n_h(z)}{\displaystyle {\int }_{M_{min}}^{M_{max}}}⟨N_h⟩n(M,z)b(M,z)\hfill \\ \hfill & \times {\stackrel{\sim }{u}}_s(k,M,z)\mathrm{d}M{]}^2(A+Bk+C{k}^2),\hfill \end{array}$$(16)

where b(M, z) is the linear bias provided by Tinker et al. (2010), while A, B and C form a scale-dependent nonlinear correction (see for example Smith et al. 2007), which reproduces the small-scale two-halo clustering of groups. Indeed, a more accurate description of the transition region where the one-halo and the two-halo contributions are almost equivalent, called quasi-linear regime, is traditionally more difficult to model (Mead et al. 2015), and requires to consider nonlinear bias (van den Bosch et al. 2013; Bhowmick et al. 2018) and exclusion effects (Tinker et al. 2005). Discrete tracers, such as galaxies and dark matter halos, cannot overlap; thus the distance r among them is always larger than the sum of their radii (Baldauf et al. 2013). This reduces the two-halo correlation function at small scales. Moreover, the two-halo term itself has a small-scale behavior that cannot be reproduced by the assumption of the model bias presented in Tinker et al. (2010), as it contains nonlinearities that, if not taken into account, would alter the two-point correlation function at the typical one-halo separation, ending up with systematic errors on the HOD parameters. Moreover, the impact of exclusion effects is limited by the abundance of groups with mass between 10⁸ M_⊙ h⁻¹ and 10¹⁰ M_⊙ h⁻¹, which have R_200c < 0.05 Mpc h⁻¹. This also allows us to count the two-halo pairs with spatial separations that are fully in the nonlinear regime. For this reason, for each redshift and cosmological model, we calibrate the nonlinear deviation from b_h²P_lin(k, z), summarized in the parameters A, B and C, by modeling the clustering of dark matter groups only, without considering the contribution of satellites.

On large scales, the shape of the correlation function is determined by ⟨N_h⟩, and essentially differs from the matter-matter correlation for a linear bias factor only. On the other hand, on small scales, the one-halo term is sensitive to ⟨N_h(N_h − 1)⟩ (Berlind & Weinberg 2002). A higher slope value, α or a lower normalization M₁, determines an increase in the number of satellites at high M. This largely boosts the number of satellite-satellite pairs in the one-halo term of the two-point correlation function. It also impacts the two-halo term, since higher masses, which are strongly biased, become more populated. In other words, for a given number of halos, we remove substructures from low-mass systems, redistributing them into more massive ones. Furthermore, higher values for γ and f_c, as in Eq. (9), steepen the satellite density profile, increasing the amplitude and the slope of the one-halo correlation function.

In Fig. 2 we show the different contributions to the two-point correlation function as a function of redshift, for a CDM cosmology and M > 10⁸ M_⊙ h⁻¹. In particular, on large scales, ξ_2h(r) is largely suppressed by the finite size of the box. This condition is called ‘integral constraint’ and requires to be forward-modeled in order to obtain unbiased results, as discussed in Sect. 4.3. On the other hand, on small scales, the predictions of the linear model are not accurate enough. In particular, at high redshift, the scale-dependent terms, B and C, are more relevant, while at z = 0 the constant A deviates more from 1. For example, at z = 2 we found A = 0.80 ± 0.04, B = 0.41 ± 0.02 and C = −0.010 ± 0.002, while at z = 0, A = 0.67 ± 0.03, B = 0.07 ± 0.01 and C = −0.004 ± 0.001. Although ξ_1h(r) always overcomes the linear ξ_2h(r) at r < 0.2 Mpc h⁻¹, it prevails on the nonlinear two-halo term only at z = 0. It is dominated by satellite-satellite pairs, meaning that our dark matter groups are largely populated by subhalos. On the other hand, the correlation function of the central-satellite pairs, ξ_cs, increases with the redshift in absolute terms and relative to its satellite-satellite counterpart, ξ_ss. Their ratio extends from a maximum of approximately 0.09 at z = 0 to 0.36 at z = 2. Given the same mass cut of 10⁸ M_⊙ h⁻¹, moving toward high z, we find a larger contribution of isolated or poorly populated halos. Conversely, when we only select subhalos with a mass greater than 10⁹ M_⊙ h⁻¹, the term ξ_cs becomes progressively more important and dominates ξ_ss for M > 10¹⁰ M_⊙ h⁻¹.

Fig. 2. From left to right: Different contributions to the two-point correlation function from z = 0 to z = 2. The blue circles represent the measurement for box 50/A for CDM cosmology considering halos with M > 108 M⊙ h−1. The corresponding errors are obtained with the jackknife method. The one-halo term is given by the dotted purple line and it is the sum of ξcs (solid orange line) and ξss (solid yellow line). The linear two-halo term (cyan dashed line) is then corrected for integral constraint, nonlinearities and halo exclusion (dashed red line). Finally, the solid black line represents the total best-fit model.

4. Results

In this section, we present the main results of our study on the halo occupation function, the subhalo distribution, and clustering. We discuss the outcome of the measurements in the CDM and alternative dark matter scenarios, and we then focus on our constraints on the HOD and density profile parameters. The Bayesian analysis was performed by adopting a Gaussian likelihood for each redshift bin

$$\begin{array}{c}\hfill ln\mathcal{L}\propto -\frac{1}{2}{\displaystyle \sum _i^{N_{\mathrm{bin}}}}{\left(\frac{m_i-d_i}{{\sigma }_i}\right)}^2,\end{array}$$(17)

where m_i is the model on the radial scale r_i, d_i is the measurement, and σ_i is the associated error.

4.1. Number counts

For each FOF group in the box at redshift z, the AIDA-TNG project gives access to the mass M_200c. Around each halo position, (x_h, y_h, z_h), we build a sphere of comoving radius R_200c^com and count the number of subhalos within the region. Then, we split the group catalog by averaging the counts into ten logarithmic mass bins. In the same way, we also measured ⟨N_h(N_h − 1)⟩.

Fig. 3 clearly separates the contribution of the central halo from that given by the associated substructures, namely the satellite halos. The applied mass selection ensures that there is at least one halo for M_200c > M_min, while subhalos are counted from M_200c > M₁. In the top panels we show the average number of pairs, ⟨N_h(N_h − 1)⟩. As explained in Sect. 3.1, this quantity enables us to determine whether the subhalo population follows a Poissonian distribution. Although the differences between cosmologies are negligible, we observe that the hypothesis is fully satisfied when M_200c > M₁, where the number of satellites begins to deviate from zero.

Fig. 3. From left to right: Average number of dark matter pairs (top panels), halo occupation function (middle panels) and residuals (bottom panels) with respect to the CDM reference level (solid blue line) as a function of M200c in box 50/A in the range 0 < z < 2 for CDM (blue circles), WDM3 (black triangles), WDM5 (crimson diamonds), SIDM1 (gold pentagons), and vSIDM (green triangles). The corresponding errors are computed with the standard deviation of the measurements for each mass bin. The dashed horizontal gray line represents the expected ratio for a Poissonian distribution. The solid horizontal gray line marks the central halo that always exists for M200c > Mmin.

In Fig. 3 we also show the occupation function ⟨N_h⟩ for each dark matter model. Alternatives to CDM are characterized by a reduction in the average number of satellites, with larger departures for the more massive systems. In WDM scenarios, this is particularly pronounced. For WDM3, the comparison with CDM shows residuals that exceed -0.5 for M_200c > 10¹¹ M_⊙ h⁻¹. A different mass selection, namely M_min > 10⁹ M_⊙ h⁻¹, significantly attenuates the effect, clarifying its strong dependence on the mass as a consequence of free streaming. For such a mass cut, only the warmest model, namely WDM1, exhibits a difference with respect to WDM3 and CDM, as we show in Fig. 4. Here, we fit Eq. (5) to quantify the dependence of α and M₁ on the redshift and mass selection of our sample; M₀, which represents only an offset in the mass of the group in the numerator, is not constrained and its presence does not significantly modify the mean values of α and M₁, although it widens their posterior distributions. For the sake of example, in Fig. 4 we show some results for boxes 50/A and 50/B. In particular, we note that for M_min > 10⁸ M_⊙ h⁻¹, we have $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=10.{42}_{-0.18}^{+0.25}$ at z = 0 and $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=10.{24}_{-0.16}^{+0.20}$ at z = 2, reflecting the fact that despite the presence of more massive structures, the average number of satellites in a given mass bin is slightly reduced in the local Universe by tidal disruption and merging phenomena. In fact, due to the hierarchical growth of cosmic structures, we first expect the formation of low-mass halos, which represent the progenitors of the systems that we count at z = 0, which can host more than 10³ subhalos. In contrast, $\alpha =1.{06}_{-0.08}^{+0.10}$ at z = 0 and $\alpha =1.{05}_{-0.07}^{+0.10}$ at z = 2, which means that ⟨N_h⟩ increases approximately with the same slope at different redshifts.

Fig. 4. Posterior distributions of the halo occupation function, ⟨Nh⟩, for different dark matter cosmologies (left panel), mass selections (middle panel), and redshifts (right panel). Different dark matter models are tested in box 50/B, with a mass selection of Mmin > 109 M⊙ h−1. The other two panels consider the CDM cosmology in box 50/A.

As also shown in Fig. 4, a higher minimum mass tends to significantly increase the value of the normalization, while the slope rises only slightly. At z = 0, for M_min > 10⁹ M_⊙ h⁻¹ we have $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=11.{47}_{-0.16}^{+0.23}$ and $\alpha =1.{17}_{-0.09}^{+0.14}$, while for M_min > 10¹⁰ M_⊙ h⁻¹ we have $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=12.{25}_{-0.17}^{+0.20}$ and $\alpha =1.{16}_{-0.12}^{+0.17}$, indicating that more massive halos tend to have relatively more massive substructures (Zehavi et al. 2005). Finally, the broadening of the posterior distribution is a consequence of the lower statistics.

Since in SIDM models the number density of dark matter halos does not differ significantly from the CDM case, as shown in Despali et al. (2025), the reduction in the average number of subhalos can be interpreted as a consequence of the more diffuse central density peak in self-interacting models, which results in satellite density profiles which are in general shallower around the core. In other words, moving toward outer regions, for a given mass, the cumulative radial count of subhalos presents a slope higher than the corresponding CDM model. However, these differences only result in a minimal increase and decrease in M₁ and α, respectively, which, however are not sufficient to discriminate between dark matter models.

4.2. Radial distribution of subhalos

In order to disentangle the contribution of the central and satellites to the clustering signal of dark matter halos, we need to study how satellites are distributed within their hosts. To do this, we split the halos into ten logarithmic-spaced mass bins, with a lower limit given by the transition mass from the central-dominated to the satellite-dominated regime of Fig. 3, at z = 0, which depends on the resolution of the boxes, and corresponds to 10¹⁰ M_⊙ h⁻¹ for 50/A. For each mass bin, we consider N_bin = 8 logarithmic-spaced radial bins. In practice, from the middle value of each mass bin we estimate the corresponding comoving radius R_200c^com, which is adopted as R_max, while R_min = R_max/(N_bin − 1), to exclude the innermost regions, where there is a lack of satellites due to mass selection, tidal disruption, and resolution effects. The counts are then averaged for every mass bin, and normalized according to the comoving volume of the i−th radial shell, $V_i=\frac{4}{3}\pi (R_{\max ,i}^3-R_{\min ,i}^3)$, obtaining ${\stackrel{\sim }{n}}_s(r)$, our stacked profile.

From measurements of the radial subhalo distribution, we can quantify possible deviations from a pure NFW model. Dealing with discrete tracers rather than a continuous dark matter density field, we compute the characteristic number density of satellites as

$$\begin{array}{c}\hfill n_s(r)=\frac{{\stackrel{\sim }{N}}_s}{M_{200c}}\rho (r),\end{array}$$(18)

where M_200c corresponds to the average in a given mass bin, ρ(r) is given by Eq. (9), and

$$\begin{array}{c}\hfill {\stackrel{\sim }{N}}_s=4\pi {\displaystyle {\int }_{R_{\min }}^{R_{\max }}}{r}^2{\stackrel{\sim }{n}}_s(r)\mathrm{d}r\end{array}$$(19)

is the normalization factor measured from the radial density profile, ${\stackrel{\sim }{n}}_s(r)$.

Figure 5 shows the radial distribution of subhalos normalized to the average value of the satellite density at a radius of R_200c, for standard ΛCDM and the alternative dark matter models. Variations among different masses are highlighted, while for r > 0.5 × R_200c the rescaled radial profiles are relatively self-similar; a feature that, however, is lost at higher redshifts because of the selection effects, which impact more when the average mass of structures is lower. In general, more massive structures are denser in terms of the number of satellites. This is because the effect of dynamical friction is progressively more relevant for massive subhalos hosted by large dark matter groups. As a consequence, they sink toward the halo center, reducing the radius of their orbits and enhancing the inner part of the density profile, while the outer region is less affected (Han et al. 2016; Chang et al. 2018). This fact incidentally explains why luminous galaxies are poor tracers of the NFW dark matter density distribution, while they match the corresponding subhalo profiles found in cosmological simulations (Watson et al. 2012; Chang et al. 2018). Conversely, since the halo mass is inversely correlated with the age, smaller hosts at z = 0 are, on average, older and have had more time to destroy subhalos located near the center. Interestingly, the situation is reversed in SIDM1 model, for which low-mass groups, with log[M_200c/(M_⊙ h⁻¹)] of 10.7 and 11.2 are centrally more peaked than the massive ones, with log[M_200c/(M_⊙ h⁻¹)] of 12.6 and 13.1, which have thermalized earlier. In fact, self-interaction between dark matter particles produces more diffuse central peaks and more extended cores for high-mass halos (Despali et al. 2025). While in vSIDM scenario the cross-section is a decreasing function of the relative velocity of particles, in SIDM1 it is constant. Therefore, moving toward the typical galaxy cluster size, this scenario is strongly constrained (Meneghetti et al. 2001; Eckert et al. 2022; Shen et al. 2022) since it could cause major departures from the ΛCDM cosmology; these are, however, ruled out by observations that found cusps at the center of galaxy clusters (Newman et al. 2013).

Fig. 5. From left to right: Radial distribution of dark matter subhalos at z = 0 for various mass bins and dark matter models. The triangles represent bins with a group mass of log[M200c/(M⊙ h−1)] = 10.7 (green), 11.2 (gray), 12.6 (gold), and 13.1 (red), respectively. The solid, colored lines reflect the gNFW best-fit model. The corresponding errors are computed with the standard deviation of the measurements, for each mass bin.

In Fig. 5 we show the best fit, calculated according to Eq. (18). In addition to gNFW, we also test the density profile model presented in Einasto (1965). Although it performs well if we analyze the full density profile up to scales larger than R_200c, it fails at small separations and for low mass halos, for which the friction force is not enough to flatten the spatial distribution of satellites within the core. We estimate the parameters of the gNFW profile, γ and f_c, choosing large, flat priors, between [ − 10, 2] and [0.1, 1.5], respectively. In Fig. 6, for a ΛCDM cosmology, we show the mass and redshift evolution of the density profile parameters, considering only mass bins where halos are present at all redshifts and for all cosmologies. In particular, γ assumes negative values, reflecting the well-established results that the halo distribution is antibiased (Ghigna et al. 1998; Springel et al. 2008), namely it is more cored toward the center than the matter density. The analytical model presented in Han et al. (2016) explains the spatial distribution of subhalos as a result of tidal stripping, generated from dynamical friction which, conversely, does not affect the dynamics of dark matter particles. Furthermore, γ increases with the mass of the group M_200c, which means that less massive systems generally present a shallower radial distribution of satellites. However, we should consider that a variation in M_200c from approximately 10¹³ M_⊙ h⁻¹ to 10¹⁰ M_⊙ h⁻¹ translates into a reduction of R_200c of about one order of magnitude; therefore we explore a shorter range of spatial separations in low-mass groups. In this case, the value of γ reflects the lack of subhalos toward more central regions, due to stripping phenomena, box resolution and selection effects. Also for this reason, the majority of subhalos are located in the outer regions of the group, which correspond to shells of larger volume.

Fig. 6. Marginalized posterior distribution of the gNFW parameters, γ (left) and fc (right), for different mass bins. Top panels: Evolution with redshift for a fixed ΛCDM cosmology. Bottom panels: Evolution in alternative dark matter models at z = 0. The errors represent the 68% percentiles.

In high-z Universe, dark matter subhalos form in an environment with a higher mean matter density, which allows a more efficient accretion process. In addiction, they are less massive and more concentrated. These features translate into a reduction of the tidal stripping. Furthermore, over time, their average number decreases as they merge. Having measured the density in comoving radial shells, only the net effect of satellite creation and destruction is considered, while the total mass of the central peak increases toward z = 0. Thus, in time the satellite density profiles become flatter. This explains the redshift evolution of the γ − M_200c relation that we see in Fig. 6. On the other hand, the concentration factor, f_c, is roughly constant for halos in a wide range of masses and redshift, which share the same deviation from the adopted c − M_200c relation. This simply describes a general shrinking in the characteristic radius r_s with respect to that of the dark matter density profile, reflecting the loss of orbital energy of the subhalos that experience dynamical friction (Chang et al. 2018).

In Fig. 6 we also see the variation of the gNFW density profile parameters in the different dark matter scenarios. Although the general mass and redshift tendencies are kept, there is no clear dependence of f_c on cosmology. Instead, WDM models prefer a higher γ, which means that the subhalos radial distribution is more concentrated and steeper around the center, as also found in the analysis of the Aquarius-WDM runs (Lovell et al. 2014), the COCO simulations (Bose et al. 2017) and zoom-in simulations of early type galaxies from the EAGLE box (Despali et al. 2020). Indeed, the free-streaming of dark matter particles determines an overall delay in the structure formation process, especially for less massive systems (Lovell 2024). Therefore, WDM subhalos of a given mass generally form later than their CDM counterparts, when the average density of the Universe is smaller, allowing more concentrated structures to grow. Conversely, in SIDM scenarios, γ is generally lower, which means that the innermost satellite density profile is shallower and less concentrated. This result has been already found in simulations by Fischer et al. (2022, 2024), as a reduction in the number of subhalos becoming progressively stronger at smaller distances to the host, and can be interpreted as a consequence of core thermalization due to particle interaction. This is in fact expected to produce a flatter underlying dark matter distribution and make satellites more sensitive to tidal disruption. For this reason, the relative number of subhalos is distributed toward the external radial shells of the FOF group. At z = 0, for M_200c > 10^12.5 M_⊙ h⁻¹, the radial distribution of the SIDM1 subhalos becomes shallower than that of the vSIDM. This once again highlights the tendency of SIDM1 to form cores even in the most massive systems, while vSIDM interactions are only significant for low-mass bound structures. The behavior of WDM and SIDM cosmologies, quantitatively appreciable already in Fig. 5, is also reflected in their clustering properties, as explained in Sect. 4.3.

4.3. The two-point correlation function

The clustering information of a sample of discrete tracers is encoded in the two-point correlation function, ξ(r). According to the definition, for a pair of elements located in volume elements δV₁ and δV₂, at a given separation r, it is expressed as an excess probability with respect to a uniform distribution,

$$\begin{array}{c}\hfill \delta P_{12}(r)=n_V^2[1+\xi (r)]\delta V_1\delta V_2,\end{array}$$(20)

where n_V is the mean tracer number density in a unit of comoving volume.

We build the random catalog by uniformly extracting the comoving coordinates of dark matter halos in a box. Since for each redshift the number of tracers is higher than 4 × 10⁵, the shot noise is limited, and the size of the random catalog is kept equal to the size of the original one, in order to speed up the computation. We measure the real-space two-point correlation function with the Landy & Szalay (1993) estimator,

$$\begin{array}{c}\hfill {\xi }_{\mathrm{obs}}(r)=\frac{DD(r)+RR(r)-2DR(r)}{RR(r)},\end{array}$$(21)

which provides an unbiased measurement of ξ_obs(r) once the normalized number of data-data, random-random, and data-random pairs, DD(r), RR(r) and DR(r) respectively, is counted. The errors are estimated using the jackknife method, considering a number of regions equal to 125 (Norberg et al. 2009). We adopt appropriate mass cuts to ensure completeness, corresponding to a lower limit of 10⁸ M_⊙ h⁻¹ for 50/A. The two-point correlation is measured into fifteen logarithmic-spaced bins between 0.05 Mpc h⁻¹ and 15 Mpc h⁻¹. The lower limit is motivated by the resolution of the simulation, whereas the upper limit is due to the finite size of the box and corresponds approximately to the scale at which the correlation function becomes negative. At these scales nonlinear effects are not negligible and the halo bias depends on r, as well as on redshift and mass. A higher-mass cut selects only more massive halos, which are more biased, and therefore provides a stronger clustering signal. Furthermore, with cosmic time, the one-halo term becomes more pronounced, implying a progressive increase of nonlinearities.

The use of cosmological simulations with size much smaller than the volume of the Universe has non-negligible consequences on the measurement of the clustering signal. From a theoretical point of view, a scale-free power spectrum satisfies the condition P(k = 0) = 0, which corresponds to imposing homogeneity of matter fluctuations on sufficiently large scales, where the cosmological principle is satisfied (Riquelme et al. 2023). According to the definition of probability, the integral of Eq. (20) over the full volume of the Universe is equal to 1. The same condition also holds if we consider the finite volume of the box, V_b. This implies that

$$\begin{array}{c}\hfill {\displaystyle \int {\xi }_{\mathrm{obs}}}(r)\mathrm{d}{V}_b=0,\end{array}$$(22)

which means that we should find a scale at which ξ_obs(r) becomes negative, which is lower than the corresponding scale for the model of ξ(r). In other words, the evaluation of the correlation function in a limited volume determines a reduction in the clustering signal. The observed correlation function is linked to the theoretical one through a window function, W(r). In Eq. (22), in order to avoid volume integration in surveys with complicated geometries, it is common to correct the model prediction by following the radial evolution of random-random pairs. This information is contained in the shape of the random catalog. Such ‘integral constraint’ condition therefore takes the form (Ross et al. 2013; Riquelme et al. 2023)

$$\begin{array}{c}\hfill {\xi }_{\mathit{IC}}(r)=\xi (r)-\frac{{\sum }_i\xi (r_i)RR(r_i)}{{\sum }_{i}RR(r_i)},\end{array}$$(23)

and spans over the range at which the pairs are counted.

In Fig. 7, we compare the measures obtained in box 50/A within the ΛCDM framework, with the results in alternative dark matter scenarios, finding that the main difference between models is located at the one-halo scales, corresponding to r ≲ 1 Mpc h⁻¹. In fact, SIDM1 and vSIDM are systematically less clustered, and the difference to CDM is only slightly affected by time evolution, since they residuals remain between −5 and −2.5 for all the redshifts considered. In particular, it is noteworthy that SIDM1 is closer to CDM at high z, indicating a late-time thermalization of the core; conversely, the vSIDM model deviates more at z = 2. In other words, vSIDM is more effective at high z, since the cross-section is higher for low-mass halos, where the velocity dispersion of the particles is smaller. On the other hand, in WDM cosmologies, the average delay in the structure formation process makes clustering differences more visible at high z, with residuals larger than 2.5 for r ≲ 1 Mpc h⁻¹ and r ≲ 0.2 Mpc h⁻¹, for WDM3 and WDM5, respectively. These results reinforce what we already found from the analysis of Sect. 4.2, according to which WDM generally produces a more peaked subhalo density distribution. As we have already seen in Sect. 3.3 the satellite-satellite pairs dominate the one-halo range. Consequently, the radial distribution of the satellites directly affects the clustering properties of the full halo population. However, the excess of clustering is also found in the two-halo term for the warmest models. Again, we can interpret this trend as a consequence of free-streaming of the particles. Indeed, the suppression of low-mass dark matter halos has a direct impact on the halo mass function, which operates as a weighting factor for the halo bias, resulting in practice in a larger effective bias of our sample.

Fig. 7. From left to right: Residuals of the two-point correlation function for alternative dark matter models for redshifts from z = 0 to z = 2. The blue, horizontal line represents the reference level for ΛCDM model. WDM3 (black triangles), WDM5 (crimson diamonds), SIDM1 (gold pentagons) and vSIDM (green triangles) are also shown.

A proper modeling of the observed two-point correlation function is well beyond the linear theory, requiring a generalization of the density profiles and the introduction of the HOD formalism. Our ultimate goal is to explore the degeneracy of HOD parameters in different cosmological frameworks. For this reason, as already described in Sect. 3.3, we calibrate the nonlinear deviation from b_h²P_lin(k, z), by modeling the clustering of dark matter groups only, without considering the contribution of satellites. In this way, we quantify the A + Bk + Ck² factor of Eq. (16), which incorporates the exclusion effects and the inaccuracies of the nonlinear bias model for each redshift and dark matter scenario, adopting their median value for the HOD characterization of the full halo populations.

Then we fit Eq. (10) to constrain the parameters M₁ and α, whose posterior distributions are not affected by the presence of M₀. We set the latter to zero, to reduce the number of parameters and thus the complexity of the model. In order to assume a physically motivated normalized density profile, ${\stackrel{\sim }{u}}_s(k,M,z)$, we adopt flat priors on the parameters f_c and γ, according to the range permitted by the posterior distributions presented in Sect. 4.2, while log[M₁/(M_⊙ h⁻¹)] and α are left free to vary in the range [9, 13] and [0.7, 1.6], respectively. For the WDM models, we use Eq. (1) to suppress the linear power spectrum, P_lin(k, z), at small scales, while for the halo mass function, n(M), we apply Eq. (3). This reduces the low-mass halo population, resulting in a higher effective bias, which determines an enhancement of the two-halo term. SIDM models, on the other hand, share the same power spectrum and mass function as ΛCDM.

In Fig. 8 we show the marginalized posterior distribution in the M₁ − α plane. The constraints on HOD parameters are in agreement with the results obtained in Sect. 4.1 with halo occupation function. In particular, for the ΛCDM cosmology, we find $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=10.{41}_{-0.12}^{+0.15}$ and $\alpha =1.{06}_{-0.04}^{+0.05}$ at z = 0; $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=9.{93}_{-0.09}^{+0.09}$ and $\alpha =0.{97}_{-0.03}^{+0.04}$ at z = 2. However, they are also more competitive. In fact, the two-halo term of the correlation function relies on ⟨N_h⟩, as in the case of number counts, while the one-halo term depends on the number of pairs ⟨N_h(N_h − 1)⟩. In addition to second-order statistics, the clustering model also includes stringent prior information from the parameters of the gNFW density profile, f_c and γ. This emphasizes the necessary synergy between these different probes for future HOD-based studies.

Fig. 8. Marginalized posterior distributions of the HOD parameters in the M1 − α plane for box 50/A and M > 108 M⊙ h−1. The 68% and 95% contours are represented for CDM (blue), WDM3 (black), WDM5 (crimson), SIDM1 (gold) and vSIDM (green) cosmological models.

Clustering of halos reveals its full potential when applied to alternative dark matter scenarios. The M₁ parameter indeed increases in the WDM and SIDM models, alongside a slight decrease in α. This reflects the fact that the average reduction in the number of satellites due to free streaming or core thermalization means that higher group masses are required for finding the same number of subhalos. This tendency was already observed in Sect. 4.1, where we analyzed the occupation functions; however, the ability to distinguish between different models was limited to box 50/B, for which only WDM1 presented a 1σ difference with CDM, for the parameter M₁. From this point of view, the study of the two-point correlation function represents a powerful probe, especially if we analyze the clustering of dark matter halos at high z. In particular, for the thermal WDM3 model we find $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=10.{83}_{-0.11}^{+0.11}$ and $\alpha =1.{04}_{-0.04}^{+0.04}$ at z = 0; $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=10.{16}_{-0.06}^{+0.07}$ and $\alpha =0.{89}_{-0.03}^{+0.03}$ at z = 2. In other words, the normalization in the number of satellites exhibits a difference greater than 2σ with respect to CDM, while at z = 2, the one on the slope is at the level of 1σ. WDM5 occupies an intermediate position, with $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=10.{53}_{-0.11}^{+0.13}$ and $\alpha =1.{06}_{-0.04}^{+0.04}$ at z = 0; $log[M_1/({\mathrm{M}}_{\odot }{h}^{-1})]=9.{95}_{-0.08}^{+0.08}$ and $\alpha =0.{93}_{-0.03}^{+0.03}$ at z = 2, and therefore requires more stringent mass cuts to achieve the same level of disagreement. Inspecting Fig. 8, it is also clear that fixed and velocity-dependent SIDM models are indistinguishable from each other. As the latter acts mainly where the velocity dispersion of the particles is lower, altering the structure and distribution of small-mass halos, we conclude that simulations with higher mass resolution are required to explore finer spatial scales in detail and to better highlight the differences between alternative dark matter models. In the future, we plan to exploit the full-physics AIDA-TNG simulations, for a subsequent work, which will take into account the fact that not all dark matter halos and subhalos are able to host galaxies, as well as the complicated relationship between the mass and the luminosity of galaxies, directly linked to their formation processes. In this way, with more complete HOD models, we aim at disentangling the effects of alternative dark matter from baryonic physics, in preparation for a full comparison with observational data from galaxy clustering.

5. Conclusions

We studied the effect of warm and self-interacting dark matter models on the abundance, the radial distribution, and the clustering properties of subhalos. For this purpose, we used a set of dark matter-only cosmological boxes with a variable size and resolution that were developed within the AIDA-TNG project. We applied the HOD formalism, which allowed us to separate the contributions of distinct populations of central and satellite halos, in order to model their occupation function and clustering. In particular, we counted the average abundances of subhalos within the comoving radius R_200c^com for different mass bins and found that the HOD parameter M₁ is higher at low redshift as a consequence of the impact of tidal stripping and merging phenomena on the average number of satellites. Alternative dark matter cosmologies predict a reduction in the number of subhalos within R_200c^com. In WDM models, particle free streaming indeed suppresses the formation of low-mass halos, while in the SIDM case, the core thermalization redistributes satellites in the outer region of dark matter groups. However, while for box 50/B for M > 10⁹ M_⊙ h⁻¹, CDM and WDM1 cosmologies exhibit an appreciable difference, in box 50/A for M > 10⁸ M_⊙ h⁻¹, the occupation function does not provide a significant way to distinguish between models, requiring higher statistics and/or mass resolution. Moreover, with the applied mass selections, we verified that for the group mass M_200c > M₁, where the average number of satellites is larger than one, they follow a Poissonian distribution for all the models considered.

We also evaluated the radial distribution of subhalos around the centers of the groups. To do this, we counted the average number density of satellites in spherical shells up to a radius equal to R_200c^com. For each mass bin and dark matter model, we fit the stacked radial density with a gNFW density profile. In this way, we were able to constrain the evolution of the gNFW parameters γ and f_c with M_200c and z. Although the latter is approximately constant and larger than one, as a consequence of dynamical friction, we found that the subhalo distribution becomes steeper for more massive systems at high-z, while it remains antibiased compared to the matter distribution. Furthermore, SIDM models generally prefer a shallower satellite radial profile as a result of the core thermalization in bounded dark matter halos, while WDM cosmologies on average admit a more cusped distribution toward the center of the groups as a result of the delayed structure formation process. This has direct consequences on the clustering properties of dark matter halos because it mainly affects the one-halo term of the two-point correlation function, where the differences between cold and alternative dark matter models are more relevant. In particular, in WDM cosmologies, the peaked subhalo distribution determines a clustering excess at r ≲ 0.2 Mpc h⁻¹, with an increasing effect toward high z. At z = 2, it indeed also extends to the two-halo scales; the free streaming of dark matter particles determines a low-mass suppression in the halo mass function, which increases the effective bias of our sample. Conversely, SIDM models are less clustered at small scales, and their residuals with respect to CDM are only slightly affected by the redshift evolution. The HOD modeling of the two-point correlation function provided constraints on M₁ and α that agree with the prediction of the halo occupation functions. In addition, it produced more stringent results, with a 2σ difference on M₁ between CDM and WDM3 from z = 0 to z = 2. This revealed the effectiveness of this cosmological probe in distinguishing between possible alternative dark matter scenarios.

In the future, we plan to conduct a study on full-physics simulations to explore the degeneracy between the properties of dark matter particles and the physics of baryons in macroscopic observables, such as abundances, density profiles, and clustering. Variations in the physics of dark matter modify the abundance and the spatial distribution of halos, with a possible direct and macroscopic impact on luminous tracers, such as galaxies. Therefore, the completeness, spatial resolution, and depth of forthcoming photometric and spectroscopic surveys, combined with new state-of-the-art simulations, will offer unique possibilities to shed light on the real nature of dark matter.

Acknowledgments

We thank Ravi Sheth for useful discussions and hints on the paper. We would like to thank the anonymous referee for the constructive comments, which have improved the quality of the manuscript. We acknowledge the financial contribution from the grant PRIN-MUR 2022 20227RNLY3 “The concordance cosmological model: stress-tests with galaxy clusters” supported by Next Generation EU and from the grants ASI n.2018-23-HH.0 and n. 2024-10-HH.0 “Attivitá scientifiche per la missione Euclid – fase E”, and the use of computational resources from the parallel computing cluster of the Open Physics Hub (https://site.unibo.it/openphysicshub/en) at the Physics and Astronomy Department in Bologna. GD acknowledges the funding by the European Union - NextGenerationEU, in the framework of the HPC project – “National Centre for HPC, Big Data and Quantum Computing” (PNRR – M4C2 – I1.4 – CN00000013 – CUP J33C22001170001). We acknowledge ISCRA and ICSC for awarding this project access to the LEONARDO supercomputer. The AIDA-TNG simulations where run on the LUMI supercomputer (Finland), thanks to the award of computing time through a EuroHPC Extreme Scale Access call.

References

All Tables

All Figures

	Fig. 1. From top to bottom: Spatial distribution of dark matter halos in box 50/B for the CDM, WDM3, and WDM1 cosmologies. The color bar shows the corresponding redshift in the range 0 < z < 2.
In the text

Fig. 2. From left to right: Different contributions to the two-point correlation function from z = 0 to z = 2. The blue circles represent the measurement for box 50/A for CDM cosmology considering halos with M > 108 M⊙ h−1. The corresponding errors are obtained with the jackknife method. The one-halo term is given by the dotted purple line and it is the sum of ξcs (solid orange line) and ξss (solid yellow line). The linear two-halo term (cyan dashed line) is then corrected for integral constraint, nonlinearities and halo exclusion (dashed red line). Finally, the solid black line represents the total best-fit model.

In the text

Fig. 3. From left to right: Average number of dark matter pairs (top panels), halo occupation function (middle panels) and residuals (bottom panels) with respect to the CDM reference level (solid blue line) as a function of M200c in box 50/A in the range 0 < z < 2 for CDM (blue circles), WDM3 (black triangles), WDM5 (crimson diamonds), SIDM1 (gold pentagons), and vSIDM (green triangles). The corresponding errors are computed with the standard deviation of the measurements for each mass bin. The dashed horizontal gray line represents the expected ratio for a Poissonian distribution. The solid horizontal gray line marks the central halo that always exists for M200c > Mmin.

In the text

	Fig. 4. Posterior distributions of the halo occupation function, ⟨Nh⟩, for different dark matter cosmologies (left panel), mass selections (middle panel), and redshifts (right panel). Different dark matter models are tested in box 50/B, with a mass selection of Mmin > 109 M⊙ h−1. The other two panels consider the CDM cosmology in box 50/A.
In the text

Fig. 5. From left to right: Radial distribution of dark matter subhalos at z = 0 for various mass bins and dark matter models. The triangles represent bins with a group mass of log[M200c/(M⊙ h−1)] = 10.7 (green), 11.2 (gray), 12.6 (gold), and 13.1 (red), respectively. The solid, colored lines reflect the gNFW best-fit model. The corresponding errors are computed with the standard deviation of the measurements, for each mass bin.

In the text

	Fig. 6. Marginalized posterior distribution of the gNFW parameters, γ (left) and fc (right), for different mass bins. Top panels: Evolution with redshift for a fixed ΛCDM cosmology. Bottom panels: Evolution in alternative dark matter models at z = 0. The errors represent the 68% percentiles.
In the text

	Fig. 7. From left to right: Residuals of the two-point correlation function for alternative dark matter models for redshifts from z = 0 to z = 2. The blue, horizontal line represents the reference level for ΛCDM model. WDM3 (black triangles), WDM5 (crimson diamonds), SIDM1 (gold pentagons) and vSIDM (green triangles) are also shown.
In the text

	Fig. 8. Marginalized posterior distributions of the HOD parameters in the M1 − α plane for box 50/A and M > 108 M⊙ h−1. The 68% and 95% contours are represented for CDM (blue), WDM3 (black), WDM5 (crimson), SIDM1 (gold) and vSIDM (green) cosmological models.
In the text