Einasto profile as the halo model solution coupled to the depletion radius

Abstract

We constrain the halo profiles outside the halo boundaries by solving for the matching profiles required by the halo model. In the halo model framework, the matter distribution in the universe can be decomposed into the spatial distribution of halos convolved with their internal structures. This leads to a set of linear equations in Fourier space which uniquely determines the optimal halo profiles for any given halo catalog. In this work, we construct three halo catalogs with different boundary definitions, and solve for the optimal profiles in each case using measurements of halo-matter and halo-halo power spectra. Our results show that for a given halo field, there is always a set of matching profiles to accurately reconstruct the input statistics of the matter field, even though it might be complex to model the profiles analytically. Comparing the solutions from different halo catalogs, we find their mass distributions inside the inner depletion radii are nearly identical, while they deviate from each other on larger scales, with a larger boundary resulting in a more extended profile. For the depletion radius based catalog, the numerical solution agrees well with the Einasto profile. Coupling the Einasto profile with the depletion catalog, the resulting halo model can simultaneously predict the halo-matter power spectra to $10\%$ and matter-matter power spectrum to $5\%$, improving over conventional models in both the interpretability and versatility.

1 Introduction

By assuming that all mass is contained within virialized objects called dark matter halos, the mass distribution of the Universe can be decomposed into the spatial distribution of halos, convolved by their internal structures. On large scales, the internal structures of halos become unimportant, and the halo model has been quite successful in tracing the large-scale structure (see e.g., Cooray & Sheth, 2002; Asgari et al., 2023, for reviews) of the Universe. On intermediate and small scales, however, an accurate halo model requires detailed modelling of the structures and boundaries of halos.

Conventionally, a halo is defined as a virialized object according to the virial radius, $R_{\rm vir}$, as motivated by the spherical collapse model (Gunn & Gott, 1972). Using cosmological simulations, the spherically averaged halo profile out to the virial scale has been found to be well described by some empirical fitting functions such as the Navarro-Frenk-White (NFW, Navarro et al., 1995, 1996, 1997) or Einasto (Einasto, 1965; Merritt et al., 2006; Navarro et al., 2004, 2010) profiles. However, the halo model constructed using such profiles truncated at $R_{\rm vir}$ does not fully match simulations, especially on the transition scale where halos start to intersect (e.g., García & Rozo, 2019; Mead et al., 2021; Zhou & Han, 2023). Some studies have made efforts to address this issue by introducing some global parameters or additional corrections on transition scales (Tinker et al., 2005; van den Bosch et al., 2013; Mead et al., 2015; Philcox et al., 2020), which further complicates the halo model. Therefore, to build a more concise and accurate halo model, more efforts are needed to quantify the matter distribution around the boundary of a halo as well as the definition of the halo boundary itself.

Physically, the virial radius is expected to only describe an equilibrium structure, while a growing halo is inevitably surrounded by a non-equilibrium region which extends beyond the virial radius. Significant gravitational influences of a halo on its satellites also start from a radius much larger than the virial radius (Ludlow et al., 2009; Bahé et al., 2013; Behroozi et al., 2014). Thus, the definition of the halo boundary should be revised to allow for a more extended profile beyond the virial radius. Many recent works have attempted to introduce new boundaries to more physically define a halo. A widely studied candidate is the splashback radius, which is defined at the first apocenter of the orbit of the an infalling particle in a growing halo (Fillmore & Goldreich, 1984; Bertschinger, 1985; Adhikari et al., 2014; Diemer & Kravtsov, 2014; Shi, 2016; Mansfield et al., 2017). In practice, it is often estimated from the steepest location in the halo density profile (More et al., 2015). The macroscopic effect of halo growth is to cause a drop in its environmental density. Accordingly, Fong & Han (2021) proposed a new halo boundary called the depletion radius which separates the growing part of a halo from the fading environment. According to continuity, the depletion radius can be found at the radius of the maximum mass infall rate.The depletion of the environment due to halo growth is shown to be responsible for the creation of a minimum in the bias profile around a halo, leading to an alternative representation of the depletion radius in the bias domain (Fong & Han, 2021; Gao et al., 2023). Some other works have attempted to define the radius of a halo using characteristics in the velocity profile (e.g., Cuesta et al., 2008; Bose & Loeb, 2021; Pizzardo et al., 2023a), phase space structure (Tomooka et al., 2020; Aung et al., 2021) or the mass profile (Pizzardo et al., 2023b). These boundary definitions characterize the halo with different physical motivations, providing us complementary insights to the structure of a dark matter halo.

Ideally, the halo profile should be self-consistent with its boundary definition so that the profile terminates at the boundary. However, because halos are intrinsically aspherical (Jing & Suto, 2002; Allgood et al., 2006; Mansfield et al., 2017; Wang et al., 2022), and more importantly, because of ongoing mergers which temporarily extend and distort the domain of a halo, the spherically average density profile can not be crudely truncated at the halo boundary. Instead, a smooth truncation is desired to better describe the halo structure outside a spherical boundary in practice.

A smoothly truncated profile leads to a question: what is the profile outside a given halo boundary? Directly measuring the halo outer profile is difficult in simulations because it is unclear how mass should be partitioned among neighboring halos. To mitigate this problem, some recent works have considered halos and a background density field separately. For example, by partitioning particles around a halo in phase space, the halo profile can be split into orbiting and infalling parts (Diemer, 2022, 2023; García et al., 2023; Salazar et al., 2024). Chen & Afshordi (2023) considered the linear density field as a background and derived the profiles of halos in excess of this background field in Fourier space. Because the density field is a convolution of the halo field with the density profiles of halos, in Fourier space the halo profiles can be solved explicitly and self-consistently from the matter and halo fields. However, because they considered the halo profile in excess of the linear density field, the resulting profiles are not directly applicable to the classical halo model and can have negative values.

In this work, we will work in the classical halo model framework and study the outer profiles of halos in light of recent developments in the halo boundary definition. Using a similar method to Chen & Afshordi (2023), we will solve for the halo profiles in Fourier space and investigate how the recovered profiles depend on the adopted halo boundary of the input catalog. Transforming the solutions back to real space, we will investigate the analytical properties of the reconstructed profiles, as well as their performances in reproducing additional large scale structure statistics when inserted back to a halo model. Through these analysis, we aim to answer a series of fundamental questions regarding halo profile and halo boundary, including i) For each given halo catalog corresponding to a given halo boundary, is there always a matching set of halo profiles to be used in the halo model? ii) With the matching profile to the halo boundary, how well does the resulting halo model work in predicting additional large scale structure statistics not used to constrain the model? iii) For different boundary choices, is there an optimal boundary to produce the best halo model, and what is the physical implication of such a choice if it exists?

This paper is organized as follows. In Section 2 we describe the method used for constraining the halo profile. In Section 3, we describe the simulation data and introduce three halo catalogs corresponding to different boundary definitions. Section 4 shows the reconstructed profiles for the three catalogs, and evaluates their performances in predicting the matter-matter power spectrum. We discuss the stability of our method on large scales and the impacts of unresolved mass in Section 6. Finally, in Section 7, we summarize the results of this paper.

2 Method

In this section we derive the set of linear equations for constraining the halo profiles in Fourier space, and explain how we account for unresolved halos and diffuse matter which distinguishes our model from some other similar methods.

2.1 Constraining halo profiles in the halo model framework

The cross power spectrum of the matter field and a halo population with mass $M$ and mean number density $n(M)$ can be expressed as

$$P_{\rm hm}(k|M) =P^{\rm 1h}_{\rm hm}(k|M)+P^{\rm 2h}_{\rm hm}(k|M),$$ (1)

$$P^{\rm 1h}_{\rm hm}(k|M) =\frac{W(k|M)}{\bar{\rho}_{\rm m}},$$ (2)

$$P^{\rm 2h}_{\rm hm}(k|M) =\frac{1}{\bar{\rho}_{\rm m}}\int n(m)P^{\rm 2h}_{\rm hh}(k|m,M)W(k|m)\mathrm{d}m,$$ (3)

where $\bar{\rho}_{\rm m}$ is the mean density of the Universe, $P_{\rm hh}^{\rm 2h}$ is the two halo term of the halo-halo power spectrum between halos of mass $m$ and $M$, adn $W(k|m)$ is the halo density profile in Fourier space. The relations of $k$-space profile $W(k|M)$ and real-space profile $\rho(r|M)$ are given by

$$W(k) =\int^{\infty}_{0}\rho(r|M)\frac{\sin(kr)}{kr}4\pi r^{2}\mathrm{d}r,$$ (4)

$$\rho(k) =\frac{1}{(2\pi)^{3}}\int^{\infty}_{0}\rho(r|M)\frac{\sin(kr)}{kr}4\pi k^{2}\mathrm{d}k.$$ (5)

In this case, $W(k)$ converges to the integrated mass $M_{\rm int}$ of the density profile when $k$ goes to 0. Equations (1) to (3) show clearly that the halo-matter power spectrum is the linear combination of the halo profiles, which implies that the halo profiles can be completely solved with enough independent constraints.

We now proceed to divide halos into $l$ mass bins to establish such constraints. Considering the halo-matter power spectrum of the $i$th mass bin at $k$, Equation (1) can be rewritten as,

$$P_{\rm hm}(k|m_{i}) =\frac{W(k|m_{i})}{\bar{\rho}_{\rm m}}$$

$$+\frac{1}{\bar{\rho}_{\rm m}}\sum^{l}_{j=1}n(m_{j})P^{\rm 2h}_{\rm hh}(k|m_{j},m_{i})W(k|m_{j}).$$ (6)

Here the integral over mass has been replaced with the summation. Equation (6) contains $l$ unknown variables $W(k|m_{j})$ at each $k$. By combining the halo-matter power spectra from $l$ mass bins, a set of solution for $W(k|m_{j})(j=1..l)$ can be completely determined.

To organize the equations into a compact form, we define the following vectors

h $$=\left(\begin{array}[]{cccc}P_{\rm hm}(k|m_{1})\\ P_{\rm hm}(k|m_{2})\\ \vdots\\ P_{\rm hm}(k|m_{l})\\ \end{array}\right),$$ (11)

w $$=\left(\begin{array}[]{cccc}W(k|m_{1})\\ W(k|m_{2})\\ \vdots\\ W(k|m_{l})\\ \end{array}\right).$$ (16)

and tensors

H $$=\left(\begin{array}[]{cccc}P_{\rm hh}(k|m_{1},m_{1})&P_{\rm hh}(k|m_{1},m_{2})&\cdots&P_{\rm hh}(k|m_{1},m_{l})\\ P_{\rm hh}(k|m_{2},m_{1})&P_{\rm hh}(k|m_{2},m_{2})&\cdots&P_{\rm hh}(k|m_{2},m_{l})\\ \vdots\\ P_{\rm hh}(k|m_{l},m_{1})&P_{\rm hh}(k|m_{l},m_{2})&\cdots&P_{\rm hh}(k|m_{l},m_{l})\\ \end{array}\right),$$ (21)

N $$=\left(\begin{array}[]{cccc}n(m_{1})&\cdots&\cdots&0\\ 0&n(m_{2})&\cdots&0\\ \vdots&\vdots&\vdots&\vdots\\ 0&\cdots&\cdots&n(m_{l})\\ \end{array}\right)$$ (26)

So that, the constraints on the halo profiles can be rewritten as

$$\textbf{{h}}=\frac{1}{\bar{\rho}_{\rm m}}(\textbf{I}+\textbf{HN})\textbf{{w}}.$$ (27)

where I is the identity matrix. With h, H, N measured from the simulation data, w can be solved directly from the above equation, independently for each $k$ mode.

2.2 Accounting for unresolved halos and diffuse matter

Theoretically, the halo mass bins in Equation (27) should cover the complete mass spectrum of halos. However, due to free streaming of cold dark matter particles(e.g., Green et al., 2005; Profumo et al., 2006; Schneider et al., 2013), and due to the finite resolution of numerical simulations, halos can only be resolved and modeled down to a certain mass limit. It is thus necessary to introduce an unresolved mass component to account for the contribution from unresolved halos as well as from some potential diffuse matter, so that Equation (27) becomes

$$\textbf{{h}}=\frac{1}{\bar{\rho}_{\rm m}}(\textbf{I}+\textbf{HN})\textbf{{w}}+\textbf{{c}}.$$ (28)

Here we have defined a vector

$$\textbf{{c}}=\left(\begin{array}[]{cccc}P_{\rm hm}^{\rm u}(k|m_{1})\\ P_{\rm hm}^{\rm u}(k|m_{2})\\ \vdots\\ P_{\rm hm}^{\rm u}(k|m_{n})\\ \end{array}\right),$$ (29)

where $P_{\rm hm}^{\rm u}$ is the cross power spectrum between halo and unresolved mass.

Zhou & Han (2023) modeled the unresolved component $\xi^{\rm u}_{\rm hm}(r|m)$ in real space by approximating the unresolved halos as mass points distributed outside the exclusion radius, so that $\xi^{\rm u}_{\rm hm}(r|m)$ can be expressed as a universal halo-halo correlation $\hat{\xi}_{\rm hh}(r)$ multiplied by the effective bias of the unresolved component, $b_{\rm unr}$, with a truncation due to halo exclusion. The expressions of these quantities are detailed in Appendix A. We will model the unresolved halo-matter correlation function following this form and convert it to $P_{\rm hm}^{\rm u}$ in Fourier space.

Equation (28) forms the basis for our solution to the halo density profile. We also compare the solutions with and without the unresolved term in Section 6, and find that the unresolved term impact the mass distribution in the halo outskirts, especially for low mass halos.

3 Simulation and halo samples

3.1 Simulation data

We use a Lambda Cold Dark Matter simulation, which is one of the CosmicGrowth simulations (Jing, 2019) run with a P³M code, to extract the data of halo and matter fields. The simulation was run in a box of side length 600$h^{-1}\mathrm{Mpc}$ containing $3072^{3}$ particles with cosmological parameters $\Omega_{\rm m}=0.268$, $\Omega_{\Lambda}=0.732$, $h=0.704$, $n_{\rm s}=0.968$, and $\sigma_{8}=0.830$.

The candidate halos are identified by the the Friends-of-Friends (FoF) algorithm with a standard linking parameter of 0.2, and then processed with HBT+ (Han et al., 2012, 2018) to identify subhalos. We define the virial mass of a halo as the mass enclosed in a sphere with a virial density according to the spherical collapse model (Bryan & Norman, 1998), and the corresponding radius is the virial radius, $R_{\rm vir}$. The depletion radii of individual halos are estimated using the relation $R_{\rm id}$=2.1$R_{\rm vir}$ (Fong & Han, 2021; Gao et al., 2023).

We construct an original sample with about $2\times 10^{6}$ distinct halos by selecting candidate objects within a mass range of $11.50<{\rm log}_{10}[m/(h^{-1}{\rm M}_{\odot})]<15.35$ at z=0. These halos are divided into seven logarithmic mass bins with equal spacing. To investigate how the definitions of the exclusion radius affect the halo outer profiles, halos are further selected according to the exclusion criteria, resulting in three halo catalogs. We put more details about the cleaning in Section 3.2.

3.2 Halo catalogs with different exclusion criteria

Zhou & Han (2023) has demonstrated that the matter field can be decomposed into some self-similar halo distributions convolved by the Einasto profiles when halos are selected according to the $R_{\rm id}$. Considering the physical picture that the $R_{\rm id}$ separate the growing halo and the depleting environment (Fong & Han, 2021; Gao et al., 2023), it is natural to choose the $R_{\rm id}$ as the exclusion radius of halos. This choice is also supported by the fact that $R_{\rm id}$ was found to coincide with the optimal exclusion radius in a flexible halo model (García et al., 2021).

Nevertheless, as we will show in Section 4 below, it is still possible to build accurate halo models using other halo catalogs than the depletion catalog, at least for reproducing the halo-matter correlation as long as appropriate halo profiles are chosen.

To investigate these different solutions, we construct three halo catalogs by selecting halos according to different exclusion radii from the original catalog. These catalogs are named as

• $R_{\rm vir}$ catalog: The exclusion radius $R_{\rm ex}$ is defined as the virial radius;

• $R_{\rm id}$ catalog: The exclusion radius $R_{\rm ex}$ is defined as the inner depletion radius;

• 3$R_{\rm vir}$ catalog: The exclusion radius $R_{\rm ex}$ is defined as three times the virial radius.

For each catalog, we remove halos which intersect with a more massive neighbor. Specifically, if the distance between a halo pair is smaller than the sum of their exclusion radii, we remove the smaller one from the catalog. We define a exclusion scale of two halo populations $r_{\rm ex}=R_{\rm ex}(m_{1})+R_{\rm ex}(m_{2})$. Within $r_{\rm ex}$, there is no pair of halos with $m_{1}$ and $m_{2}$ due to exclusion. The remaining halos form a cleaned catalog with statistics such as the halo mass function and halo-halo correlation different from those in the original catalog. In practice, FoF halos hardly overlap in the virial region because the linking parameter is optimized for disentangling halos according to the virial radius. Thus, the resulting statistics of the original and $R_{\rm vir}$ catalogs are nearly identical.

Figure 1 illustrates how halos are distributed in different catalogs. In the left panel, all candidate halos are “isolated” and contained in the $R_{\rm vir}$ catalog since they do not overlap in their virial regions. As the exclusion radius increases in the right panels, more and more halos are removed from the catalog, and halos become more sparse. Correspondingly, it is natural to expect a steeper outer profile when halos are selected with a smaller exclusion criteria, and vice versa, to avoid double counting masses from halos.

Figure 1: The halo distributions in three catalogs. The projected map of a small region in the simulation box with a thickness of 3$h^{-1}{\rm Mpc}$ is shown, with halos marked by solid circles according to different boundary definitions from left to right. The dashed circles mark removed halos in each catalog according to the exclusion criteria. There are several solid circles overlapping because of the projection effect.

3.3 Measuring the power spectra

The auto and cross-power spectra are computed from the three-dimensional matter density field and halo fields in seven mass bins using pylians (Villaescusa-Navarro, 2018). We use the cloud-in-cell (CIC) scheme to assign particles to the mesh. The simulation box with side length $L=600h^{-1}\mathrm{Mpc}$ is divided into $N^{3}=1024^{3}$ voxels so that the smallest and largest $k$ modes accessible should be determined as $k_{\rm min}=2\pi/L$ and $k_{\rm max}=\pi N/L$ where $k_{\rm max}$ is estimated as the Nyquist frequency.

For a Gaussian field, the statistical error on the power spectra can be estimated as (Hamilton, 1997; Bernardeau et al., 2002)

$$\Delta P_{\rm aa} =\left(\frac{2}{N_{\rm mode}}\right)^{1/2}\left[P_{\rm aa}+P^{\rm shot}_{\rm a}\right],$$ (30)

$$\Delta P_{\rm ab} =\left(\frac{2}{N_{\rm mode}}\right)^{1/2}\sqrt{\frac{P^{2}_{\rm ab}+(P_{\rm aa}+P^{\rm shot}_{\rm a})(P_{\rm bb}+P^{\rm shot}_{\rm b})}{2}},$$ (31)

where $P_{\rm aa}$ and $P_{\rm ab}$ are the auto and cross power spectra. $P^{\rm shot}$ is the shot noise. $N_{\rm mode}(k)\simeq 2\pi/(k\sqrt{\Delta kV})$ is the number of $k$ modes in a bins centered at $k$, $\Delta k$ is the bin width, and $V$ is the volume of the box.

In practice, the raw power spectrum $P_{\rm r}(k)$ obtained from fast Fourier transform (FFT) is a smoothed version of the true power spectrum due to the window function effect from mass assignment scheme, and suffers from the aliasing effect due to discrete sampling of the field (Baugh & Efstathiou, 1994; Jing, 2005). To correct for these effects, we also compute the power spectrum from the Fourier transform of the correlation function,

$$P^{\prime}(k)=\int_{r_{\rm min}}^{r_{\rm max}}\xi(r)\frac{\sin(kr)}{kr}4\pi r^{2}\mathrm{d}r,$$ (32)

where $r_{\rm min}=0.02h^{-1}{\rm Mpc}$ and $r_{\rm max}=90h^{-1}{\rm Mpc}$. The correlation function $\xi(r)$ is directly measured from the simulation for $0.02\sim 20h^{-1}{\rm Mpc}$ and extrapolated to $90h^{-1}{\rm Mpc}$ using the linear correlation function computed with colossus (Diemer, 2018). With $r_{\rm min}$ much smaller than the grid size of the FFT, the Fourier transform of the correlation function is closer to the true power spectrum on small scale, but may lose some large-scale information due to our limited $r_{\rm max}$ in the transform. We thus combine the two to give our final measurement of the power spectrum as

$$P(k)=\left\{\begin{array}[]{rcl}P_{\rm r}(k)&&k<0.5h{\rm Mpc}^{-1}\\ P^{\prime}(k)&&k\geq 0.5h{\rm Mpc}^{-1}\end{array}\right.$$ (33)

Figure 2 shows the raw and corrected power spectra for $k<k_{\rm max}$. We find that the FFT estimation of the halo-matter power spectrum is significantly biased for $k>1h\mathrm{Mpc}^{-1}$. This is mostly due to the window function effect which smoothes out the halo-matter correlation on small scale. After correction, the power spectrum is higher. For the matter-matter power spectrum, the raw power spectrum is affected by numerical effects near the Nyquist frequency. We do not correct the halo-halo power spectra because the exclusion scale $k_{\rm ex}=1/r_{\rm ex}$ is smaller than $1h\mathrm{Mpc}^{-1}$ for most halo pairs, beyond which the halo-halo power spectrum is unimportant as we will discuss in Section 6.

Figure 2: The raw power spectra measured using pylians (Villaescusa-Navarro, 2018) (dotted curves) and the corrected power spectra from Equation (30) (solid curves). $P_{\rm hm}$ is the cross power spectrum of the matter field and the halo population with mass $10^{11.50}<(m/h^{-1}{\rm M}_{\odot})<10^{12.05}$.

To concentrate on the transition region we are interested in, in the following we will focus on the power spectra in the range of $0.1<k<3h\mathrm{Mpc}^{-1}$, which is a narrower range than $k_{\rm min}<k<k_{\rm max}$. This $k$ range corresponds to a radial range from several to about sixty $h^{-1}\mathrm{Mpc}$, containing information about halo profiles in the transition region.

García & Rozo (2019) and Zhou & Han (2023) have illustrated how the exclusion scheme impacts the statistics of halo populations in real space. Here, we show some statistics of different halo catalogs in Fourier space. Figure 3 shows the halo-matter power spectra $P_{\rm hm}$s and the halo-halo power spectra $P_{\rm hh}$s of three mass bins. For a given halo catalog, with increasing halo mass $P_{\rm hm}$ and $P_{\rm hh}$ are higher and $P_{\rm hh}$ truncates at a larger scale. This is because the halo bias $b$ and the exclusion scale $r_{\rm ex}$ are increasing functions with respect to halo mass $m$. Comparing the measurements from different catalogs for a given mass bin, $P_{\rm hm}$ and $P_{\rm hh}$ are lower with increasing exclusion radius because halos are less clustered in a catalog with larger exclusion radii. On small scales, $P_{\rm hm}$s converge because the exclusion scheme hardly affects the small-scale matter distribution around halo centers. $P_{\rm hh}$ truncates at larger scales with increasing exclusion radius. Across mass bins, we find the difference between $P_{\rm hm}$s of different catalogs decreases with increasing halo mass, since the exclusion scheme significantly affects the statistics of low-mass halos but the impacts become weak for high-mass halos.

Figure 3: Comparison of population statistics from different halo catalogs. Top: the halo-matter power spectra $P_{\rm hm}(k|m)$ for logarithmic mass bins $[11.50,12.05]$, $[12.60,13.15]$, and $[13.70,14.25]$. Bottom: halo-halo power spectra $P_{\rm hh}(k|m,m^{\prime})$. $m^{\prime}$ is in the logarithmic mass bin $[11.50,12.05]$ and $m$ are in the same mass bins as top plots.

4 Numerical solutions for the halo density profiles

4.1 Matching halo profiles in Fourier space

We measure the halo mass function, the halo-halo power spectra, and the halo-matter power spectra from the simulation, and then solve Equation (28) to obtain the halo profiles. As the large-scale statistics of halos depend on the exclusion criterion, the solution to Equation (28) varies with halo catalogs. We refer to these solutions as the matching profiles to each halo catalog.

Figure 4 shows the matching profiles in Fourier space. For the depletion catalog, Zhou & Han (2023) found that the Einasto profile works well for constructing a halo model in real space. We thus fit the halo profiles in the depletion catalog within $R_{\rm id}$ using the Einasto formula in real space and invert them into Fourier space as references.

On small scales of $k>1h\mathrm{Mpc}^{-1}$, the matching profiles follow the corresponding Einasto profiles well in all catalogs for halos with mass $M>10^{12.05}h^{-1}M_{\odot}$. This is consistent with the fact that the inner profiles of halos are barely affected by the exclusion criteria.

On scales $k<1h\mathrm{Mpc}^{-1}$, the profiles are more complicated. As $k\rightarrow 0$, $W(k)$ approaches the integrated mass of the density profile. It is thus expected that this asymptotic amplitude of $W(k)$ will increase as a larger exclusion radius is adopted, which is indeed the case in Figure 4 except for the lowest mass bin. In the $R_{\rm vir}$ catalog, the optimal profiles are generally lower than the Einasto fits, while they are higher than the references in the 3$R_{\rm vir}$ catalog. By contrast, the profiles in the $R_{\rm id}$ catalog still follow the Einasto model well. For each profile, we have also shown the integrated mass up to the corresponding halo radius in each catalog, which is lower than the integrated mass as the profiles all extend beyond the exclusion radii.

In the lowest-mass bin, however, the matching profiles no longer follow the above expectations. In the $R_{\rm vir}$ catalog, the matching profile of the lowest-mass bin is significantly higher than the Einasto fit and has not converged to a constant towards lower $k$. This can be largely attributed to improper modeling of the unresolved term (see Appendix A) in this catalog when solving for the matching profile. The lowest mass bin has masses closest to the unresolved halos and are thus the most degenerate with the unresolved component. We discuss the influence of the unresolved component further in section 6.

Another notable feature is that there are some fluctuations on scales $k<1h\mathrm{Mpc}^{-1}$, especially in low-mass bins. According to Equation (30), the measurements of the power spectrum on larger scales have higher uncertainties. These uncertainties can be propagated to the solution via Equation (28). Low mass halos have a lower $W(k)$, resulting in a much larger relative fluctuation.

In Figure 5 we also compare the matching profiles in the $R_{\rm id}$ catalog to the NFW profiles, which are more extended than the Einasto profiles and do not converge to a finite mass on large scale. It can be seen that the NFW fits agrees with the matching profiles on small scale, but over-predicts them on large scale.

Figure 4: The matching profiles in Fourier space for the three catalogs. Solid curves are the profiles solved from Equation 28, while the dotted curves are the Einasto fits to the real-space profiles in the $R_{\rm id}$ catalog (identical across three panels). Different colors correspond to different mass bins. For each profile, the arrow on the left marks the total mass enclosed in its exclusion radius (i.e., $M_{\rm vir}$, $M_{\rm id}$ or $M({\rm<3R_{\rm vir}}$)).

Figure 5: Comparison of the NFW profile, Einasto profile, and matching profile for the $R_{\rm id}$ catalog. The NFW profiles are fitted in real-space up to $R_{\rm id}$ and then transformed to $k$-space.

We verify that the solutions can accurately reproduce the halo-matter power spectra by inserting the matching profiles back into Equation (28), and found residuals less than $10^{-12}$. Thus we conclude that, for a given halo catalog, there is indeed a set of matching profiles to accurately reconstruct the mass distribution around halos, at least in the halo-matter correlation.

4.2 Matching profiles in real space

To gain more physical insights into halo profiles, we transform the matching profiles into real space. As our Fourier space solutions only cover the scales $0.1<k<3h\mathrm{Mpc}^{-1}$, extrapolations to both the high-$k$ and low-$k$ ends are needed before we can transform them to real space. According to the behaviors discussed above, we extrapolate the profiles using the Einasto fits on scales $k>3h\mathrm{Mpc}^{-1}$, and using constant values of $W(k=0.1h\mathrm{Mpc}^{-1})$ for $k<0.1h\mathrm{Mpc}^{-1}$. The inverse Fourier transform (IFT) is evaluated over the $k$ range of $0<k<10^{3}h\mathrm{Mpc}^{-1}$. Note that there are some potential numerical effects in the IFT. For instance, the IFT results can be affected by the missing high-$k$ modes due to the finite $k$ range used, resulting in noises at high frequency. In addition, the interpolation of k-space profiles can lead to false signals in the IFT results. To reduce noises caused by these effects, the final results are further smoothed using the Savitzky–Golay algorithm. We also invert the $k$-space Einasto profiles to real space by integrating over the same $k$ range to avoid numerical effects in the comparison.

Figure 6: The matching profiles for the three catalogs in real space. In each panel, the solid and dashed curves are the matching profiles and the reference Einasto profiles in real space, respectively. The upper triangles, lower triangles, and stars mark the locations of $R_{\rm vir}$, $R_{\rm id}$, and 5$R_{\rm vir}$, respectively.

Figure 6 shows the matching profiles in real space. For the $R_{\rm id}$ catalog, the profiles follow the Einasto profiles within 5$R_{\rm vir}$. On scales $r>5R_{\rm vir}$, the matching profiles deviate from the Einasto profiles, and are affected by numerical perturbations. In the $R_{\rm vir}$ catalog, the matching profiles follow the Einasto profiles within $R_{\rm id}$ for all mass bins. Outside $R_{\rm id}$ they drop faster than the Einasto profiles and are also perturbed due to numerical effects. The exception comes from the lowest-mass bin whose profiles are higher than the Einasto profiles on scales $r>R_{\rm id}$. Similarly, the matching profiles to the 3$R_{\rm vir}$ catalog follow the Einasto profiles within $R_{\rm id}$, but become higher outside $R_{\rm id}$ except for the lowest mass bin. These behavior are all consistent with our previous discussions in frequency space.

Combining the three panels, it is very interesting to notice that the matching profiles inside $R_{\rm id}$ are nearly identical across the catalogs, while the outer profiles become more extended with increasing exclusion radii.

In the $R_{\rm id}$ catalog, the profiles follow the Einasto fits well out to a very large scale. This result agrees with the findings of Zhou & Han (2023), in which the halo-matter correlation functions are accurately modeled when coupling Einasto profiles to the $R_{\rm id}$ catalog. In this work, we directly solve the profile matched with the $R_{\rm id}$ catalog and find that the optimal profile can be analytically modeled with the Einasto formula, which automatically explains the choice of the halo profile in the Depletion-radius-based Halo Model (here after DHM) of Zhou & Han (2023). The deviations outside the 5$R_{\rm vir}$ do not challenge this conclusion since the profile on very large scales is unimportant compared with the 2-halo term.

5 Towards an explicitly verifiable and versatile halo model for the large scale structure

The numerical solutions in section 4 illustrate that a matching halo profile exists for each halo catalog in accurately reproducing the halo-matter power spectrum. The remaining question is then whether and how these solutions may be used to construct accurate analytical halo models. To this end, it is desirable to have analytical prescriptions for each of the model components involved, including the halo profile, halo-halo correlation, halo mass function and the unresolved component which were mostly handled numerically rather than analytically in section 4. Establishing analytical prescriptions for these components is not always straightforward due to the complexity and non-universality of the components in some catalogs.

For the $R_{\rm id}$ catalog, however, we have showed that the Einasto profile is a good analytical model for the matching profile. In fact, such a profile has already been adopted in Zhou & Han (2023) when constructing the real-space halo model based on the depletion catalog, along with analytical recipes for all other components. It is thus straight-forward to translate the model of Zhou & Han (2023) to Fourier space for the prediction of the power spectrum. Now, we proceed to check how well this model works.

To construct the model, the profile parameters and the statistics of halos are generated using the fitting formula in Zhou & Han (2023). The Einasto profiles, the halo-halo correlations, and the unresolved halo-matter correlations are then converted to Fourier space and inserted into Equation (28) to predict the halo-matter power spectra. The predictions are shown in the left panel of Figure 7. The error bars show the standard errors estimated using Equation (31). We find that the DHM performs well in predicting the halo-matter power spectra, achieving about 10$\%$ accuracy (within statistical errors). This performance agrees with that of DHM in real space.

Ideally, a perfect halo model can fully match the entire density field of the universe, thus capable of predicting multiple statistics of the halo and matter field simultaneously. In Figure 7 we show that the same DHM can also accurately predict the total matter power spectrum. Following the idea in Section 2.2, we modify the classical halo model to include the contribution from the unresolved mass, and compute the matter-matter power spectrum as detailed in Appendix B. One can see that the DHM achieves 5$\%$ accuracy in a $k$ range that covers the transition of 1-halo and 2-halo terms, by directly summing up the various components.

For comparison, Figure 8 shows the results from three other models. First, we construct a “classical” halo model adopting commonly-used recipes for its components, including a 2-halo term in proportion to the linear power spectrum, a halo mass function fitted with Sheth & Tormen (2002), and an Einasto profile truncated at the virial radius with profile parameters specified in Diemer & Joyce (2019) and Gao et al. (2008). We also compare the DHM with some fine-tuned models for predicting the power spectrum, including halofit (Smith et al., 2003; Takahashi et al., 2012) and hmcode (Mead et al., 2015, 2016, 2021).

Focusing on the classical model, we find it performs poorly in this region. This under-prediction suggests a defect in this model, which is caused by the virial-truncated profile and a simplistic 2-halo term. According to the results in section 4, a matching profile is extended and have a smooth truncation. Using a virial-truncated profile will miss the contribution of the halo outskirts, leading to the underprediction of the matter-matter power spectrum. In addition, the 2-halo term also deviates from the linearly-biased linear power on small scale. On the other hand, DHM not only uses an extended Einasto profile to accurately model the mass distribution in the halo outskirts, but also carefully considers the halo exclusion and unresolved mass in the 2-halo term, improving the accuracy on transition scales.

Compared with halofit and hmcode, DHM reaches a comparable accuracy. However, it is important to realize that the former two models are “implicit” or pseudo halo models, while DHM is an explicit and physical halo model. In DHM, a halo population based on a physical boundary definition is explicitly identified and modeled, which facilitates further studies on the corresponding halos. By contrast, halofit is a fitting method without involving any halo property, and hmcode pertains to a population of “effective halos” whose properties differ from those measured in simulations. Moreover, DHM is more interpretable. The model parameters in DHM have clear physical meanings and can all be measured from the halo population except for one parameter of the unresolved component which can still be derived from first principle. On the other hand, the parameters in halofit and hmcode are less interpretable because they do not correspond to any physical halo population. Finally, the physical nature of DHM enables it to predict multiple statistics simultaneously, including (but not limited to) the total matter and halo-matter power spectra as we have shown. By contrast, it is challenging to extend halofit and hmcode to predict the power spectra of other tracers because they use some ad-hoc fixes in modeling the matter-matter power spectrum.

Figure 7: Left:The halo-matter power spectra from the simulation and the DHM. The errors are estimated using Equation (31). The bottom panel shows the ratios between the simulated and model power spectra. The shaded areas mark the $10\%$ difference level. right: The matter-matter power spectrum from the simulation and the DHM. The red curve with error bars are simulated power spectrum. The model prediction is decomposed into different components including the resolved 1-halo and 2-halo terms, the halo-unresolved (hu), unresolve-unresolved (uu) and the Poisson terms as detailed in Appendix B. The bottom panel shows the ratios between the simulated and model power spectra. The shaded areas mark the $5\%$ difference level.

Figure 8: Ratios between the simulated and predicted power spectra using different models. Red curve is computed using virial-truncated profiles and a linearly-biased 2-halo term, as detailed in Section 5. Orange, purple and blue curves are computed using halofit (the version of Takahashi et al., 2012), hmcode, and DHM. Error bars represent the standard errors estimated using Equation (30). The shaded areas mark the fractional errors less than 0.05.

6 Discussion: profile uncertainty due to the unresolved term

The modeling of the unresolved term relies on an extrapolation of the halo-halo correlation to the diffuse mass limit, in addition to a model for the integrated mass over the halo profile before solving for the matching profile. In this work, we have modeled these components following the treatment of Zhou & Han (2023), which is consistent with the depletion catalog but may not work for the other catalogs.

To show the influence of this modeling uncertainty, in Figure 9 we compare the solutions to Equation (28) with and without the unresolved term for the depletion catalog. The solutions are barely changed for the high mass bins, while that for the lowest mass bin is significantly affected. This is not difficult to understand as the unresolved term should be most degenerate with halos close to being unresolved. Without the unresolved term, the profile of the lowest mass bin is significantly overestimated, to compensate for the missing mass and power. Similar overestimation of the matching profile can be observed in Figure 4 for the lowest mass bin in the $R_{\rm vir}$ catalog, due to an underestimated unresolved term as detailed in Appendix A.

Figure 9: Same as Figure 4, but comparing the profile solutions of the $R_{\rm id}$ catalog with and without the unresolved term, $P^{\rm u}_{\rm hm}(k)$, in Equation (28), as labeled.

The significance of the unresolved term can also be understood from its relative contribution to the power spectrum in each halo mass bin, as shown in Figure 10. For the highest mass bin, the power spectrum is dominated by the 1-halo term over the entire $k$-range covered. For lower mass halos, the 2-halo term always dominates on the large scale. As the halo mass decrease, the contribution from the unresolved term increases, while the 1-halo term becomes less and less important. It is thus natural to expect that a variation in the unresolved term will have a significant influence on the halo profile of the low mass halos, while it barely affects that of high mass halos.

Figure 10: Decompositions of the halo-matter power spectra for the $R_{\rm id}$ catalog. In each panel, the thick grey curve shows the measured power spectrum for a given mass bin as labeled, while the blue solid curve shows the model prediction, with dashed curves showing the contributions from the 1-halo, 2-halo and unresolved terms respectively. The black solid curve shows the Einasto model as in Figure 4.

7 Summary and conclusions

In this work, we solve for the halo density profiles in Fourier space from measurements of the halo-matter and halo-halo power spectra in a cosmological simulation. This is possible because the matter field can be modeled by the halo field convolved with the halo density profiles. The Fourier space solution thus avoids ambiguities in partitioning matter among neighboring halos, enabling us to derive the complete halo profile out to large scale.

We have applied the method to three halo catalogs with different boundary definitions, including the virial radius, $R_{\rm vir}$, the depletion radius, $R_{\rm id}$, and $3R_{\rm vir}$, to study the relations between halo profile and boundary definition, and their implications for halo model. Our main findings are as follows.

• For each of the halo catalogs, a set of matching halo profiles can always be found to accurately reconstruct the input halo-matter power spectra. This implies there can be multiple ways of decomposing the matter field into halos, and the matter distribution around the halos can always be well described with a halo model, at least in terms of the halo-matter power spectrum.

• The matching profiles vary in different halo catalogs due to different population properties. Halos defined with a more extended boundary also have more extended outer profiles and larger integrated masses over the profiles. The profiles all extend beyond the exclusion radius used to define the halo catalogs, reflecting ongoing mass accretion and mergers which extend the mass distribution around halos.

• The matching profiles are nearly identical within $R_{\rm id}$ across the three catalogs. This indicates that the depletion region characterized by $R_{\rm id}$ is a universal feature around halos.

• For the $R_{\rm id}$ catalog, the matching profiles are well described by the Einasto profiles over the scale covered by our measurements. This supports the choice of the Einasto profile in the depletion halo model of Zhou & Han (2023). For the $R_{\rm vir}$ catalog, the matching profiles drop faster than the Einasto profile outside the depletion radius.

• Proper modeling of unresolved halos and diffuse matter is important for extracting the profiles around low mass halos. Ignoring or underestimating the unresolved term can result in overestimation of the low mass profiles.

The agreement between the numerical solution and the Einasto profile model for the depletion catalog enables us to construct an analytical halo model for this catalog, which is the Fourier space equivalence to the depletion halo model in Zhou & Han (2023). This model can predict the halo-matter power spectra to 10% accuracy across scales and halo masses. Moreover, the same model also predicts the matter-matter power spectrum to 5% accuracy. This illustrates the advantage of an explicit and physical halo model, that the same model can simultaneously predict multiple statistics of the halo and matter fields. For comparison, existing halo models for predicting the matter power spectrum such as halofit and hmcode can be regarded as implicit or effective models. These models do not directly correspond to an identifiable halo population, and therefore the model ingredients can not be directly verified with simulations. As a result, the predicting power of these implicit models is usually limited to the total matter power spectrum which is used for training the model.

In this work, we have only constructed and verified the DHM for a single snapshot at $z=0$. In future works, we will extend our model to different redshifts and cosmologies, by studying the redshift evolution and cosmology dependence of the model ingredients. This may eventually lead to a precise, unified and physical understanding of structure formation using halos as building blocks.

Appendix A Modeling the unresolved term

According to Zhou & Han (2023), the halo-halo correlation functions follow a universal shape. The unresolved term can be modeled by generalizing the universal halo correlation down to the diffuse matter limit, expressed as

$$\xi^{\rm unr}_{\rm hm}(r|M)=b_{\rm unr}b_{\rm hh}(M)\hat{\xi}_{\rm hh}(r)\mathrm{H}[r-R_{\rm id}(M)],$$ (A1)

where $\hat{\xi}_{\rm hh}(r)$ is the universal halo correlation that can be modeled by a power-law, and ${\rm H}(x)$ is the Heaviside step function which is unity at $x>0$ and 0 otherwise. The effective bias $b_{\rm unr}$ can be estimated using local mass conservation,

$$b_{\rm unr}\approx 1-\frac{1}{\bar{\rho}_{\rm m}}\int^{m_{\rm max}}_{m_{\rm res}}M_{\rm int}n(m)b_{\rm hh}(m)\mathrm{d}m,$$ (A2)

where $m_{\rm max}$ and $m_{\rm res}$ are the mass limits of resolved halos, and $M_{\rm int}=\int\rho(\textbf{{r}})\mathrm{d}^{3}\textbf{{r}}$ is the mass integral of the halo profile, which is not necessarily equal to a general mass label $m$. Estimating $b_{\rm unr}$ according to Equation (A1) requires knowing the halo profile a priori. For the $R_{\rm id}$ catalog, we have shown that the Einasto formula can model the halo profile up to the transition scale, so that $b_{\rm unr}$ of the $R_{\rm id}$ catalog can be directly estimated. For the $R_{\rm vir}$ and $3R_{\rm vir}$ catalogs, the Einasto formula performs poorly in modeling the outer profiles, as shown in Section 4.2. Despite this, we still adopt the integrated mass of the Einasto profile when evaluating Equation (A1) for all three catalogs for simplicity. Meanwhile, the halo bias and mass function are fitted separately using the parametric functions in Zhou & Han (2023) for each catalog, as these quantities vary in different catalogs. As a result, we expect $b_{\rm unr}$ to be underestimated for the $R_{\rm vir}$ catalog but overestimated for the $3R_{\rm vir}$ catalog, due to over/under estimation of the integrated mass in the two catalogs respectively.

Appendix B Fourier version of the depletion-radius-based halo model

To derive the matter-matter power spectrum with the DHM, we write the matter overdensity field $\delta_{\rm m}(\textit{{x}})$ as

$$\delta_{\rm m}(\textit{{x}})=\frac{1}{\bar{\rho}_{\rm m}}\sum_{j}n(m_{j})\delta_{\rm h}(\textit{{x}}|m_{j})\circledast\rho(\textit{{x}}|m_{j})+\delta^{\rm u}_{\rm m}(\textit{{x}}),$$ (B1)

where $n(m)$, $\delta_{\rm h}(\textit{{x}}|m)$, and $\rho(\textit{{x}}|m)$ are the mean number density, the overdensity field and the profile of a halo population with mass $m$, respectively, and $\circledast$ is the convolution operation. The term $\delta^{\rm u}_{\rm m}(\textit{{x}})$ represents the contribution of the unresolved mass. Considering the Fourier transform of Equation (B1), we have

$$\delta_{\rm m}(\textit{{k}})=\frac{1}{\bar{\rho}_{\rm m}}\sum_{j}n(m_{j})\delta_{\rm h}(\textit{{k}}|m_{j})W(\textit{{k}}|m_{j})+\delta^{\rm u}_{\rm m}(\textit{{k}}).$$ (B2)

The matter-matter power spectrum is

$$P_{\rm mm}(\textit{{k}}) =\frac{1}{\bar{\rho}_{\rm m}^{2}}\sum_{i}\sum_{j}n(m_{i})n(m_{j})W(\textit{{k}}|m_{i})W(\textit{{k}}|m_{j})P_{\rm hh}(\textit{{k}}|m_{i},m_{j})$$

$$+\frac{2}{\bar{\rho}_{\rm m}}\sum_{j}n(m_{j})W(\textit{{k}}|m_{j})P_{\rm hu}(\textit{{k}}|m_{j})$$

$$+P_{\rm uu}(\textit{{k}}),$$ (B3)

where $P_{\rm hu}(\textit{{k}}|m_{j})$ and $P_{\rm uu}(\textit{{k}})$ are the cross power spectra of halo-unresolved and unresolved-unresolved.

The 1-halo term arises from $P_{\rm hh}(\textit{{k}}|m_{j})$ when two halos are identical. For unresolved low mass halos, we can approximate them as point masses with $W(k|m\rightarrow 0)\approx M_{\rm int}$. We thus break the 1-halo term into resolved and unresolved halos as

$$P^{\rm 1h}_{\rm mm}(k)=\frac{1}{\bar{\rho}_{\rm m}^{2}}\sum_{j}^{\rm resolved}n(m_{j})W^{2}(k|m_{j})+\frac{1}{\bar{\rho}_{\rm m}^{2}}\int_{m_{\rm fs}}^{m_{\rm res}}M_{\rm int}^{2}n(m)\mathrm{d}m,$$ (B4)

where $m_{\rm fs}$ and $m_{\rm res}$ are the freestreaming mass and the minimum mass of resolved halos considered by the model. The second term represents the shot-noise from point-mass halos. For the power spectrum measured from simulations, the diffuse matter are resolved by discrete particles and thus should also contribute a shot noise term. In the following we will not distinguish the two contributions but model them with a single constant $C$.

The 2-halo term involves contributions from the halo-halo, halo-unresolved, and unresolved-unresolved terms, expressed as

$$P^{\rm 2h}_{\rm mm}(k) =\frac{1}{\bar{\rho}_{\rm m}^{2}}\sum_{i}\sum_{j\neq i}n(m_{i})n(m_{j})W(k|m_{i})W(k|m_{j})P^{\rm 2h}_{\rm hh}(k|m_{i},m_{j})$$

$$+\frac{2}{\bar{\rho}_{\rm m}}\sum_{j}n(m_{j})W(k|m_{i})P^{\rm 2h}_{\rm hu}(k|m_{j})$$

$$+P^{\rm 2h}_{\rm uu}(k).$$ (B5)

Note here the unresolved terms are contributed by both unresolved halos and diffuse matter. The power spectra $P^{\rm 2h}_{\rm hu}(k|m_{j})$ and $P^{\rm 2h}_{\rm uu}(k)$ can be obtained from Fourier transforms of the correlation functions specified in Zhou & Han (2023),

$$\xi^{\rm 2h}_{\rm hu}(r|m_{j}) =b_{\rm unr}b(m_{j})\hat{\xi}_{\rm hh}(r){\rm H}[r-R_{\rm ex}(m_{j})],$$ (B6)

$$\xi^{\rm 2h}_{\rm uu}(r) =b^{2}_{\rm unr}\hat{\xi}_{\rm hh}(r).$$ (B7)

Replacing the summations with integrations, the final expression of the matter-matter power spectrum in the depletion halo model is

$$P_{\rm mm}(k) =P^{\rm 1h}_{\rm mm}(k)+P^{\rm 2h}_{\rm mm}(k)$$ (B8)

$$P^{\rm 1h}_{\rm mm}(k) =\frac{1}{\bar{\rho}_{\rm m}^{2}}\int_{m_{\rm res}}^{m_{\rm max}}n(m)W^{2}(k|m)\mathrm{d}m+C$$ (B9)

$$P^{\rm 2h}_{\rm mm}(k) =\frac{1}{\bar{\rho}_{\rm m}^{2}}\int_{m_{\rm res}}^{m_{\rm max}}\int_{m_{\rm res}}^{m_{\rm max}}n(m_{1})n(m_{2})W(k|m_{1})W(k|m_{2})P^{\rm 2h}_{\rm hh}(k|m_{1},m_{2})\mathrm{d}m_{1}\mathrm{d}m_{2}$$

$$+\frac{2}{\bar{\rho}_{\rm m}}\int_{m_{\rm res}}^{m_{\rm max}}n(m)W(k|m)P^{\rm 2h}_{\rm hu}(k|m)\mathrm{d}m$$

$$+P^{\rm 2h}_{\rm uu}(k).$$ (B10)

where $m_{\rm res}$ and $m_{\rm max}$ represent the mass limits of resolved halos.

An alternative derivation of the above equations can be obtained from Fourier transforming the matter-matter correlation,

$$\xi_{\rm mm}(r)=\int f(m)\xi_{\rm hm}(r|m)\circledast u_{\rm h}(r|m)\mathrm{d}m,$$ (B11)

where $\xi_{\rm hm}(r|m)$ is modeled in Zhou & Han (2023) and the integration continues down to the diffuse matter limit.

In the above model, we have not considered the existence of subhalos inside each halo, which could contribute additional small-scale power to the 1-halo term on top of the smooth density profile (Hezaveh et al., 2016; Diaz Rivero et al., 2018; Díaz Rivero et al., 2018). A complete treatment will have to involve the distribution and internal structure of subhalos (Han et al., 2016; He et al., 2023). For scales much larger than the sizes of subhalos, however, we can approximate subhalos as point masses, so that their contributions can also be approximated as a Poisson term, which can be absorbed into the parameter $C$ above. As shown in Section 5, the above model performs well for $0.1<k<3h\mathrm{Mpc}^{-1}$ with 5$\%$ accuracy, with a best-fitting parameter $C=16.04$ for our adopted $m_{\rm res}=10^{6}h^{-1}{\rm M}_{\odot}$ and $m_{\rm max}=10^{16}h^{-1}{\rm M}_{\odot}$. Nevertheless, one should keep in mind that the matter-matter power spectrum does not converge to a constant with $k$ increasing, and a scale-dependent correction term is needed if a higher precision is required. We will extend our model with subhalo terms in future works.