camb计算物质功率谱

今天利用camb计算宇宙的dark matter物质功率谱遇到一个很费解的问题，函数get_matter_power_interpolator()两个参数官方文档是这样写的：

hubble_units – if true, output power spectrum in (({\rm Mpc}/h)^{3}) units, otherwise ({\rm Mpc}^{3})
k_hunit – if true, matter power is a function of k/h, if false, just k (both ({\rm Mpc}^{-1}) units)

这个参数hubble_units很好理解，如果为True，那么输出的功率谱是带h单位的功率谱，也就是$({\rm Mpc}/h)^3$；
k_hunit很令人费解，为什么不管True还是False，单位都是${\rm Mpc}^{-1}$，如果为True则k/h；

具体理解

为了理解这点可以从带哈勃常数h单位来理解：
假设$h = 0.7$，那么：

长度：$1 {\rm Mpc} = 0.7 {\rm Mpc}/h$

这里的0.7是怎么计算的呢？可以将变换后的数值设为x，那么：

$1 {\rm Mpc} = x {\rm Mpc}/h$

消除相同项目，因此$x = 1\cdot h = 0.7$。总之，单位引入h的意义是不管h怎么变化，换算成实际长度都会不变。

同理对于波数$k$：

波数：$k\ h/{\rm Mpc} = x\ 1/{\rm Mpc}$

因此：$x = k h$

如果左右换过来转换：

波数：$k\ 1/{\rm Mpc} = x\ h/{\rm Mpc}$

因此：$x = k/h$，这也就符合官方文档$k/h$的描述。可以推理官方文档的意思是，如果k_hunit=True，那么计算的出来的函数(比如PK.P(rz_center, k))，中有关$k$的项包含单位，官方文档“不管True还是False，单位都是$Mpc^{-1}$”意思是内部不管什么都是不带h单位计算，使用者只需要关注外部带不带h单位即可。

例子：

import camb
import numpy as np
import matplotlib.pyplot as plt
from astropy.cosmology import Planck18
from camb import model, initialpower
# 创建CAMB参数实例
from matplotlib.ticker import FuncFormatter

from rsd_ps.catalog import get_h5

h = Planck18.h
rz_center = 0.3

k = np.linspace(1.e-2, 1, 10) # 假设单位 h/Mpc
print(k.min(), k.max())

# 设置参数
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.66, ombh2=0.022, omch2=0.12, mnu=0.06,meffsterile=0.0,nnu=3.044, omk=0, tau=0.05)
pars.set_dark_energy(w=-1.0,wa=0)
pars.InitPower.set_params(As=2e-9, ns=0.965, r=0)
pars.set_for_lmax(lmax=4000, lens_potential_accuracy=4)

pars.Transfer.kmin = 0

pars.Transfer.kmax = 1

# 设置精度参数
pars.AccuracyBoost = 2
pars.lSampleBoost = 1
pars.lAccuracyBoost = 1

# 设置其他参数
pars.WantTensors = False
pars.DoLensing = True
pars.set_matter_power(redshifts=[rz_center])
# pars.NonLinear = model.NonLinear_both
pars.NonLinear = 3


# 计算结果
results = camb.get_results(pars)

# 物质功率谱
# results.calc_power_spectra(pars)
PK = results.get_matter_power_interpolator(k_hunit=True, hubble_units=True)  # 功率谱单位(Mpc/h)^3
ps_matter = PK.P(rz_center, k)

print(ps_matter)

PK_onh = results.get_matter_power_interpolator(k_hunit=False, hubble_units=True)  # 功率谱单位(Mpc/h)^3
ps_matter_noh = PK_onh.P(rz_center, k*h)

print(ps_matter_noh)

输出：

[15368.24577508  3142.99422918  1245.50941985   730.35516333
   510.35717433   394.2712478    323.50403882   275.64542495
   240.71770316   213.7227539 ]
[15368.24577508  3142.99422918  1245.50941985   730.35516333
   510.35717433   394.2712478    323.50403882   275.64542495
   240.71770316   213.7227539 ]

参考

https://www.astro.ljmu.ac.uk/~ikb/research/h-units.html