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