Protons losing all their energy after Neutron decay

[antoniocapanema]

Dear authors,

While performing 1D simulations of UHECR propagation, I encountered a strange behaviour of the code whenever neutrons are produced. For simplicity, consider the emission of $N$ neutrons at some far away source. One would expect that all neutrons decay into protons and we observe $N$ protons at the Earth. However, the observer at the Earth is detecting $\sim 0.95N$ protons. This does not happen whenever I inject protons instead of neutrons (it can possibly still happen - albeit much less often - if neutrons are produced in $p\gamma \to n\pi^+$, followed by the decay of that neutron).

Here's an example code where this happens:

from crpropa import *

electrons = False
photons = False
neutrinos = False
antinucleons = True

sim = ModuleList()
sim.setShowProgress(True)
sim.add(SimplePropagation())
sim.add(Redshift())
sim.add(MinimumEnergy(0.01*TeV))

source = Source()
source.add(SourceParticleType(nucleusId(1, 0)))
source.add(SourcePowerLawSpectrum(0.01*EeV, (10**3.5)*EeV, -1))
source.add(SourceUniform1D(100*Mpc, redshift2ComovingDistance(8)))
source.add(SourceRedshift1D())

sim.add(ElectronPairProduction(CMB(), electrons))
sim.add(ElectronPairProduction(IRB_Gilmore12(), electrons))
sim.add(PhotoPionProduction(CMB(), photons, neutrinos, electrons, antinucleons))
sim.add(PhotoPionProduction(IRB_Gilmore12(), photons, neutrinos, electrons, antinucleons))
sim.add(PhotoDisintegration(CMB(), photons))
sim.add(PhotoDisintegration(IRB_Gilmore12(), photons))
sim.add(NuclearDecay(electrons, photons, neutrinos))

# Detects all particles at all points during the trajectory
obs = Observer()
obs.add(ObserverDetectAll())
output = TextOutput('test-all.txt', Output.Event1D)
output.disableAll()
output.enable(Output.SourceEnergyColumn)
output.enable(Output.CurrentEnergyColumn)
output.enable(Output.CurrentIdColumn)
output.enable(Output.SourcePositionColumn)
obs.onDetection(output)
obs.setDeactivateOnDetection(False)
sim.add(obs)

# Detects all particles at z=0
obs2 = Observer()
obs2.add(ObserverPoint())
output2 = TextOutput('test-all2.txt', Output.Event1D)
output2.disableAll()
output2.enable(Output.SourceEnergyColumn)
output2.enable(Output.CurrentEnergyColumn)
output2.enable(Output.CurrentIdColumn)
output2.enable(Output.SourcePositionColumn)
obs2.onDetection(output2)
sim.add(obs2)

# Detects nuclei at all points during the trajectory
nuclobs = Observer()
nuclobs.add(ObserverDetectAll())
nuclobs.add(ObserverPhotonVeto())
nuclobs.add(ObserverElectronVeto())
nuclobs.add(ObserverNeutrinoVeto())
nucloutput = TextOutput('test-nucl.txt', Output.Event1D)
nucloutput.enable(Output.SourcePositionColumn)
nucloutput.enable(Output.CreatedPositionColumn)
nucloutput.enable(Output.CreatedIdColumn)
nucloutput.enable(Output.CreatedEnergyColumn)
nuclobs.onDetection(nucloutput)
nuclobs.setDeactivateOnDetection(False)
sim.add(nuclobs)

# Detects nuclei at z=0
nuclobs2 = Observer()
nuclobs2.add(ObserverPoint())
nuclobs2.add(ObserverPhotonVeto())
nuclobs2.add(ObserverElectronVeto())
nuclobs2.add(ObserverNeutrinoVeto())
nucloutput2 = TextOutput('test-nucl2.txt', Output.Event1D)
nucloutput2.enable(Output.SourcePositionColumn)
nuclobs2.onDetection(nucloutput2)
sim.add(nuclobs2)

sim.run(source, 100, True)

This led me to explore neutron-injection scenarios further and check what was happening during propagation (check for example the output files with ObserverDetectAll()). It seems that the proton is created from the neutron's decay, but it soon loses 100% of its energy and ends up with exactly zero energy. Thus, it's not detected by the observer at the Earth (since the minimum energy of the simulation is 0.01 TeV).

What could be causing this 100% energy loss? Activating/deactivating some interaction modules, I found that this bug(?) stops happening when ElectronPairProduction is not included in the simulation. So it seems to be caused somehow by Bethe-Heitler pais production right after the proton is produced.

This also leads to loss of nuclei at z=0 when running simulations with heavier nuclei. E.g.: sometimes, He nuclei emission, only two protons will arrive at the Earth. My guess is that in the decay channel $He^4 \to ... \to 2p+2n$, the protons that should come from the decay of the 2 neutrons are also missing the the same reason.

System:

• Operating system: Ubuntu Linux 22.04
• Computer (remotely accessed)
• Python version: 3.10.12
• Swig version: 4.0.2
• Cmake version: 3.22.1

[mohller]

I haven't run your script but it looks correct. The simulated situation and the results you describe are consistent. The cosmic rays simulated span a minimum of 100 Mpc and extend to ~4Gpc, which cover the characteristic distances for pair production and photodisintegration losses (which your simulation includes). Therefore, there should be a significant amount of energy loss for the cosmic rays. I am not sure how you obtain the 0.95 value, but it seems you are not accounting for the energy losses, but just counting the particles. However, the results you describe seem reasonable for the simulated scenario.

[antoniocapanema]

Thanks for your prompt reply.

My apoligies, perhaps I was not clear in explaining what the problem was. Let me explain it in another way: If I inject 10 protons using the above code, I get 10 protons at the Earth every single time, as expected. If I inject 10 neutrons using the above code, one would expect 10 protons at the Earth. Instead, I only get only 9 protons (sometimes 8, sometimes the full 10) in the output file. I found that these protons are missing because, after the neutrons decay into the protons, some protons will lose all of their energy and go literally to $E=0$.

Take a look at an example of this happening below. In this example, I track a neutron since its emission. Note the last two lines, where the proton from the decay of that neutron is produced at 3.197 EeV and then goes straight to zero energy in the next line. What's causing this total energy loss of the proton?

D               ID              E                E0              X0              ID1             E1              X1
1.00000E-04     1000000010      3.20152E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03
1.09929E-02     1000000010      3.20151E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03
2.18857E-02     1000000010      3.20150E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03
3.27785E-02     1000000010      3.20149E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03
4.36712E-02     1000000010      3.20148E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03
5.45639E-02     1000000010      3.20148E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03
6.54565E-02     1000000010      3.20147E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03
7.63490E-02     1000010010      3.19705E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03
2.92830E+02     1000010010      0.00000E+00      3.20152E+00     6.24176E+03     1000000010      3.20152E+00     6.24176E+03

[mohller]

The difference in propagation distance in the last two lines is of ~300 Mpc which is larger than the energy loss length for protons of such energy ~3 EeV (see for example figure 2 in this paper). The energy loss is consistent with what it is expected. You can reduce the propagation step in order to resolve the energy loss. To do this, use SimplePropagation(minstep, maxstep) and provide suitable values for the parameters minstep, maxstep.

[antoniocapanema]

I see! It seems that decreasing the maxstep in SimplePropagation() has indeed resolved the issue. Thank you so much for your help (and for the reference)!