Checking (and fixing) transfer functions

[jmsull]

Starting a closer comparison of Bolt perturbations with CLASS transfer functions at the same cosmology at high (z~3000) and low (z=0) redshifts.

[jmsull]

Already a quick check shows qualitative agreement on many of the species, but big problems with a few of them.

Here are both transfer functions on the same plot (split up a bit to make them less painful to look at):

High redshift:

(initially) relativistic species monopoles (dipoles show similar agreement):

(initially) nonrelativistic species:

Low redshift:

(initially) relativistic species monopoles (dipoles show similar agreement):

(initially) nonrelativistic species:

These are just to get an idea of where the biggest problems are at the moment. I will list a few:

• Something is clearly wrong in the massive neutrino normalization in a redshift-dependent way (as they become nonrelativistic), I will look into this.
• The shape for the photon monopole is not at all right - not familiar enough with it atm to know what exactly the problem is? Maybe more ell_max?
• Low-k normalization is generally slightly off at high redshift.

There are a couple of normalization factors I threw in there that should be checked more carefully (but are listed in the legend). k has been scaled to k/(b.H0*c/100) to be in h/Mpc.

Further smaller quantitative disagreements abound but it is unclear (to me) how much is resolution-dependent.

---

Here is the code I used to generate the CLASS transfer functions (just to have a record if a mistake,I will put the saved CLASS transfer files in a future commit) - the Bolt code is in scripts/plot_perts_k.jl (after merging recfast so YHe agrees):

bolt_names = ['h','Omr','Omb','Omm','ns','Y_p']
bolt_cosmo = [0.7,5.042e-5,0.046,0.224,1.0,0.24]
cosmo_dict = dict(zip(bolt_names,bolt_cosmo))
h = cosmo_dict['h']
N_ncdm=1
m_ncdm=0.06 
N_ur=3.046-N_ncdm
cosmo_dict['Omcdm'] = cosmo_dict['Omm']#we should change the name for Omm
ceng = Class()
ceng.set({'h':h,
          'omega_b':cosmo_dict['Omb']*h**2,
          'omega_cdm':cosmo_dict['Omcdm']*h**2,
          'n_s':cosmo_dict['ns'],
          'gauge':'newtonian',
          'N_ur': N_ur,
          'N_ncdm': N_ncdm,
          'm_ncdm': m_ncdm,
          'YHe': cosmo_dict['Y_p']
          })
ceng.set({'output': 'mTk, vTk','P_k_max_1/Mpc':10.0,'z_max_pk':30000.})
ceng.compute()
transfers = ceng.get_transfer(z)
ceng.struct_cleanup()
ceng.empty()
np.savetxt('./class_tf_xm8.dat',np.array([tf for tf in transfers.values()]).T)

[jmsull]

Forgot to update the ionization history according to new recfast format (whoops) - replaced the plots in the above comment, but discussion is basically the same.

[xzackli]

Wow, that photon monopole is very alarming. I'll work on figuring out what's going on there.

[jmsull]

To better diagnose the problem I added a script to plot perturbations at a single k as a function of time (scripts/plot_perts_x.jl) and compare to CLASS (code at the bottom).

Here is a sample of the evolution for modes for massive neturinos at k = 0.3 (unlabeled) and 0.03 h/Mpc (sufficiently subhorizon to be interesting):

Clearly something bad is happening with the plotted massive neutrino perturbations near x=-5 (z~150), which is around when they become non-relativistic. The evolution seems to have roughly the correct shape - the massive neutrinos are certainly feeling their mass and the perturbation grows as they become non-relativistic.

I can recover something pretty close to the CLASS evolution (almost as close as for the other species) if I just set the scale factor used in the normalization of the neutrino perturbations (post-proceessing step, not part of the hierarchy evolution) to be small enough - basically removing the impact of the mass (in epsilon) in the q integral to get the integrated monopole. These are the "relnorm" lines on the above plot.

It may be that since we chose the normalization at the initial time to be wrt the massless neutrinos, this is a choice that should be kept fixed as time evolves to be consistent with CLASS? I will dig around in class_public and see if I can find how they set their normalizations. Then it should be clearer whether this is a physics problem or convention problem.

I have added the CLASS transfer functions and perturbations as a function of x in test/data/ - let me know if you want them to go somewhere else.

Also, for reference here are the rest of the perturbations using the same format as above (it is easier to see the photon monopole problem here):

k=0.03 h/Mpc:

(initially) relativistic species monopoles:

(initially) non-relativistic species monopoles (plus Phi):

k=0.3 h/Mpc:

(initially) relativistic species monopoles:

(initially) non-relativistic species monopoles (plus Phi):

Code for generating CLASS perturbations:

#Copied from class_public notebook - many_times - dropping the things needed to make their beautiful figure for brevity
#i.e. damping scale etc

z_max_pk = 46000 #30000.#      # highest redshift involved 
k_per_decade = 400     # number of k values, controls final resolution
k_min_tau0 = 40.       # this value controls the minimum k value in the figure (it is k_min * tau0)
P_k_max_inv_Mpc =1.0   # this value is directly the maximum k value in the figure in Mpc
tau_num_early = 2000   # number of conformal time values before recombination, controls final resolution
tau_num_late = 200     # number of conformal time values after recombination, controls final resolution
tau_ini = 10.          # first value of conformal time in Mpc

# Cosmological parameters and other CLASS parameters
bolt_names = ['h','Omr','Omb','Omm','ns','Y_p','Neff']
bolt_cosmo = [0.7,5.042e-5,0.046,0.224+0.046,1.0,0.24,3.046] 
cosmo_dict = dict(zip(bolt_names,bolt_cosmo))
cosmo_dict['Omcdm'] = cosmo_dict['Omm']-cosmo_dict['Omb']
common_settings = {
                   'output':'mTk, vTk',
                   'h':cosmo_dict['h'],
                   'omega_b':cosmo_dict['Omb']*cosmo_dict['h']**2,
                   'omega_cdm':cosmo_dict['Omcdm']*cosmo_dict['h']**2,
                   'n_s':cosmo_dict['ns'],
                   'gauge':'newtonian',
                   'N_ur': cosmo_dict['Neff']-1,
                   'N_ncdm': 1,
                   'm_ncdm': 0.06,
                   'YHe': cosmo_dict['Y_p'],
                   'P_k_max_1/Mpc':10.0,
                   'z_max_pk':z_max_pk,
                   'recfast_Nz0':z_max_pk*2.,
                   'recfast_z_initial':z_max_pk+1.,
                   'k_per_decade_for_pk':k_per_decade,
                   'k_per_decade_for_bao':k_per_decade,
                   'perturb_sampling_stepsize':'0.05',
                  }

# call CLASS 
M = Class()
M.set(common_settings)
M.compute()

# define conformal time sampling array
times = M.get_current_derived_parameters(['tau_rec','conformal_age'])
tau_rec=times['tau_rec']
tau_0 = times['conformal_age']
tau1 = np.logspace(np.log10(tau_ini),np.log10(tau_rec),tau_num_early)
tau2 = np.logspace(np.log10(tau_rec),np.log10(tau_0),tau_num_late)[1:]
tau2[-1] *= 0.999 # this tiny shift avoids interpolation errors
tau = np.concatenate((tau1,tau2))
tau_num = len(tau)

background = M.get_background() # load background table

background_tau = background['conf. time [Mpc]'] # read conformal times in background table
background_z = background['z'] # read redshift
background_z_at_tau = interp1d(background_tau,background_z)

# check and inform user whether intiial arbitrary choice of z_max_pk was OK
max_z_needed = background_z_at_tau(tau[0])
if max_z_needed > z_max_pk:
    print( 'you must increase the value of z_max_pk to at least ',max_z_needed)
    () + 1  # this strange line is just a trick to stop the script execution there
else:
    print( 'in a next run with the same values of tau, you may decrease z_max_pk from ',z_max_pk,' to ',max_z_needed)

pkeys = ['k (h/Mpc)', 'd_g', 'd_b', 'd_cdm', 'd_ur', 'd_ncdm[0]', 'd_tot', 'phi', 'psi', 't_g', 't_b', 't_cdm', 't_ur', 't_ncdm[0]', 't_tot']

zz = np.zeros_like(tau) 
k_num = len( M.get_transfer(zz[0])['k (h/Mpc)']) #get number of k pts
#get  all results
perts = np.zeros((len(pkeys),tau_num,k_num))
for i in range(tau_num):
    zz[i] = background_z_at_tau(tau[i])
    one_time = M.get_transfer(zz[i]) # transfer functions at each time 
    for s in range(len(pkeys)):
        perts[s,i,:] = one_time[pkeys[s]][:] 

#get the slice first for a particular k
#this is lazy but just computing the full tau-k grid then slicing, also reversing
def get_kslice(ktarget, eps=0.001):
    k = perts[0,0,:]
    itarget = np.where((k>ktarget-eps) & (k<ktarget+eps))[0][0]
    return perts[:,:,itarget][:,::-1]  

zz=zz[::-1] #flip things around to go from early to late times
xx = np.log(1/(1+zz))

#save them
sliced_pert_1 = get_kslice(0.3)
sliced_pert_2 = get_kslice(0.03)

np.savetxt('./class_px_kp3.dat',np.vstack([xx,sliced_pert_1]))
np.savetxt('./class_px_kp03.dat',np.vstack([xx,sliced_pert_2]))

[jmsull]

After a long journey, I figured out the mistake I made with the neutrino normalization - which in retrospect is a little silly. I thought I was matching some convention by dividing the massive neutrino perturbation by the massless neutrino background energy density since this gives the right result in perturbations at early times and in the initial conditions, as well as in the transfer output for massive neutrinos when compared to CLASS.

However, just sticking with the exact version of the MB hierarchy is what I should have done, because at late times the integrated distribution function changes to reflect the nonrelativistic transition of the neutrinos, which is what was broken using the earlier ad hoc normalization.

The Einstein equations then receive contributions that look different from the other species, but come from a straightforward manipulation of the MB equations (updated this in the notes, I think correctly). There is no need for extra normalization factors.

The neutrinos now make the transition to being nonrelativistic in quite nice agreement with CLASS (few percent).

When plotting the integrated neutrino perturbations (not during solving), the perturbations do now have to be divided by the normalization (the background density) - I introduced a convenience spline for the background density we can take out later if necessary. There is something a little weird going on with this though for the different background vs perturbation integrals, I'll open an issue for it.

On a related note, I added new CLASS perturbations (in x) and transfer functions (in k) to test/data for the case where no fluid approximation is used (which takes a very long time in CLASS). The plot shows the CLASS fluid approximation as well in orange. CLASS uses a fluid approximation by default and I was seeing a factor of 2 error in the perturbations, but this goes away with an apples to apples comparison.

---

On the photon multipole issue - this appears related to \ell_max.

When I bump up \ell\gamma to 1000 I get the following plot for k=0.03 h/Mpc. Solid is CLASS, green dashed is \ell\gamma=1000, orange dashed \ell_\gamma 10 (similar to before). Compared with before we resolve more oscillations in time before getting the divergent behavior at late x, which is an improvement from the earlier plot above. Using higher \ell_nu we resolve the later time massless neutrino oscillations better as well. Given that MB says they used ell max of 2000 for massless particles I am not sure this is surprising?

---

Finally, using Rodas5 seems to be okay now (no errors, though still slower than KenCarp4). I'll close issue #40

For all the perturbations, there are still several percent deviations even at very early times, I guess this is the next thing to look into. I don't believe it is due to the neutrinos (at least I am more confident it is not an error in implementation). Part of it may be that the k mode values for CLASS and Bolt are off by a few percent. Also, adding multipoles seems to have some effect on the agreement (e.g. \Phi). But I guess it is more than either of those (maybe due to early time evolution, I think CLASS mentions they evolve \Psi for numerical stability instead of using the constraint equation, can check this by computing CLASS perturbations at higher z).

[xzackli]

This is really impressive work. The resolution to the Rodas solver mystery is very satisfying, although it is still surprising to me that the previous neutrino normalization was in the region of stability for KenCarp4, but outside that of Rodas.

I wonder if CAMB/CLASS are using something higher-order for the photon boundary condition. Great find for the agreement at low cutoff -- I think this points to the culprit being some sort of spurious oscillations caused by the upper boundary condition.

[jmsull]

Yeah I tried to visualize the regions of stability when this issue first came up but gave up relatively quickly. I remember doing this for BDF methods in a class but the methods here are a bit more complicated =)

Re the photon BCs:

I need to read section 4 of the 2nd class paper more carefully but it appears CLASS switches on the radiative streaming approximation (RSA) for photons as the default after a certain time to mitigate the \ell max issue.

CAMB also does a late radiation truncation (in matter domination) which is described in Section XI of the CAMB notes, and actually appears to be the method prposed in CLASS (RSA). The notes mention this reduces the errors related to small scale oscillations due to the boundary (ell max) as you mention.

We can open an issue for adding this (or something like it)?

[jmsull]

Previous issues aside I started on the linear matter power.

Here is a first version of linear theory power and ratio (in scripts/first_plin.jl, still a little messy). In the ratio you can see the agreement is ~5% and scale-dependent. For the matter power CLASS uses gauge-independent variables and I have applied the gauge transformation to the transfer functions used for the matter power spectrum.

Here are some plots of transfer function ratios between Bolt and CLASS at z=0 where the differences can be seen more easily. Note the x-axes are not the same. There seems to be some broad overestimate of 5% in Bolt at all k, and then some smaller scale issues with the shape (which look due to baryon wiggles and I guess are due to issues with earlier evolution). The shape of the difference in transfers is similar to the one for plin as above so I would think fixing this should also fix plin. Of course the photon issue is visible here as well (I am using ell_max=100), and a similar issue is seen in the massless neutrinos (though not as bad). These perturbations are in Newtonian gauge, before the gauge change. I will try and resolve the scale-independent offset, which I guess may be related since there is no scale-dependent offset in plin.

[codecov-commenter]

Codecov Report

Merging #44 (19af79f) into main (6bee19a) will increase coverage by 2.97%. The diff coverage is 29.41%.

@@            Coverage Diff             @@
##             main      #44      +/-   ##
==========================================
+ Coverage   43.75%   46.72%   +2.97%     
==========================================
  Files           7        7              
  Lines         544      672     +128     
==========================================
+ Hits          238      314      +76     
- Misses        306      358      +52     

Impacted Files	Coverage Δ
src/Bolt.jl	100.00% <ø> (ø)
src/perturbations.jl	0.00% <0.00%> (ø)
src/spectra.jl	0.00% <0.00%> (ø)
src/ionization/ionization.jl	22.03% <25.00%> (+22.03%)	:arrow_up:
src/background.jl	86.36% <88.46%> (-2.21%)	:arrow_down:
src/ionization/recfast.jl	91.78% <100.00%> (+12.09%)	:arrow_up:
src/util.jl	98.38% <0.00%> (+12.90%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 6bee19a...19af79f. Read the comment docs.

[xzackli]

The RSA implementation actually looks very simple, and I think I can get to an implementation and tests by the end of this week. I've been busy coaxing some papers (needed for my thesis) to publication, so I might be a little slow on this.

[jmsull]

Large-scale/early-time mismatches in the above plots are now gone (resolved by fixing ICs #52). I will provide updated versions once RSA #48 or something like it is working.

[jmsull]

Last update possible before fixing RSA - wanted to see what the effect of fixing ICs was on Pk and the large-scale transfer function deviations seen above. Here I have turned off the radiation evolution at the RSA switch, even though it is not being evolved at all (which is wrong) it is less wrong than what it was doing before in the hierarchy due to the boundary condition issue.

Pk and all the non-relativistic transfers at z=0 are now in sub-percent agreement (minus one noise spike in \Phi around first bao wiggle, and high-k issues due to no RSA and low number of multipoles - using 50,50,20 here) with CLASS!

[jmsull]

At long last it looks like there is generally percent-level agreement on perturbations at all x(!). At least for the k mode (k~0.03 h/Mpc) I have been using for comparison the previous issues in the disagreement on the perturbations as a function of x are essentially resolved. Higher k modes need higher ell_max for the massless species.

After adding RSA (#48), adding reionization (#51), and fixing a missing He term in the tau' function (972fb7a), the late-time evolution of the massless monopoles and dipoles (and baryon perts) agrees well with CLASS. I just pushed a cleaned up version of the plot_perts_x script that uses a quickly wrapped rsa verson of boltsolve - using this I produced the following residual plot:

ell_max for neutrinos at 50 seems sufficient for percent-level accuracy (the residual wiggles are I think due to UFA in CLASS). The caveat to percent level accuracy comes in for the photons. The photon monopole is percent-level accurate except for whewhene there are spurious wiggles prior to RSA. I am not sure this is a Bolt problem - it may be that using the default value of CLASS ell_max is too low and it is producing this feature (and similar features are shown in the CLASS paper in their Fig. 6). I don't think the large-amplitude oscillation is physical, all the wiggles should be contained in a decaying envelope like we see in Bolt:

With these changes the linear power also looks better - the tau' issue was likely causing the large residual for the wiggles. The high-k drop at 2% ish seen in a previous version of this plot goes away when relativistic ell_max's are increased (I tried 500).

To achieve a similar level of agreement with CLASS that CLASS attains with CAMB (<0.1% in CLASS 3 paper), we would need to remove this ~0.5% offset seen on all scales. Beyond this such a comparison will involve lots of fiddling with precision parameters - I think this can wait for later?

[jmsull]

Spent some time enjoying the opportunity to learn about Julia's typing system - was able to restore the differentiability through recfast and ionization history that was breaking before - can get a derivative out of plin! Haven't checked the CMB example but the derivative with ForwardDIff should work there through recfast now also.

[jmsull]

CAMB comparison on pk for completeness - not sure what this 0.5% deviation is between CAMB and CLASS, but I have not been careful with the tolerances etc.

[jmsull]

@xzackli It looks like you already did some of the work here for a more precise comparison - I will come back to this using your precision parameters (but with neutrinos) later =)

[xzackli]

@jmsull oh god I wrote that notebook as a first year grad student and people keep contacting me about it. Warning: It's a bit out of date because Antony updated the CAMB interface to be nicer.

[xzackli]

This is more of a summary for myself, to wrap my head around the massive amount of work in this PR! This PR adds massive neutrinos and achieves a matter power spectrum within 0.5% of CAMB/CLASS.

Massive neutrino monopole normalization seemed to be off. Something really bad happened to the perturbation at x = -5, or z = 150. This was happily resolved by not dividing the massive neutrino perturbation by the massless neutrino background energy density, following MB.

There was also a problem with the ICs of the metric perturbations. Phi != psi in the presence of neutrinos. I still don't quite understand the symptom involving Neff. Again, following MB fixes the problem.

The photon monopole showed a deviation from the CLASS perturbation at around x = -6. It can be fixed by increasing the truncation multipole \ell_\gamma. We thought this might also be fixed by RSA, which would fix unphysical boundary conditions reflecting off of the truncation boundary. Indeed, it was! Currently the code overwrites the solution with the RSA solution, but eventually this could be incorporated using some sort of DifferentialEquations.jl hook.

Finally, this PR adds reionization to the ionization history, as well as fixing some problems with differentiability from pushing through the dual types. It has split the ionization history functions a little -- this could use some cleanup in the future.

With these changes, P(k) shows a 0.5% deviation with CLASS. I'd be interested in seeing these comparison plots using the precision parameters called CLASS:02 in CLASS Paper III, because those are the parameters which achieve 0.1% agreement with CAMB (under parameters CAMB:06).

[xzackli]

Something that isn't totally clear to me -- can you turn off the UFA in CLASS for applies-to-applies comparison? What do you mean by "CLASS uses a fluid approximation by default", is this referring to UFA?

[jmsull]

Something that isn't totally clear to me -- can you turn off the UFA in CLASS for applies-to-applies comparison? What do you mean by "CLASS uses a fluid approximation by default", is this referring to UFA?

Turning off UFA to compare the massless neutrino perturbations is something we should definitely do to follow up on the teal wiggles in this plot - I believe this can be done via the CLASS precision file. To really get apples to apples on everything it may be good to turn off TCA in CLASS as well, since we are not using that.

No, I mean they use a different fluid approximation for the massive neutrinos at late times (as described in CLASS IV). I turned this on and off in the plot the quote you gave is attached to.

[xzackli]

I gave the code one more look over, I'm merging it!!