lesgourg / class_public

Public repository of the Cosmic Linear Anisotropy Solving System (master for the most recent version of the standard code; GW_CLASS to include Cosmic Gravitational Wave Background anisotropies; classnet branch for acceleration with neutral networks; ExoCLASS branch for exotic energy injection; class_matter branch for FFTlog)
225 stars 282 forks source link

numerical instability with newtonian gauge? #498

Open alexsobotka opened 1 year ago

alexsobotka commented 1 year ago

Hello all!

I came across an interesting issue and I was wondering if anyone else has seen something similar.

As I understand it, the precision parameter _start_large_k_at_tau_h_over_tauk sets when CLASS begins the evolution of a perturbation mode. I am interested in a decaying particle similar to decaying dark matter (dcdm) and would prefer to initialize mode evolution before the decay. However, I discovered that sufficiently small values of _start_large_k_at_tau_h_over_tauk (i.e. initializing a mode earlier in time) leads to weird results in Newtonian gauge.

I am attaching some plots demonstrating this for basic LambdaCDM with a minimal assumption for ncdm. These plots were created with a newly installed un-modified CLASS v3.2 and the input settings:

For Newtonian gauge {'k_output_values': 0.003, 'gauge': 'newtonian', 'N_ncdm':1, 'm_ncdm': 0.06, 'N_ur': 2.0308, 'output':'tCl, mPk'}

Newtonian_wCls2

For synchronous gauge {'k_output_values': 0.003,
'N_ncdm':1, 'm_ncdm': 0.06, 'N_ur': 2.0308, 'output':'tCl, mPk'}

Synchronous_wCls2

In these figures, the dark blue line with _start_large_k_at_tau_h_over_tauk = 0.07 is the default value set by CLASS. As you can see, the Bardeen potential, Phi, has some very weird behavior in Newtonian gauge if _start_large_k_at_tau_h_over_tauk is sufficiently small so that the mode evolution starts well before horizon entry. The vertical dotted lines in the top panel show when the mode evolution begins. To my knowledge, there is no physical reason why we should see Phi growing in Newtonian gauge. More importantly, this has a very strong impact on the TT spectrum (bottom panel).

If the calculations are done in Synchronous gauge, this problem no longer occurs (second figure).

As far as I can tell, this behavior is true for all k modes. Additionally, the problem persists (albeit not as egregiously) even without the inclusion of ncdm:

{'k_output_values': 0.003, 'gauge': 'newtonian', 'output':'tCl, mPk'}

I’m perfectly happy working in Synchronous gauge, but this seemed to be a pretty significant deviation so I wanted to see if this has been addressed before. Thank you so much for your help!