CRPropa / CRPropa3

CRPropa is a public astrophysical simulation framework for propagating extraterrestrial ultra-high energy particles. https://crpropa.github.io/CRPropa3/
https://crpropa.desy.de
GNU General Public License v3.0
68 stars 67 forks source link

Problematic definition of the spectral index in PlaneWaveTurbulence #387

Closed LeanderSchlegel closed 2 years ago

LeanderSchlegel commented 2 years ago

Dear all,

we noticed a problematic definition respective a bug regarding the spectral index treatment in PlaneWaveTurbulence and opened this issue to ask you if we see it correctly and formally report the issue. If we should change it, it should only be a small modification of the handling of the spectral index in the source code.

Definition of the issue Using SimpleTurbulenceSpectrum just producing a powerlaw, the default value for the spectral index of 5/3. seems not to result in the expected -11/3. spectrum as shown in the following plot, that reweights the observed spectrum with +11/3. Instead the argument of -1/3. yields the -11/3. spectrum.

spectrum

Minimal example to reproduce the issue Attached you find an example code to reproduce the plot, based on a script by @antoniusf, that generates the turbulent field, samples it onto a grid to perform a FFT gaining a spectrum that gets reweighted with +11/3.

example.txt

Thank you!

rafaelab commented 2 years ago

Hi @LeanderSchlegel, Maybe @adundovi can tell a bit more about the conventions he adopted for the definition of the power spectrum.

adundovi commented 2 years ago

Hi @LeanderSchlegel, thanks @rafaelab,

yes, it is true that there could be a mistake in calculating the correct index, I have noticed the same in PlaneWaveTurbulence, but didn't have time to play around with it to double check it. The internal conversion should convert it from 5/3 to -11/3.

antoniusf commented 2 years ago

Hi, just a quick note from me: The problem appears to be in TurbulenceSpectrum::energySpectrum, where the energy spectrum is computed as

This is in contrast to the formula given by Tautz and Dosch (2013), which as you know PlaveWaves is implemented after, who explicitly state that for their formula, Kolmogorov turbulence is achieved by setting s to 5/3:

Note the missing "+1" in the exponent – this explains why the correct spectral index is obtained when passing in -1/3 instead of 5/3.

While it's true that the grid turbulence models need a spectral index of 11/3 due to the way they're implemented (and this is what the extra "+1" achieves when given 5/3), the TD13 model does not need this extra correction and hence outputs the wrong spectrum when passed an 11/3 spectrum. (Which again is what energySpectrum produces when given a spectral index of 5/3.)

adundovi commented 2 years ago

Hmm, OK, I'll take a look how to make it consistent since the same spectrum interface should apply to both, the grid and the plane-wave methods.

adundovi commented 2 years ago

@antoniusf @LeanderSchlegel, O.K., I agree that we should stick to the second formula as it is a text-style definition (I also use it in my 2020 turbulence paper, Eq. 20) and add the "+1" correction in GridTurbulence. I am really sorry for not reviewing and replying sooner to this issue. I cannot commit the change in next few days, so if you are able, please be free to do so.

LeanderSchlegel commented 2 years ago

Thanks for your ok! We prepared a small fix changing the +1 and found that it would be uniform to also change the argument of energySpectrum in TurbulentField.h from khat to k, putting the multiplication with l_bendover into the function. We tested that it works and could prepare a pullrequest in the coming days :)