CHAM: a fast algorithm of modelling non-linear matter power spectrum in the sCreened HAlo Model
Bin Hu Xue-Wen Liu Rong-Gen Cai
astro-ph.CO

















    ł/γμν∇_μ∇_νℋ𝒵 et al. i.e.  e.g. etc.bhu@bnu.edu.cn^1Department of Astronomy, Beijing Normal University, Beijing, 100875, China
        ^2CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
        
        Chinese Academy of Sciences, P.O. Box 2735, Beijing 100190, China 
        ^3School of Physical Sciences, University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing 100049, China
        We present a fast numerical screened halo model algorithm () for modeling non-linear power spectrum for
        the alternative models to ΛCDM. This method has three obvious advantages.
        First of all, it is not being restricted to a specific dark energy/modified gravity model.
        In principle, all of the screened scalar-tensor theories can be applied.
        Second, the least assumptions are made in the calculation.
        Hence, the physical picture is very easily understandable.
        Third, it is very predictable and does not rely on the calibration from N-body simulation.
        As an example, we show the case of Hu-Sawicki f(R) gravity.
        In this case, the typical CPU time with the current parallel Python script (8 threads) is roughly within 10 minutes.
        The resulting spectra are in a good agreement with N-body
        data within a few percentage accuracy up to k∼1 h/ Mpc.
        CHAM: a fast algorithm of modelling non-linear matter power spectrum in the sCreened HAlo Model
            Rong-Gen Cai^2,3
            December 30, 2023
        ===============================================================================================
        
        
        
        
        
        
        Several up-coming large scale structure surveys, such as Euclid [<http://sci.esa.int/euclid>], LSST [<http://www.lsst.org>], WFIRST [<https://wfirst.gsfc.nasa.gov>], are aiming to measure the matter power spectrum range  from 0.1 to 10  Mpc/h up to 1% accuracy. One of the theoretical obstacles is how to model the non-linear power spectrum on these scales. By means of high precision simulation, we are able to
        do this modelling with the requested accuracy within the ΛCDM paradigm. However, another scientific issue associated with  these activities, is to understand the nature of late-time cosmic acceleration, which is normally interpreted by dark energy (DE) or modified gravity (MG) models.
        
        A number of alternative models to the ΛCDM  have inflated over the past decades.
        The predictions on the background expansion, from most of the viable DE/MG models, are hard to be distinguished from the standard scenario.
        Motivated by this fact, one of methods, is to focus on the parametrization of perturbation dynamics with given background cosmology.
        For the linear perturbation, the effective field theory approach <cit.>, provides us a uniformed parametrization for scalar-tensor
        type DE/MG models. The resulted linear Einstein-Boltzmann solvers, such as  <cit.>[<http://www.eftcamb.org/>] are being
        developed extensively.
        
        As for the non-linear part, several different methods are being studied, such as N-body simulation <cit.> and the hybrid scheme of the Lagrangian perturbation theory with N-body simulation <cit.>.
        Take N-body simulation as an example, although this approach could provide an accurate result, it is too much expensive to construct the template for all the DE/MG models.
        On the other hand, based on the current observations, the fifth force has to be shielded on the small scales, such as the solar system.
        This results in the fact that, most of DE/MG models can be categorised into a few types, according to different screening mechanism, such as chameleon <cit.>, Vainshtein mechanism <cit.>, etc.
        
        This motives us to unify the non-linear spectrum modelling via the screening mechanism. Recently, Ref. <cit.> proposed a generic parametrization of the modified gravity.
        In this parametrization, all the modified gravity effects are encoded into the scale-dependent gravitational constant
         
            G_ eff/G_ N   =    A + ∑_i^N_0 B_i∏_j^N_i b_ij
                   ·(r/r_0ij)^a_ij{[1+(r_0ij/r)^a_ij]^1/b_ij - 1 } ,
        
        where A, B_i, r_0ij, a_ij, b_ij are the screening parameters.
        With this great simplification, in this paper, we are aiming to present a fast and reasonably accurate algorithm for the non-linear spectrum in the DE/MG models.
        
        
        
        In details, our method is the extended halo model <cit.> with the scale-dependent gravitational constant given by Eq. (<ref>).
        The flow chart of our algorithm (namely, , denotes for sCreened HAlo Model) is presented in Fig.<ref>.
        In the literature, the halo model for the DE/MG models has been extensively studied, such as <cit.> for f(R) gravity. 
        Compared with them, our method has three advantages. First of all, this algorithm is not being restricted to a specific model.
        Basically, all of the screened scalar-tensor theories can be applied. Secondly, the least assumptions are made in the calculation.
        For example, instead of using the conventionally Sheth-Tormen mass function <cit.>, we solved the halo distribution function with a moving
        barrier in terms of the excursion set formalism <cit.>.
        In particular, the scale dependence of critical density, δ_c(r), is attributed not only to the MG effect in the process of spherical collapse, but
        also to the ellipsoidal collapsing effect <cit.>.
        Last but not the least, our model is very predictive due to the fact that we have very limited parameters.
        
        
        
        
        
        As a demonstration, we consider Hu-Sawicki f(R) gravity model <cit.>, which can satisfy the background ΛCDM expansion history
        and evade the solar system tests:
        
            f(R)=-m̅^2c_1(R/m̅^2)^n/c_2(R/m̅^2)^n+1, m̅^2≡κ^2ρ̅_m0/3,
        
        where m̅ refers to the present Compton mass of the extra scalar field.
        In the quasistatic but linear regime, the Poisson equation reads
        
            k^2 Ψ = -4 π G_ N( 4  3 - 1  31  k^2λ_C^2 + 1) a^2  δρ_ m .
        
        It is straightforward to see that below the Compton wavelength (kλ_C≫1) of the extra scalar field, the gravitational constant
        is enhanced by a factor 4/3. Above this scale, the General Relativity is recovered.
        
        In the non-linear high density regime, the fifth force carried by the scalar field is shielded by the chameleon mechanism.
        Considering a spherical over density regime, this mechanism is very similar to the static electrodynamics phenomenon.
        The scalar charge is only distributed on the surface, hence the thinner the surface is, the more significantly
        the fifth force is screened
        
            G_ eff/G_ N≈ 1 + 1/3min(3x-3x^2+x^3,1),  x≡Δ R/R_ TH ,
        
        where x denotes the surface thickness, which can be parametrized schematically as <cit.>
            Δ R/R_ TH≈ - C_1 r [(1+C_2 r^-3)^1/(α-1)-(1+C_3r^-3)^1/(α-1)],
        
        where C_1,C_2,C_3 and α can be read from the model parameters.
        
        The above description can be concluded by Fig.<ref>.
        Below the averaged screening scale (r_01) the MG effect is shielded.
        Above the Compton wavelength (r_02) the General Relativity is also recovered.
        Between these two scales, the gravitational constant is enhanced by a factor 4/3.
        
        
        
        
        
        
        With the above model setup, the process of spherical collapse is accordingly modified by the scale dependent
        gravitational constant. For simplicity, we firstly study the top-hat spherical collapse.
        The Lagrangian radius can be solved in terms of  y≡ (r/r_i - a/a_i) <cit.>
            y”+H'/Hy'=-1/2Ω_ma^-3-2Ω_Λ/Ω_ma^-3+Ω_Λy
             -1/2Ω_ma^-3/Ω_ma^-3+Ω_Λ·
                  G_ eff(r)/G_N(a/a_i+y)[(1/ya_i/a+1)^3(1+δ_i)-1] ,
        
        where the prime is the derivative with respect to log a and a_i, r_i, δ_i denote the initial time, radius
        and density of the corresponding Euclidean regime.
        
        By the energy conservation law, as long as the initial (negative) potential energy dominates over the (positive) kinetic energy,
        the initial over density patch, will always decouple from the background Hubble flow and finally collapse into a
        virialized object. The collapsing threshold is proxied by the extrapolated linear density δ_c. In this work, we adjust
        the initial value by asking the over density patch collapse at a=1. Hence, we integrate the following linear density
        equation to the present time
         
            δ”_m(r,a) + [2-3/2Ω_m(a)]δ'_m-3/2G_ eff(r)/G_N·Ω_m(a)δ_m =0,
        
        with G_ eff(r)/G_N described by Fig.<ref>. 
        We shall emphasize that unlike the ΛCDM case, the critical density in the modified gravity models is generally scale
        dependent due to the scale dependent gravitational constant. This can also be seen from Fig.<ref> (right vertical axis).
        When the gravitational constant restores the Newtonian value, the spherical collapse critical density recovers δ_sc=1.676.
        When G_ eff is enhanced, δ_sc reaches 1.692.
        
        
        
        
        
        
        Following the flow chart Fig.<ref>, we use the excursion set formalism <cit.> to compute the probability f(S) of forming a virialized object
        with given linear matter fluctuation
        
            S(r)≡∫ d^3k|W̃(kr)|^2P_L(k) .
        
        At this step, let us consider a more realistic model, such as the ellipsoidal collapsing process. 
        It will also introduce a scale dependent
        critical density, such as <cit.>
            δ_c = √(a)δ_sc[1+β (aν)^-α],  ν≡δ^2_sc/S,
        
        where α=0.615,β=0.485 and δ_sc is given by the spherical MG collapse.
        Here we assumed that the effects of modified gravity and ellipsoidal collapse can be treated separately.
        With the above moving barrier, the solution of distribution probability is given by <cit.>
            f(S)    =    g(S) + ∫^S_0 dS'f(S')h(S,S'),
        
        in which
        
            g(S)    ≡   [δ_c/S-2dδ_c/dS]P(δ_c,S),
            
            h(S,S')    ≡   [2dδ_c/dS-δ_c-δ'_c/S-S']P(δ_c-δ'_c,S-S'),
        
        where P(δ, S) is the normalized Gaussian distribution. This equation could be integrated numerically on an
        mesh with equal spacing on: S_i=iΔ S with i=0,1,⋯,N and Δ S=S/N.
        
        
        
        
        
        
        
        Armed with the distribution probability f(ν), we present the halo mass function in Fig.<ref>,
        where the data points are from N-body simulation in the literature <cit.>.
        We can see that  prediction (blue) is systematically better than the Sheth-Tormen formalism (red),
        in particular, in the low mass range. As for the high mass range, the two prescriptions behave statistically similar due to the
        large scatters of the N-body data. 
        
        The linear bias in  is derived by using the peak-background split approach. In the large cell
        limit, the mass inside a cell M is much larger than the typical halo mass m, this leads to the number of haloes with mass m
        inside the mass cell M can be approximated by
        
            b_ L(M_vir)   =    1-∂ln n/∂δ_sc.
        
        The derivative is realised by using Richardson four step interpolation to narrow the numerical error.
        
        As for the density profile, we use the NFW form <cit.>
            ρ(r) = ρ_s/(r/r_s)(1+r/r_s)^2,
        
        where r_s and ρ_s are characteristic radius and density, which can be parametrized via the
        concentration parameter
         <cit.>
            c_vir(M_vir)     =     9 [ M_*(M_vir)/M_vir]^0.13,
        
        where M_* is defined via σ(M_*)=δ_sc. We shall emphasize again that the scale
        dependence of δ_c leads to a different characteristic mass M_* at different scales.
        
        
        
        Now, we can assemble the mass function, linear bias as well as the concentration into the spectrum.
        In the halo model, all mass is within individual halos. Hence, the correlation function is made out of two terms. One 
        is the halo-halo correlation function which describes the density correlation between two halos on the large scale. 
        The other is one halo term which describes the density correlation between two points inside one halo. 
        This term shall dominate on the small scale. In Fourier space, the matter power
        spectrum can be described as
        
            P_ mm(k)     =     I^2(k) P_ L(k) + P^1h(k),
        
        with
        
            P^1h(k) =  ∫ dln M_vir n_ln M_virM_vir^2/ρ̅_m ^2| y(k,M_vir) |^2 
                  I(k) = ∫ d ln M_vir n_ln M_virM_vir/ρ̅_m y(k,M_vir) b_ L,
        
        where y(k,M) is the Fourier transform of NFW density profile and normalized as lim_k→0 y(k,M) = 1.
        
        Following this algorithm, we present our final spectrum at z=0 in Fig.<ref>.
        The red and blue curves are Sheth-Tormen and  prescriptions, respectively.
        The same as the mass function results, in order to
        verify our prediction, we compare the spectrum results with N-body simulation data from <cit.>. In details,
        we use the cosmology with (Ω_m=0.24, Ω_Λ=0.76, H_0=73, n_s=0.958,σ_8=0.8).
        The linear power spectrum is output from  Hu-Sawicki f(R) module <cit.>.
        From Fig.<ref>, we can clearly see that the  spectrum results agree with N-body data within a
        few percentage accuracy up to k∼1 h/ Mpc.
        We shall emphasize here that, unlike the halofit philosophy, in this model,  only has two parameters, namely
        r_01, r_02, which denote  the averaged screening and Compton wavelength scales, respectively.
        
        In conclusion, we developed a fast numerical halo model algorithm for modelling non-linear power spectrum for
        the alternative models to ΛCDM. As an example, we show the case of Hu-Sawicki f(R) gravity. The typical
        CPU time with the current parallel Python script[ repository: <https://github.com/hubinitp/CHAM>] (8 threads) is roughly within 10 minutes.  The resulting spectra are in a good agreement with N-body
        data within a few percentage accuracy up to k∼1 h/ Mpc. More importantly, this method is very predictive
        and it does not ask for the calibration from N-body simulation. 
        We believe this method can be widely used in several aspects of data analysis, such as covariance matrix, parameter estimation,  etc.
        
        
        
        We thank Hou-Jun Mo, Zu-Hui Fan, Lucas Lombriser for helpful discussion.
        BH is supported by the Beijing Normal University Grant under the reference No. 312232102,
        Chinese National Youth Thousand Talents Program and the Fundamental Research Funds for
        the Central Universities under the reference No. 310421107. XWL and RGC are supported by  the National Natural Science Foundation of China Grants No.11690022,  No.11435006 and No.11647601, and by the Strategic Priority Research Program of CAS Grant No.XDB23030100 and by the Key Research Program of Frontier Sciences of CAS.
        
        
        apsrev4-1
