Wrong ICS spectrum of an electron

[alexkorochkin]

Dear Developers,

I was testing the EM module and found a bug with inverse Compton scattering of an electrons. The simplest case of ICS of an electron for which the spectrum of secondaries can be easily calculated is the upscattering of CMB photons by monochromatic electrons while all other backgrounds are turned off.

For this case I calculated the theoretical prediction for the ICS (red dashed lines on the plot) an compared it with the output from CRPropa (green solid lines) for 5 different initial energies of electrons: 1e9, 1e10, 1e11, 1e12 and 1e13 eV. Simple estimates of a typical secondary photon energy are shown with gray vertical lines (Esec = 4/3 E_cmb (E_electron/m_electron)**2, E_cmb=6.6e-4 eV). CRPropa predicts systematically higher energy of secondaries which can't be attributed to the statistical fluctuation.

After checking the EMInverseComptonScattering.cpp file, I realized that it works passably and the error is in the precomputed tables located at CRPropa/build/data/EMInverseComptonScattering: 'cdf_CMB.txt' and 'rate_CMB.txt'. The clear error in 'rate_CMB.txt', when the interaction rate depends on energy in Thomson regime, increased my suspicions. So I recalculated the files for the CMB myself and put it instead of original ones. The result after correction are shown with blue lines, the problem disappears.

The error in precomputed tables arises from bad integration with too big steps. The only change that I made was using smaller integration step sizes. Here I attached the new table 'cdf_CMB.txt' computed for the energies 1e9-3e13 eV. I used constant Compton scattering cross section valid in this energy range.

I didn't check but I expect that the same error exists for all other backgrounds however the effect of bad integration for EBL is smaller because the EBL bumps are wider so the integration with big steps gives more precise result.

Alexander Korochkin

CRPropa3-3.1.7

import crpropa from crpropa import *

simulation = ModuleList() simulation.add(SimplePropagation(1pc, 1kpc)) CMB_instance = CMB() simulation.add(EMInverseComptonScattering(CMB_instance, True, 0, 0.01))

Ee=1e0GeV # energy of electrons minThreshold = MinimumEnergyPerParticleId(1eV); minThreshold.add(11, Ee-1e-8GeV); # kill electron after first interaction minThreshold.add(-11,Ee-1e-8GeV); minThreshold.add(22, 1*eV); simulation.add(minThreshold)

obs = Observer() obs.add(ObserverPoint()) txtout = TextOutput("CRpropa_ICS_1_GeV.txt", Output.Event1D) txtout.set(Output.WeightColumn, True) txtout.set(Output.SourceIdColumn, False) txtout.set(Output.SourceEnergyColumn, False) txtout.set(Output.CreatedIdColumn, True) obs.onDetection(txtout) simulation.add(obs)

source = Source() source.add(SourcePosition(Vector3d(10,0,0)*Mpc)) source.add(SourceParticleType(11)) source.add(SourceEnergy(Ee))

simulation.run(source, 100000, True, True)

cdf_CMB_corrected.txt

[lukasmerten]

Hi @alexkorochkin Thanks for your detailed bug report. I will take a look at it later this week. When I can confirm it, we will include a fix for the upcoming CRPropa 3.2 version. You said the bug was connected to a too sparse integration of the interaction rate, wasn't it? So which supporting points or which integration routine did you use?

[rafaelab]

Thanks for the report, @alexkorochkin. @lukasmerten: Maybe this was introduced in 3.1.6 together with the introduction of custom photon fields? At some point I compared the implementation with the theoretical prediction and other codes, and they showed good agreement. This makes me think that the bug could have been introduced in changes to the scripts in CRPropa-data (which, btw, are not up-to-date).

[alexkorochkin]

@lukasmerten, I used the simplest trapezoidal method with 125 logarithmic bins per decade for s integration in range 1e1-1e13 eV^2 and 250 logarithmic bins per decade for the integral with respect to E_cmb in range 1e-5-1e-2 eV.

[rafaelab]

Hi @alexkorochkin,

I spent quite some time working on this issue in the past few days. There was certainly something weird going on with the rates in the Thomson regime. This was indeed just the integration steps. I'm preparing a pull request to fix this now. However, I cannot fully reproduce your results.

• Using the file you provided our results match.
• I tried increasing the integration steps and I still don't get a peak at the same position as you, not even with 10000 log bins per decade.

How did you calculate the rates? Did you use the scripts from the our CRPropa3-data repository? I am under the impression that in the CDF file you sent a different distribution for the CMB was used (maybe its derivative?).

After improving the integration steps and calculating the average and peak energies of the distribution I get the expected result. What you refer as the peak in your issue report is actually the average of the distribution, not its peak. Under this assumption things seem fine, although the integration steps indeed needed to be fixed. Check out the plots below for 100 GeV and TeV electrons.

[alexkorochkin]

Hi @rafaelab,

I calculated rates using my own script, because I didn't know CRPropa scripts are available online. In my previous post, by "simple estimate of a typical secondary photon energy" \<E> I meant the expected value of E given the distribution dN/dE:

\<E>=integral( E(dN/dE)dE )/integral( (dN/dE)*dE )

From theory we know that \<E>=4/3 (E_electron/m_electron)^2. For 1 TeV electrons \<E> ~= 3.2 GeV. This value we can compare with it's estimate from the CRPropa simulation. To avoid any ambiguities with the data I took the simulation points directly from your plot for 1 TeV electrons. More precisely, I took two curves - solid black ("corrected") and dotted blue ("issue corrected"). It seems to me that you have a typo in the y axis label, there probably should be EdN/dE instead of E^2 dN/dE. So, in further analysis, I thought that you actually plot E\dN/dE. Files with the points I attached below. Using these points I numerically calculated the integrals which I wrote above for \<E>. Python script with the integrals I also attached below.

Results: \<E>("issue corrected") = 3.15 GeV \<E>("corrected") = 4.74 GeV

One can see that the estimate of E with the table that I provided in my first post \<E>("issue corrected") is very close to the theoretical value and \<E>("corrected") is much higher. In any case, you can compare the full distribution dN/dE from the original CRPropa with the theoretical one to make sure that they do not match.

I'm a bit confused by your orange vertical line ("corrected - average"). How do you calculate it?

If increasing number of integration steps doesn't help I don't know where this bug comes from.

Issue_script.txt CRPropa_issue_corrected.txt CRPropa_corrected.txt

[LeanderSchlegel]

Hi @alexkorochkin,

we tested the table for the photon density of the CMB recently and thought perhaps it could add to the discussion. We found that in CRPropa 3.1.7 this table had a bug that got fixed in the current version. In the attached plot the density tables in 3.1.7 and the current master (the upcoming CRPropa 3.2) are plotted as well as a self generated table by sampling the now analytically available CMB() function and a python implementation of that same formula. While 3.1.7 shows a different curve, all others lay exactly on top of each other.

CMB_test.pdf

[rafaelab]

Hi @alexkorochkin

During the development of new features for CRPropa 3.2, the TabularPhotonField class was added. I believe there was an error in the CMB implementation inheriting from it: the integral or derivative of the field was being used. This bug is obvious by analysing @LeanderSchlegel's plot. In fact, if you take the CMB peak wavelength (not frequency) of 1.063 mm, you obtain results for the plot I sent earlier.

In my previous comment, I used the correct CMB distribution for the tests, but the one used internally in the simulation was incorrect, hence the difference between our peaks. I didn't expect the CMB to have been affected by any changes, so I had completely ignored this possibility at the time.

This problem is unique to version 3.1.7. It was not there in 3.1.6 and earlier versions. It has already been corrected in the development version. For your comparisons, I encourage you to use the current development version (which is very close to the one we'll use in the next release, 3.2), or the 3.1.6 release, not the 3.1.7. Also the rates for inverse Compton at low energies (<100 GeV) were problematic.

If you are benchmarking your code against CRPropa, could you please make sure you use the correct version (not 3.1.7) and let us know if the problem remains? I ran a quick test and your predictions are in agreement with CRPropa's.

You can just git clone the current version in the repository or wait for the release of 3.2 (next week, probably).

Once you confirm that the problem was fixed, I will close this issue.

[rafaelab]

Hi @alexkorochkin As I said, I can reproduce the correct results with the data you sent. I'm still struggling to identify the source of the discrepancy, if it is the CRPropa-data repository, or if it is the CRPropa code itself. The findings of @LeanderSchlegel make it even more puzzling to me because I cannot find an obvious change that could have broken things between 3.1.6 and 3.1.7. Have you found something specific? If yes, please report to us through this issue so that we can further investigate. Many thanks!

[alexkorochkin]

Hi @rafaelab

I think I identified the origin of the problems in CRPropa. I checked CRPropa-data rep and focused on the script which calculates interaction rate and cdf (interactionRate.py). I found that if I replace calc_rate_s function with my own, the output tables become the same as the tables I calculated myself without CRPropa-data scripts (see attached plot for the inverse Compton scattering table for 100 GeV electrons). Moreover I found that the problem with original calc_rate_s function is one lost term when calculating cdf. I've attached the pdf file with the detailed description of the problem with this term. interaction_rate.pdf

I checked that the problem exists for all available CRPropa versions (3.1.5, 3.1.6, 3.1.7, master) and also affects ICS on EBL. The wrong cmb table that @LeanderSchlegel found is not used in ICS.

Could you please cross-check my results? If I'm right then almost all interactions in CRPropa have wrong energy distribution of secondaries because interactionRate.py is used by many other interaction scripts.