Repository: skypyproject/skypy
Pull request #445 (2021-03-03 16:10:06)

BUG: Fixed wavenumber definition and amplitude scaling in Eisenstein & Hu model
## Description

The Eisenstein & Hu approximation for the linear matter power spectrum now assumes the input wavenumber vector to be in units of Mpc instead of Mpc/h (See  [ADR-2](https://github.com/skypyproject/skypy/blob/fad2b7624c87efdbbddf23e249480006ad8c224c/docs/arch/adr-02.md))

The overall amplitude was also missing factors of h. This has now been fixed such that multiplying `eisenstein_hu` by the square of `growth_function_carroll` gives a power spectrum that closely matches CAMB and CLASS.

This addresses issue #442

```python
from astropy.cosmology import Planck18
from matplotlib import pyplot as plt
import numpy as np
from skypy.power_spectrum import camb, eisenstein_hu, growth_function_carroll

k = np.logspace(-3, -1, 50)
z, A_s, n_s = 0, 2.2E-9, 0.97

p_camb = camb(k, z, Planck18, A_s, n_s)
p_eh = eisenstein_hu(k, A_s, n_s, Planck18)
d = growth_function_carroll(z, Planck18)

plt.plot(k, p_camb, label='CAMB')
plt.plot(k, p_eh * d * d, label='Eisenstein-Hu')
plt.loglog()
plt.legend()
plt.show()
```
Returns a power spectrum which has the correct wavenumbers and amplitude
![download](https://user-images.githubusercontent.com/34517350/109840622-32fb1900-7c40-11eb-9c41-5c94207f7b0d.png)
 
#

## Checklist
- [x] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/main/CONTRIBUTING.rst)
- [x] Write unit tests
- [x] Write documentation strings
- [x] Assign someone from your working team to review this pull request
- [x] Assign someone from the infrastructure team to review this pull request



Patch:
diff --git a/skypy/power_spectrum/_eisenstein_hu.py b/skypy/power_spectrum/_eisenstein_hu.py
index a1c8cea..7796c5c 100644
--- a/skypy/power_spectrum/_eisenstein_hu.py
+++ b/skypy/power_spectrum/_eisenstein_hu.py
@@ -6,6 +6,7 @@ formula for the matter power spectrum.
 
 from astropy.utils import isiterable
 import numpy as np
+from astropy import constants
 
 
 __all__ = [
@@ -54,8 +55,8 @@ def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02):
     >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
     >>> A_s, n_s = 2.1982e-09, 0.969453
     >>> transfer_with_wiggles(wavenumber, A_s, n_s, Planck15, kwmap=0.02)
-    array([9.92144790e-01, 7.78548704e-01, 1.29998169e-01, 4.63863054e-03,
-       8.87918075e-05])
+    array([9.84445487e-01, 6.77428507e-01, 8.21122028e-02, 2.44250638e-03,
+       4.41412308e-05])
 
     References
     ----------
@@ -84,7 +85,6 @@ def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02):
         raise ValueError("Tcmb0 for input cosmology must be non-zero if"
                          "wiggles = True")
 
-    ak = wavenumber * h0
     om0h2 = om0 * h0**2
     ob0h2 = ob0 * h0**2
     f_baryon = ob0 / om0
@@ -126,8 +126,8 @@ def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02):
     beta_b = 0.5 + f_baryon + (3 - 2 * f_baryon) * np.sqrt((17.2 * om0h2)**2
                                                            + 1.0)
 
-    q = ak / (13.41 * k_eq)
-    ks = ak * sound_horizon
+    q = wavenumber / (13.41 * k_eq)
+    ks = wavenumber * sound_horizon
 
     T_c_ln_beta = np.log(np.e + 1.8 * beta_c * q)
     T_c_ln_nobeta = np.log(np.e + 1.8 * q)
@@ -143,11 +143,11 @@ def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02):
         (1 - T_c_f) * f(T_c_ln_beta, T_c_C_alpha)
 
     s_tilde = sound_horizon * (1 + (beta_node / ks)**3)**(-1 / 3)
-    ks_tilde = ak * s_tilde
+    ks_tilde = wavenumber * s_tilde
 
     T_b_T0 = f(T_c_ln_nobeta, T_c_C_noalpha)
     T_b_1 = T_b_T0 / (1 + (ks / 5.2)**2)
-    T_b_2 = alpha_b / (1 + (beta_b / ks)**3) * np.exp(-(ak / k_silk)**1.4)
+    T_b_2 = alpha_b / (1 + (beta_b / ks)**3) * np.exp(-(wavenumber/k_silk)**1.4)
     T_b = np.sinc(ks_tilde / np.pi) * (T_b_1 + T_b_2)
 
     transfer = f_baryon * T_b + (1 - f_baryon) * T_c
@@ -190,8 +190,8 @@ def transfer_no_wiggles(wavenumber, A_s, n_s, cosmology):
     >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
     >>> A_s, n_s = 2.1982e-09, 0.969453
     >>> transfer_no_wiggles(wavenumber, A_s, n_s, Planck15)
-    array([9.91959695e-01, 7.84518347e-01, 1.32327555e-01, 4.60773671e-03,
-       8.78447096e-05])
+    array([9.84321214e-01, 6.83446969e-01, 8.15258290e-02, 2.42055117e-03,
+       4.36725897e-05])
 
     References
     ----------
@@ -210,7 +210,6 @@ def transfer_no_wiggles(wavenumber, A_s, n_s, cosmology):
     ob0 = cosmology.Ob0
     h0 = cosmology.H0.value / 100
     Tcmb0 = cosmology.Tcmb0.value
-    ak = wavenumber * h0
     om0h2 = om0 * h0**2
     f_baryon = ob0 / om0
 
@@ -218,8 +217,8 @@ def transfer_no_wiggles(wavenumber, A_s, n_s, cosmology):
         np.log(22.3 * om0h2) * f_baryon**2
     sound = 44.5 * np.log(9.83 / om0h2) / \
         np.sqrt(1 + 10 * (f_baryon * om0h2)**(0.75))
-    shape = om0h2 * (alpha + (1 - alpha) / (1 + (0.43 * ak * sound)**4))
-    aq = ak * (Tcmb0 / 2.7)**2 / shape
+    shape = om0h2 * (alpha + (1 - alpha) / (1 + (0.43 * wavenumber * sound)**4))
+    aq = wavenumber * (Tcmb0 / 2.7)**2 / shape
     transfer = np.log(2 * np.e + 1.8 * aq) / \
         (np.log(2 * np.e + 1.8 * aq) +
          (14.2 + 731 / (1 + 62.5 * aq)) * aq * aq)
@@ -272,8 +271,8 @@ def eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap=0.02, wiggle=True):
     >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
     >>> A_s, n_s = 2.1982e-09, 0.969453
     >>> eisenstein_hu(wavenumber, A_s, n_s, Planck15, kwmap=0.02, wiggle=True)
-    array([6.47460158e+03, 3.71610099e+04, 9.65702614e+03, 1.14604456e+02,
-       3.91399918e-01])
+    array([2.99126326e+04, 1.32023496e+05, 1.80797616e+04, 1.49108261e+02,
+           4.53912529e-01])
 
     References
     ----------
@@ -282,7 +281,8 @@ def eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap=0.02, wiggle=True):
     .. [3] Komatsu et al., ApJS, 180, 330 (2009)
     """
     om0 = cosmology.Om0
-    h0 = cosmology.H0.value / 100
+    H0_c = (cosmology.H0 / constants.c).to_value('Mpc-1')
+
     if wiggle:
         transfer = transfer_with_wiggles(wavenumber, A_s, n_s, cosmology,
                                          kwmap)
@@ -290,8 +290,8 @@ def eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap=0.02, wiggle=True):
         transfer = transfer_no_wiggles(wavenumber, A_s, n_s, cosmology)
 
     # Eq [74] in [3]
-    power_spectrum = A_s * (2 * wavenumber**2 * 2998**2 / 5 / om0)**2 * \
-        transfer**2 * (wavenumber * h0 / kwmap)**(n_s - 1) * 2 * \
-        np.pi**2 / wavenumber**3
+    power_spectrum = A_s * (2 * (wavenumber / H0_c)**2 / 5 / om0)**2 * \
+        transfer**2 * (wavenumber / kwmap)**(n_s - 1) * 2 * \
+        np.pi**2 / (wavenumber)**3
 
     return power_spectrum


Test patch:
diff --git a/skypy/power_spectrum/tests/test_eisenstein_hu.py b/skypy/power_spectrum/tests/test_eisenstein_hu.py
index fa670c8..e080659 100644
--- a/skypy/power_spectrum/tests/test_eisenstein_hu.py
+++ b/skypy/power_spectrum/tests/test_eisenstein_hu.py
@@ -37,13 +37,13 @@ def test_eisenstein_hu():
                                       wiggle=True)
     pk_eisensteinhu_nw = eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap,
                                        wiggle=False)
-    pk_cosmosis_w = np.array([6.47460158e+03, 3.71610099e+04, 9.65702614e+03,
-                              1.14604456e+02, 3.91399918e-01])
-    pk_cosmosis_nw = np.array([6.47218600e+03, 3.77330704e+04, 1.00062077e+04,
-                              1.13082980e+02, 3.83094714e-01])
+    pk_pre_w = np.array([2.99126326e+04, 1.32023496e+05, 1.80797616e+04,
+                         1.49108261e+02, 4.53912529e-01])
+    pk_pre_nw = np.array([2.99050810e+04, 1.34379783e+05, 1.78224637e+04,
+                         1.46439700e+02, 4.44325443e-01])
 
-    assert np.allclose(pk_eisensteinhu_w, pk_cosmosis_w)
-    assert np.allclose(pk_eisensteinhu_nw, pk_cosmosis_nw)
+    assert np.allclose(pk_eisensteinhu_w, pk_pre_w)
+    assert np.allclose(pk_eisensteinhu_nw, pk_pre_nw)
 
     # Test for failure when wavenumber <= 0
     negative_wavenumber_scalar = 0