Two issues about Halo Model and Matter Power Spectrum

[CaoYe2024]

(1) Lines form 108 to 116 in pyccl/halos/halo_model.py

def _integrate_over_mf(self, array_2):

∫ dM n(M) f(M)

i1 = self._integrator(self._mf * array_2, self._lmass)
return i1 + self._mf0 * array_2[..., 0]

def _integrate_over_mbf(self, array_2):

∫ dM n(M) b(M) f(M)

i1 = self._integrator(self._mf * self._bf * array_2, self._lmass)
return i1 + self._mbf0 * array_2[..., 0] 

self._mf is the halo mass function, self._bf is the halo bias, and array_2 is halo profile. self._integrator = simpson (from scipy.integrate import simpson) i1 term has an extra part which should not be there. Shouldn't these two functions only return “i1”?

(2) In pyccl/halos/pk_2pt.py

p_of_k_a (:class:~pyccl.pk2d.Pk2D): a Pk2D object to be used as the linear matter power spectrum. If None, the power spectrum stored within cosmo will be used.

The linear matter power spectrum used in PyCCL seems to be incorrect. We compare the result at z = 0 with Planck 2018 Linear-theory Matter Power Spectrum, it is very different from the Planck 2018.

[damonge]

(1) The additional term is there to incude the contribution from masses below the minimum mass for which our interpolators work. This is relevant for some calculations. See here

(2) I would be surprised if the matter power spectrum in CCL is incorrect, as it is one of the most heavily validated quantities in the package. If you have a code that reproduces the issue, we could take a look.

[CaoYe2024]

(2) I would be surprised if the matter power spectrum in CCL is incorrect, as it is one of the most heavily validated quantities in the package. If you have a code that reproduces the issue, we could take a look.

Pk.pdf I plotted the matter power spectrum at z=0, given the same cosmological parameters. The blue is calculated by Pyccl and the orange is calculated by CLASS, the orange is the same as Planck. Can you check it?

[damonge]

Hi. My guess is that there is a mismatch of h units. CCL uses Mpc, not Mpc/h. It's hard to know without the code that generated this though.

[CaoYe2024]

Hi. My guess is that there is a mismatch of h units. CCL uses Mpc, not Mpc/h. It's hard to know without the code that generated this though.

I don't think it's the mismatch of h units, I've converted them to the same unit. I wrote a simple code to compare the spectra.

import numpy as np
from classy import Class
import pyccl as ccl
import matplotlib.pyplot as plt

k_min_h_by_Mpc = 1.e-5
k_max_h_by_Mpc = 1.e5
zmin  = 0.
zmax  = 10.
nzmax = 501

cosmo = ccl.CosmologyVanillaLCDM()
Pm_pyccl = cosmo.parse_pk(None)

cosmo_class=Class()
cosmo_class.set({'output': 'mPk','z_max_pk': zmax,'P_k_max_h/Mpc': k_max_h_by_Mpc,'z_pk': '0.0',
    'n_s': cosmo['n_s'],'h': cosmo['h'],'Omega_b': cosmo['Omega_b'],'Omega_cdm': cosmo['Omega_c'], 
    'sigma8':cosmo['sigma8']})
cosmo_class.compute()
z = np.linspace(0, zmax, num=nzmax)
# Get power spectrum P(k=l/chi(z),z) from cosmological module
kmin_in_inv_Mpc = k_min_h_by_Mpc * cosmo_class.h()
kmax_in_inv_Mpc = k_max_h_by_Mpc * cosmo_class.h()
nk=512
k = np.logspace(-4, 4,nk)
pk = np.zeros(( nzmax , nk))
for index_z in range(nzmax):
    for index_k in range(nk):
        pk[index_z,index_k] = cosmo_class.pk(k[index_k], z[index_z])

a = 1./(1+z[::-1])
Pm_classy = ccl.pk2d.Pk2D(lk_arr=np.log(k),a_arr=a, pk_arr=pk[::-1,:], is_logp=False)

plt.plot(k/cosmo['h'],Pm_pyccl(k,1)* cosmo['h']**3)
plt.plot(k/cosmo['h'],Pm_classy(k,1)* cosmo['h']**3)
plt.xlim(0.002,6)
plt.ylim(1,200000)
plt.ylabel('$P_{m}(k)\ [(h^{-1}Mpc)^{3}]$')
plt.xlabel('$k\ [hMpc^{-1}]$')
plt.semilogx()
plt.semilogy()

[damonge]

Hmmm, I'm not sure what to tell you. The code below (which is just a very mildly modified version of yours) returns the following result:

import numpy as np
from classy import Class
import pyccl as ccl
import matplotlib.pyplot as plt

zmin  = 0.
zmax  = 10.
nzmax = 501

cosmo = ccl.CosmologyVanillaLCDM()
Pm_pyccl = cosmo.parse_pk(None)

cosmo_class=Class()
cosmo_class.set({'output': 'mPk','z_max_pk': zmax,'P_k_max_h/Mpc': 1E4, 'z_pk': '0.0',
    'n_s': cosmo['n_s'],'h': cosmo['h'],'Omega_b': cosmo['Omega_b'],'Omega_cdm': cosmo['Omega_c'], 
    'sigma8':cosmo['sigma8']})
cosmo_class.compute()

z = np.linspace(0, zmax, num=nzmax)
nk=512
k = np.logspace(-3, 3, nk)
pk = np.zeros(( nzmax , nk))
for index_z in range(nzmax):
    for index_k in range(nk):
        pk[index_z,index_k] = cosmo_class.pk(k[index_k], z[index_z])
a = 1./(1+z[::-1])
Pm_classy = ccl.pk2d.Pk2D(lk_arr=np.log(k),a_arr=a, pk_arr=pk[::-1,:], is_logp=False)

plt.plot(k/cosmo['h'],Pm_pyccl(k,1)* cosmo['h']**3, 'r-', label='CCL')
plt.plot(k/cosmo['h'],Pm_classy(k,1)* cosmo['h']**3, 'k--', label='Classy')
plt.xlim(0.002,6)
plt.ylim(1,200000)
plt.semilogx()
plt.semilogy()
plt.legend()

[CaoYe2024]

Can you tell me which version of PyCCL you are currently using?

[damonge]

The latest public one (3.something). I don't think anything has changed in terms of P(k) for a while now, so older versions should return the same answer, I think.

[CaoYe2024]

The latest public one (3.something). I don't think anything has changed in terms of P(k) for a while now, so older versions should return the same answer, I think.

I updated pyccl to version 3.0.2 using conda, then I got the same result as you. I reinstalled version 2.8.0, but the previous error did not reoccur. I think it is possible that some file in my previous pyccl was tampered with. Of course, it also updated some other packages while updating pyccl with conda, possible that the versions of previous package were too old.I'm not sure what the problem is. Thank you very much！