[BUG] growth_rate, currently implemented following Hamilton (2001) eq. 4, is very far from the CAMB growth rate

[MinhMPA]

Describe the bug
Currently, growth_rate() in the growth model GrowthFactor implements the linear growth rate $f(a)$ following eq. 4 of Hamilton (2001). However, this particular eq. does not yield an accurate growth rate when compared to $d\ln D/d\ln a$ from CAMB and the growth model CambGrowth.

My suggestion, if I may, and I can make a PR for this, is to either a/ use Hamilton (2001) eq. 5 instead, or b/ simply $f(a)=\Omega_m(a)^\gamma$ where $\gamma=0.55$.

To Reproduce
Code snippet reproduce the behavior:

import matplotlib.pyplot as plt
%matplotlib inline

from hmf import MassFunction,Transfer

transfer_CambGrowth=Transfer(z=0.0, lnk_min=-18.420680743952367, lnk_max=9.903487552536127, dlnk=0.05, transfer_model='CAMB', transfer_params=None, takahashi=True, growth_model='CambGrowth', growth_params=None, use_splined_growth=True)
transfer_GrowthFactor=Transfer(z=0.0, lnk_min=-18.420680743952367, lnk_max=9.903487552536127, dlnk=0.05, transfer_model='CAMB', transfer_params=None, takahashi=True, growth_model='GrowthFactor', growth_params=None, use_splined_growth=True)

z=np.linspace(0.,100.,256)
a=1./(1.+z)
D_CambGrowth=transfer_CambGrowth.growth.growth_factor(z)
f_CambGrowth=np.gradient(np.log(D_CambGrowth),np.log(a))
f_GrowthFactor=transfer_GrowthFactor.growth.growth_rate(z)

fig, ax = plt.subplots(2, sharex=True, sharey=False, gridspec_kw={'hspace': 0., 'wspace':0., 'height_ratios':[8, 2]}, figsize=(6.4, 4.8))
ax[0].semilogx(a,f_GrowthFactor,c='b',ls='-',label=r'GrowthFactor, Hamilton Eq. (4)')
ax[0].semilogx(a,f_CambGrowth,c='r',ls='--',label=r'CambGrowth')
ax[0].legend(frameon=False)
ax[0].set_ylabel(r'$f(a)$')
#ax[0].set_ylim(0.65,1.)
ax[1].semilogx(a,f_GrowthFactor/f_CambGrowth)
ax[1].set_ylabel(r'$f_{\mathrm{GF}}/f_{\mathrm{CG}}$')
ax[1].set_xlabel(r'$a$')

Expected behavior
This is the comparison between either Hamilton eq. 4 or eq. 5, against the CAMB numerical derivative $d\ln D/d\lna$.

Additional context
This is related to, but not the same issue, #196.

[steven-murray]

Hmmm, this is interesting. In principle, Eq (4) should be more accurate than Eq (5), right? Is your Eq (4) plot including the fix in #197?

[MinhMPA]

Yes. The fix in #197 is included. If I were not to include it, the growth rate from GrowthFactor would have been mis-normalized, about 2 orders of magnitude higher than the growth rate from CambGrowth.

I am actually not familiar with Eq(4) in Hamilton 2000. So I went back to check the paper and Hamilton actually cited Lahav et al. 1991. Checking the latter, I think the proper expression for growth rate, should be their Eq(9).

Here is what I get if I implement Eq(9) in Lahav et al. 1991 instead:

[MinhMPA]

@steven-murray Actually, we can also follow the definition $f\equiv\frac{d\ln D}{d\ln a}=\frac{dD/da}{D}a$*, like so

def growth_rate(self, z):
  a = np.exp(self._lna)
  Da = self.growth_factor(-1.+(1./a))
  Dadot_spline = _spline(a, Da).derivative()
  az = 1./(1.+z)
  return Dadot_spline(az)/self.growth_factor(z)*az

Then I would still get a growth rate $f$ much closer to CAMB, although some numerical errors might accumulate through the antiderivative() inside growth_factor and the derivative() here.

In summary, I'm really skeptical about Hamilton Eq. (4). Maybe I'm being dense but I just can't derive Hamilton Eq. (4) straight from the definition of $f$.

*One can also write this as $f=1+\frac{dg/da}{g}$.

[MinhMPA]

Another technical detail, independently from all this, that makes me wonder is, whether picking the reference scale k=1.0---when trying to extract the growth factor from CAMB, as done currently in hmf---really makes sense. If we instead used a much lower reference k, naively I'd expect we get a cleaner growth factor since those modes would not be strongly affected by the linear transfer function $T(k)$.

[steven-murray]

Thanks @MinhMPA -- I think you're right about the CAMB reference scale. I'm not sure why I originally set it at k=1.0, but lower sounds like a better idea (as a default, anyway). I'm still a little nervous that implementing f in four different ways is giving us four different answers. This would be fine if I understood why the differences existed, and so we could pick a best default. But I'll have to think about this more systematically.