> 1% difference between class v2.9.4 and class v3.2.0 for lensing result at l>3000? Loss of accuracy?

[borisbolliet]

Hi,

I was upgrading class version for some data analysis that requires high accuracy at high ell.

I find significant differences for lensed TT, TE and EE spectra and for lensing potential power spectrum, when I run both versions in the exact same settings.

The settings are:

common_settings = {# LambdaCDM parameters
                   'h':0.67810,
                   'omega_b':0.02238280,
                   'omega_cdm': 0.1201075,
                   'A_s':2.100549e-09,
                   'tau_reio': 0.05430842,
                   'n_s': 0.9619} 

M = Class()
M.set(common_settings)
M.set({'output':'tCl,pCl,lCl,mPk',
       'modes':'s',
       'lensing':'yes',
       'l_max_scalars':l_max_scalars,

       'non linear': 'hmcode',
       'eta_0' :  0.603,
       'c_min' : 3.13,
       'N_ur' : 3.046,

       'accurate_lensing': 1, #necessary to get accurate high ell lensed

       'k_max_tau0_over_l_max':15., #necessary to get accurate high ell lensed
       'perturbations_sampling_stepsize' : 0.05, #necessary to get accurate high ell lensed
      })

M.compute()

Note that perturbations_sampling_stepsize is called perturb_sampling_stepsize in v2.9.

See plots attached. [PP: cl^phi-phi]

Class v2.9 matches CAMB better (at the 0.5% level down to l=9000 in those settings for all spectra, see #386 for comparison between CAMB and class v2.9.4 on cl^phi-phi). Hence it seems that Class v3 is less accurate in those settings.

Am I missing something?

If correct, what are the settings I need to set to class v3 to recover class v2.9 or CAMB results at better than 1% accuracy?

Thanks a lot! Boris

[lesgourg]

Thank you @borisbolliet , this is really useful, we will try to find some time to investigate this asap!

[borisbolliet]

One difference between class v2.9.4 and v3 is that the recombination treatment is done by recfast in v2.9 and hyrec in v3. This certainly explains the difference for the red curve (thanks Yacine Ali-Haimoud and @jcolinhill for noting that). I will check. However this can't be the reason why the lensing power spectrum, or the lensing of the cmb, is so different.

[schoeneberg]

Dear all, dear @borisbolliet ,

I was recently made aware of this issue, and I started investigating (given that I spearheaded a lot of the transition from 2.9.4 to 3.0.0 I think that it's only fair I look into issues like this). First, I want to report that this has little to do with the actual parts of the code I touched for the transition (or recombination, or any background, thermodynamics, or perturbation quantity). Indeed, if you set 'nonlinear':'none' and 'recombination':'recfast' (as per Yacine and @jcolinhill ), you will find that the difference is now vanishingly small. This leads me to conclude that the issue is in the way that halofit is used between 2.9.4 and the current v3. In particular, from reading through the relevant git diff, it seems like halofit itself didn't have any major change in this period. This can be confirmed by looking at the small difference in the nonlinear P(k) at the permille level (see figure below) As such, I will need some more time to dig deep into the nl_corr_density term that CLASS uses in order to track down this bug. Thanks for letting us know!

[schoeneberg]

Dear @borisbolliet (dear @lesgourg ),

I finally tracked down the "bug", except that it isn't really a but, but rather an issue of not increasing the precision high enough (at least that is my conclusion for now). When I adjust the parameter P_k_max_1/Mpc to 600, I get the following comparison Since that difference is sub-permille, I would claim excellent agreement between class v2.9.4 and class v3.

It appears the difference in your case was caused by the Halofit correction being computed or not computed at corresponding high redshift. Note that to compute accurately the lensed spectra for such high L, we are deep in the non-linear regime, even at relatively early times where the CMB lensing kernel is relevant. As such, it is very important to have good computations from e.g. halofit at high redshift and high k, but with the default k_max this isn't possible. I think the idea was to increase k_max by increasing k_max_tau0_over_l_max, but this is the k_max with a very high k-sampling which is even linear, so it's not quite the correct parameter for getting accurate non-linear predictions. For that, one should raise instead P_k_max_1/Mpc to some reasonably high value (e.g. 500 or 1000 or higher), in order to have well converged integrals for the nonlinear correction. Otherwise, the nonlinear correction will only be computed at relatively low redshift, and the precise redshift up until which it is computed differs from v2.9.4 to v3, leading to corresponding differences in the C_L^PP and the lensed CL at very high L.

Since halofit was not designed to accurately predict the high-redshift high-wavenumber regime, I would additionally like to include a great word of caution against using these results with too much certainty: In this regime the answer is probably not well determined, both from lack of calibration of halofit to simulations with super high k and high z, but also from non-inclusion of other physical effects (e.g. baryonic effects), which would probably dominate in any case.

In conclusion, Thank you for reporting this issue. The fix is to raise "P_k_max_1/Mpc" to some high level (~500) at which non-linear corrections can accurately be computed. I also caution it's dangerous to use halofit in this regime.

Let me know if there's anything else.

[jcolinhill]

Hi @schoeneberg , thanks for the investigation! Indeed this is a topic that Fiona, Mat, and I wrote a paper about: https://arxiv.org/abs/2103.05582

I suppose the puzzle is primarily why the default behavior changed so much from v2.9 to v3.

Also, for the analysis of current ACT and SPT data, such settings already have an impact (cf. Appendix A of https://arxiv.org/abs/2109.04451 ). Perhaps it would be good to increase the numerical convergence in the next release?

(I agree the question of the overall accuracy of halofit, etc is a separate issue.)

Thanks, Colin

[schoeneberg]

Dear @jcolinhill , thanks for weighing in! Let me answer a few points you raised:

• I think the default behavior did not change much (at all?) between v2.9 and v3, but if you attempt to compare non-well-converged results, you will almost always find large differences. I am instead rather surprised by #386 agreeing so well with CAMB despite in principle not ensuring convergence with the chosen precision settings. So far, beyond corrections in the thermodynamical quantities (which are now more consistent/accurate) I only found differences in the precise time steps used for the halofit computations, so I am not entirely sure why the difference is so drastic beyond just small errors being super-amplified by lack of numerical convergence in either case. I'd be happy about any input from somebody else investigating this, though :smile:
• I am not sure the point about ACT and SPT requirements is entirely true -- the final lensed CL are only used up until L ~ 4000, and the impact of the divergence in the PP spectra beyond L ~ 5000 is hopefully not felt by them, otherwise we should seriously put some doubt on any analysis relying too closely on the high-L data. If the impact really is large (which remains to be proven), then increasing the default precision is a possible step, but given the question of whether even a well converged halofit makes any sense in this regime, it's not a trivial decision.
• This is especially true since the runtime grows quite a lot if we start to include very large k ranges. One doesn't notice it for only massless neutrinos (as in these investigations), but as soon as massive neutrinos are present, there are many many oscillations to be followed. As such going up to k ~ 500/Mpc by default would make the code prohibitively slow even in many unrelated cases for normal users, who do not require such accuracy. All of this to potentially only "solve" a problem that is in any case plagued by a much larger uncertainty (of which nonlinear modeling to use).
• It's safe to say, though, that there should be more warnings about this, e.g. when the user requests lensed CL at very high L.

I am not sure what the best treatment here is (and the decision in the end is up to @lesgourg ). Probably it would be nice to accelerate the code in the large-k regime somehow, or to at least dynamically detect when more precision is required.

For now, there is a working hotfix :wink:

[borisbolliet]

Ah I am just seeing this ! Excellent, thank you so much @schoeneberg for taking time and finding the solution :) Closing the issue now :)