Incorrect Lensing Example on Documentation Page?

[angelinapartenheimer]

Hello,

I have tried to reproduce the example provided at https://crpropa.github.io/CRPropa3/pages/example_notebooks/galactic_lensing/lensing_cr.v4.html but find very different results.

I have generated a magnetic lens using the backtracking simulation and createlens.py code provided at https://crpropa.github.io/CRPropa3/pages/example_notebooks/galactic_backtracking/galactic_backtracking.v4.html#Backtracking-to-Generate-a-Lens on the documentation page. I have not modified the provided code in any significant way.

I am running the exact code provided on the documentation page:

for i in range(1000):
    particleId = crpropa.nucleusId(1,1)
    energy = 10 * crpropa.EeV
    galCenter = crpropa.Vector3d(-1,0,0)
    momentumVector = crpropa.Random.instance().randFisherVector(galCenter, 200)
    M.addParticle(particleId, energy, momentumVector)

lens = crpropa.MagneticLens(''../lens_data/lens128/lens.cfg'')
lens.normalizeLens()
M.applyLens(lens)

crMap = np.zeros(196608)
for pid in M.getParticleIds():
    energies = M.getEnergies(int(pid))
    for i, energy in enumerate(energies):
        crMap += M.getMap(int(pid), energy * crpropa.eV )

healpy.mollview(map=crMap, title='Lensed')

Both my test and the documentation example begin with a cluster of protons coming from the direction of the galactic center at (x, y, z) = (-1, 0, 0).

Figure 1: An unlensed cluster of 10 EeV protons coming from the direction of the galactic center. This cluster is used in the example on the documentation page and in my lens test. Note that the protons must pass through the galactic center to reach the observer at the solar system.

The documentation page shows most of the 10 EeV protons passing through the galactic center and reaching the observer with moderate lensing effects:

Figure 2: The lensed cluster of 10 EeV protons shown on the documentation page. These protons are coming from the direction of the galactic center, which means they must pass through the galactic center to reach an observer at the solar system. Most protons in this example are still able to reach the observer.

Using my lens, however, almost no 10 EeV protons are able to reach the observer after passing through the galactic center:

Figure 3: The effect of my lens on 10 EeV protons passing through the galactic center. Almost no protons are able to pass through the strong coherent magnetic fields at the center of the galaxy and reach the observer at the solar system.

To demonstrate that my lens is able to handle a cluster of 10 EeV protons in other cases, here is an example of the same cluster of protons coming from the direction opposite to the galactic center, with my lens applied to it. The difference in resolution compared to the documentation comes only from my using nside = 128 instead of the default nside = 64. While this example is not shown on the documentation page, it shows that my lens is able to handle a cluster of 10 EeV protons without any issues with normalization/binning:

Figure 4: My lens applied to a cluster of 10 EeV protons approaching from the direction opposite the galactic center (x, y, z) = (1, 0, 0). While this example is not shown on the documentation page, it demonstrates that my lens is able to handle a cluster of 10 EeV protons without any normalization or binning issues. Note that these protons do not have to pass through the galactic center to reach the observer, and therefore do not encounter the strong coherent magnetic fields present there.

The example provided on the documentation page seems physically inaccurate. A simulation of a 10 EeV (or lower) protons passing through the center of the galaxy shows that the proton experiences very strong deflections, and is highly unlikely to reach an observer on the other side of the galactic center. I have provided some examples of these trajectories: in my plots, the red dot is the galactic center, and the green dot is the solar system observer. The proton approaches the galactic center along the positive x-axis.

Figure 5: Trajectories of 2, 5, and 10 EeV protons approaching the galactic center. The protons encounter the strong coherent magnetic fields at the galactic center (red dot) and are fully deflected away from an observer at the solar system (green dot).

The heavy deflections of the <10 EeV proton are due to the regular component of the galactic magnetic field. Near the center of the galaxy, for r<5 kpc, the ring, disk, and other regular magnetic field components of the JF12 model have values ranging from 0.1-4 muG. For a 10 EeV proton in a 1 muG magnetic field, the gyroradius is ~10.8 kpc, and in a 4 muG magnetic field it is ~2.7 kpc. This is sufficient to deflect the proton by 90 degrees within the 10 kpc range of the galactic center where these strong coherent magnetic fields exist.

Because of this, I believe the example on the documentation page indicates 10 EeV protons passing through the galactic center with insufficient deflection. The example is irreproducible given the lens building code provided on the documentation page. Could you give a check and confirm this finding?

Thank you,
Angelina Partenheimer

[lukasmerten]

Hi @angelinapartenheimer

Thanks for the detailed description of your issue.

I am not sure if that's part of the problem but the color scales of your figs. 2. and 3. are very different; maximal value 4e-4 in fig 2 and 5 in fig 3.

When you generated your lenses did you use a Galactic magnetic field model including a random component or only a coherent large scale field. In the case of a mixed (regular/coherent + turbulent) magnetic field I would expect that particles can also propagate through the Galactic center. At that stage it is only an educated guess from my side but I currently have no time to investigate this in more detail.

However, we plan to update the documentation and examples in the weeks and will definetely take a closer look at this issue then.

Lukas

[angelinapartenheimer]

Hi @lukasmerten

The color scales on my plot are different from the ones listed on the documentation page because I did not apply normalizeLens(), while the documentation page did. This was my oversight, I apologize for being inconsistent. However, I know this is not the root of the problem, as applying normalizeLens() does not change the appearance of the skymap. I reran the code with normalizeLens() applied and attached the new plots below. As you can see, the output is the same (aside from slight differences, as the random cluster I am lensing is slightly different every time). I will admit that my maximum pixel value is still much higher than that shown on the documentation page. I am normalizing the lens and plotting with code copied and pasted from the example at https://crpropa.github.io/CRPropa3/pages/example_notebooks/galactic_lensing/lensing_cr.v4.html, as you can see:

lens = crpropa.MagneticLens(‘../lens_data/lens128/lens.cfg')
lens.normalizeLens()

M.applyLens(lens)

crMap = np.zeros(196608)
for pid in M.getParticleIds():
    energies = M.getEnergies(int(pid))
    for i, energy in enumerate(energies):
        crMap += M.getMap(int(pid), energy * crpropa.eV )

hp.mollview(map=crMap, title='Lensed: 10 EeV, passing through galactic center', cmap=‘jet')

I’m not sure why the results are different, but this is an issue I am also hoping to resolve.

I generated my lens using both random and coherent components of the magnetic field. In fact, the code I am using to set up the magnetic field is exactly the same as that shown in the example at https://crpropa.github.io/CRPropa3/pages/example_notebooks/galactic_backtracking/galactic_backtracking.v4.html. The only difference is that the documentation example uses initTurbulence(), which is deprecated. Instead, I just use B.randomTurbulent(). Based on reading the source code, this seems to automatically implement the same magnetic field as in the backtracking run on the documentation page (see lines 102-114 at https://github.com/CRPropa/CRPropa3/blob/master/src/magneticField/JF12Field.cpp.) So you can compare, here is the code I used to generate my magnetic field:

def backtrackingRun(logE=18., nside=128):
    
    # magnetic field setup
    B = JF12Field()
    seed = 1703202123
    B.randomStriated(seed)
    B.randomTurbulent(seed)

    #simulation setup
    sim = ModuleList()
    sim.add(PropagationCK(B, 1e-4, 0.1 * parsec, 100 * parsec)) 

I know the turbulent component of the field is up and running because I can turn it off and see a visible difference in the lens I generate: the overall lensing effect is the same, but the blob of particles is less scattered. However, if you think there is something wrong with my magnetic field setup, or if the parameters I am using are somehow incorrect, please let me know!

I will continue checking the documentation page for the updates.

Thank you,

Angelina Partenheimer