f_z inconsistencies in camb interface

[jessmuir]

I just wanted to flag some inconsistencies in the growth_parameters section of the camb_interface.py.

Currently, the interface gets sigma_8(z) and fsigma_8(z) from CAMB output, but then computes f(z) by doing a spline derivative of the linear matter power spectrum at the smallest available k. As a result, fsigma8(z)/sigma8(z) != f(z).

I think there are two things going on here.

First is that CAMB computes fsigma8 by evaluating
sigma8^2_vd/sigma8, where sigma8^2_vd is a sigma8-like quantity computed by integrating over the linear cross-power spectrum between delta and velocity. If you then divide sigma8(z) out from that to isolate f(z), you may get a slightly different results than just picking k=kmin and doing a derivative of P_dd(kmin,z), since the CAMB fsigma8 quantity will be slightly sensitive to scale-dependent effects (e.g. from neutrino mass).

I imagine there could be reasons for someone to want the spline derivative f(z) rather than the sigma_vd-based one, but it might be a good idea to flag in the documentation that these are not the same quantity and make sure that modules that need growth rate e.g. for RSD calculations (like the projection module) use the CAMB-based quantities.

Second, I think there might be a bug in how the camb interface computes or stores f(z). I'm not sure what's happening, but when I compare the f_z output by a test sampler to spline derivatives I've done of the linear P(k,z) output, there are some pretty dramatic mismatches. In the plot below I'm comparing outputs from a test sampler run. In it, the line labeled "cosmosis interface" is what is stored in the f_z.txt file. The black "CAMB" line shows what happens if I divide the contents of fsigma_8.txt by sigma_8.txt, and the "manual" lines show my calculations based on the tabulated linear matter power spectrum.

N.b. To do the manual calculations I followed the camb interface's compute_growth_rates pretty closely so I'm not sure where this difference is coming from:

plin=np.loadtxt("p_k.txt")
k = np.loadtxt("k_h.txt")
z_pow = np.loadtxt("z.txt")
pkmin = plin[:,0]
kmin = k[0]
a_pow = 1/(1+z_pow)
loga = np.log(a_pow[::-1])
Dkmin = np.sqrt(pkmin/pkmin[0]).squeeze()[::-1]

logD_spline_min = InterpolatedUnivariateSpline(loga, np.log(Dkmin) )
f_spline_min = logD_spline_min.derivative()
fmin = f_spline_min(loga)[::-1].squeeze()

While I'm generating plots, this is just to show that the manual spline interpolation and the CAMB fsigma8/sigma8 calculation are (very) slightly different. This difference becomes larger if for larger neutrino masses.

[joezuntz]

Hi Jessie.
Thanks for looking so carefully at this.

For your first plot - can I check if you're using the latest version of the standard lib? We merged some fixes here: https://github.com/joezuntz/cosmosis-standard-library/pull/45/files

Fot the different definitions of f(z) - thanks so much, that's really useful. Maybe I should make two outputs, one using the camb definition (which seems to come from Planck for "general models") and one using the traditional version?

[jessmuir]

Hi Joe,

Just checked, and I was in fact using an older version of cosmosis-standard-library. Oops! Pulling the update and rerunning makes the cosmosis f_z match the CAMB fsig8/sig8 calculation, which appears to be what is implemented in the merged fixes you linked.

For the different definitions, maybe having the options for both outputs could be useful, either outputting the spline derivative version to some other tabulated quantity (f_z_atkmin?); or having the option to switch from one to the other. That's could probably be a low priority task, though: My guess is that in most cases the CAMB output is what'll be most relevant, and should be the default, since using the linear velocity power spectrum means it'll be capturing the effect relevant for RSD.

[jessmuir]

Related to this: Inside the projection module cosmosis-standard-library/structure/projection/project_2d.py, f(z) gets computed again separately, using the spline approach. The effect of this is small, but ideally it would be preferable to switch the projection module over to reading f(z) and D(z) from the datablock rather than recomputing them, to help make sure everything stays consistent.

I think the only place that would need edits occurs inside the function set_nonlimber_splines