Error in background.radial_comoving_distance ?

[MickaelRigault]

Hello guys,

I'm sorry but I have a stupid question. I have the feeling that there is an issue of h in the function background.radial_comoving_distance: The doc says it is returning distances in Mpc/h while I think it returns it in Mpc * h.

Here is an example comparing to astropy.cosmology

import jax_cosmo as jcosmo
from astropy.cosmology import Planck15 as astropy_planck15
jc_planck15 = jcosmo.Planck15()

# jax-cosmo
z = jnp.asarray([0.5, 2.5])
dist_mpch = jcosmo.background.radial_comoving_distance(jref_cosmo, a=jcosmo.utils.z2a(z))
# so
dist_mpc_expected = dist_mpch * jc_planck15.h  # Mpc/h * h => Mpc
dist_mpc_devided = dist_mpch / jc_planck15.h  # Mpc/h / h => Mpc/h^2

# astropy
dist_mpc_astropy = astropy_planck15.comoving_distance(z).value # retuns in Mpc

print(f"astropy: {dist_mpc_astropy}")
print(f"expected dist_mpch * jc_planck15.h: {dist_mpc_expected}")
print(f"dist_mpch / jc_planck15.h: {dist_mpc_devided}")

astropy: [1945.56126208 5971.73020615]
expected dist_mpch * jc_planck15.h: [ 893.26105 2744.334  ]
dist_mpch / jc_planck15.h: [1946.6506 5980.625 ]

[MickaelRigault]

See also:

import jax.numpy as jnp
import numpy as np
import jax_cosmo as jcosmo
import matplotlib.pyplot as plt

from astropy.cosmology import Planck15 as astropy_planck15
jc_planck15 = jcosmo.Planck15()

z = jnp.linspace(0.001, 1.5, 100).astype("float32")

fig, ax = plt.subplots()

dist_mpch = jcosmo.background.radial_comoving_distance(jc_planck15, a=jcosmo.utils.z2a(z)) # in Mpc/h (??) or Mpc*h (!!)

ax.plot(z, dist_mpch / jc_planck15.h, label="dist_mpch / jref_cosmo.h", ls=":", lw=5)
ax.plot(z, dist_mpch * jc_planck15.h, label="dist_mpch * jref_cosmo.h", ls="--", )
ax.plot(z, astropy_planck15.comoving_distance(np.asarray(z)).value, label="astropy comoving_distance", ls="-")
ax.legend()

ax.set_ylabel("dist [Mpc]")
ax.set_xlabel("redshift")

[jecampagne]

Hello had you a look at this notebook which is accompanying the Jax-Cosmo paper?

import pyccl as ccl
from jax_cosmo import Cosmology, background

ccl.spline_params.A_SPLINE_NLOG=4*1024
ccl.physical_constants['T_CMB'] = 2.726  # as jax-cosmo
cosmo_ccl= ccl.Cosmology(
    Omega_c=0.3, Omega_b=0.05, h=0.7, sigma8 = 0.8, n_s=0.96, 
    Neff=0,Omega_g=0,
    transfer_function='eisenstein_hu', matter_power_spectrum='halofit')

cosmo_jax = Cosmology(Omega_c=0.3, Omega_b=0.05, h=0.7, sigma8 = 0.8, n_s=0.96,
                      Omega_k=0., w0=-1., wa=0.)
...
# Test array of scale factors
a_min = 0.01
a = np.linspace(a_min, 1.,100)

chi_ccl = ccl.comoving_radial_distance(cosmo_ccl, a)
chi_jax = background.radial_comoving_distance(cosmo_jax, a, log10_amin=np.log10(a_min),steps=5*1024)/cosmo_jax.h

[MickaelRigault]

But background.radial_comoving_distance is supposed to be in Mpc/h ?
(docs says: Radial comoving distance in [Mpc/h] for a given scale factor)

So, to get distances in Mpc you should *multiply by cosmo_jax.h, not divide by it.

In your example (as in mine actually) indeed you find the good chi_jax by dividing which means that your background.radial_comoving_distance returns distances in Mpc * h

So my bug report

[hsimonfroy]

Hello,

Correct me if I am wrong, but if $r$ is in $[\text{Mpc}]$, then $r \times h$ is indeed in $[\text{Mpc}/h]$. This is because $h$ is implicitly in its own unit $[h^{-1}]$, so we have $r [\text{Mpc}] \times h [h^{-1}] = r \times h [\text{Mpc}/h]$. For example, BAO constrains $r_d h [\text{Mpc}/h]$.

Exactly the same way that $3 [\text{km}] \times 1000 [1000^{-1}] = 3000 [\text{km}/1000] = 3000 [\text{m}]$

Then if r_mpcph is in Mpc/h, we indeed should divide by h to get r_mpc = r_mpcph / h in Mpc. So the implementation seems fine.

[MickaelRigault]

You mean that “cosmo_jax.h" is in unit of h^-1 ?
If that is true, that this is a non-sense to me.

Though:

import jax_cosmo as jcosmo
jc_planck15 = jcosmo.Planck15()
jc_planck15.h

0.6774

So h= H0/100, then... h is not in h^-1 unit, right ?

It's very confusing...