Caustic-like Structures in UHECR Flux after Propagation in Turbulent Intergalactic Magnetic Fields

We start with a simplified model of a source emitting protons isotropically in all directions. The entire volume around the source is filled with the homogeneous turbulent magnetic field. Let us also assume that the parameters of the magnetic field and the energy of the particles are chosen in such a way that the deflections after propagating a distance of several correlation lengths are small, but not negligible. In this model, one would expect that the distribution of UHECR arrival positions on a sphere of radius $D>\lambda_{C}$ would remain isotropic, since the deflections of particles initially aimed in different directions are statistically independent in a random magnetic field. However, as we show below, this is true only if the radius of the sphere is much larger than the correlation length. On the other hand, if the radius of the sphere is on the order of several correlation lengths, the distribution of particles on the sphere becomes essentially anisotropic at intermediate angular scales.

The physical reason for the appearance of inhomogeneities during propagation in a turbulent field can be understood analytically. In [6] it is shown that in the case of initially parallel proton beam, inhomogeneities in the particle distribution are caused by fluctuations of the magnetic field rotor along the trajectories of neighboring particles. Here we follow the logic of [6] and derived the modified formula for the flux amplification for the case of a divergent particle beam.

Consider two particles shifted by the vectors $\Delta\vec{x}_{i}$, $i=1,2$ relative to the reference trajectory $\vec{x_{0}}$. If the distance between particles is less than the correlation length, the evolution of the displacement vector $\Delta\vec{x}_{i}$ follows Equation 4.1 from [8]:

$$\frac{\mathrm{d}^{2}\Delta\vec{x}_{i}}{\mathrm{d}t^{2}}=\frac{Zec}{E}\left(\frac{\mathrm{d}\vec{x}_{0}}{\mathrm{d}t}\times\Delta\vec{B}_{i}(\vec{x}_{0})+\frac{\mathrm{d}\Delta\vec{x}_{i}}{\mathrm{d}t}\times\vec{B}(\vec{x}_{0})\right)$$

(2)

which is a direct consequence of the Lorentz equation. Here $Ze$ is the charge of the particle, $c$ is the speed of light and $\Delta\vec{B}_{i}(\vec{x}_{0})\equiv\vec{B}(\vec{x}_{0}+\Delta\vec{x}_{i})-\vec{B}(\vec{x}_{0})$. To calculate the change in the flux, consider two particles displaced relative to the reference in the perpendicular directions $\Delta y$ and $\Delta z$. The change in flux is inversely proportional to the change in area $A$ subtended by these displacement vectors. For an initially parallel beam of particles, the area between the particles turns out to be linearly dependent on the rotor of the magnetic field [6]:

$$A(D)=A_{0}\left(1-\frac{Ze}{E}\int\limits_{0}^{D}(D-s)\,(\mathrm{rot}\vec{B}\cdot\mathrm{d}\vec{s})\right)$$

(3)

where $s$ is a parameter along the trajectory.

In the case of a diverging beam, the equation (3) should be modified due to different beam geometry, which determines different initial conditions. We have used these initial conditions and have derived the equation for the diverging beam:

$$A(D)=A_{0}\left(1-\frac{Ze}{E}\int\limits_{0}^{D}s\left(1-\frac{s}{D}\right)(\mathrm{rot}\vec{B}\cdot\mathrm{d}\vec{s})\right)$$

(4)

It should be noted that the equations are valid only for small amplification $(A(D)-A_{0})/A_{0}\ll 1$ and small deflections of particles.

Equation 4 provides a direct prediction for the flux amplification that can be tested numerically. Indeed, if $N_{0}$ is the mean number of particles that have passed through a given pixel, then the number of particles $N$ that actually hit a given pixel corresponds to the amplification of the flux. Rewriting equation 4 in terms of the number of particles, we arrive to a relation that can be verified directly in numerical simulations:

$$\frac{N-N_{0}}{N_{0}}=\frac{1}{R_{L}}\left[\frac{\int\limits_{0}^{D}s\left(1-\frac{s}{D}\right)(\mathrm{rot}\vec{B}\cdot\mathrm{d}\vec{s})}{B}\right]\,,$$

(5)

where $R_{L}$ is a Larmor radius.

We use publicly available codes CRbeam [9] and CRPropa [10, 11] for cosmic ray propagation. All interactions and cosmological expansion of the Universe are turned off. The particles were emitted isotropically and propagated until they reached a sphere of a given radius $D$. For each particle, we store its initial direction, final direction, and final position on the sphere. The resulting output is processed with the package healpy [12, 13].

The use of two codes makes it possible to model a turbulent magnetic field in two different ways. In CRbeam the magnetic field is generated as a sum of plane waves with random phases and directions [14] and its exact value is recalculated on the fly before each next step during particle propagation. In CRPropa, on the contrary, the field is precalculated on the grid and its value is set by interpolation between grid points. In both codes, we generate a magnetic field with a Kolmogorov spectrum and a minimum scale 100 times smaller than the maximum scale.

In all our simulations, CRbeam and CRPropa produced consistent results. Given the fact that the generation of magnetic fields and the propagation of particles in these codes are carried out in completely different ways, this increases the reliability of the results.

For the first test, we set the magnetic field strength to $B=1$ nG, $\lambda_{C}=1$ Mpc and the proton energy to $E=10$ EeV. The radius of the final sphere has been chosen to be $D=5$ Mpc. The results are shown in the Figures 2 and 3. Upper panel of the Figure 2 shows the number of particles in each pixel of the healpy skymap. The average number of particles $N_{0}$ expected in each pixel is 1600 and the color scale is adjusted to the range of 800 - 2400 particles, which corresponds to a relative deviation of 0.5 from the mean. The lower panel shows the calculation of the integral of the rotor of the magnetic field, i.e. the pixel value corresponds to the r.h.s. of equation 4, calculated along the radius of the sphere directed towards the center of the given pixel. The color scale is also adjusted to the range -0.5 - 0.5, allowing one-to-one comparison of the graphs. It can be seen that the positions of particle density fluctuations coincide with the fluctuations of the integral of r.h.s of equation 4.

This is also clearly seen in Fig.3 where the relative amplification of the number of particles in a pixel is shown as a function of the integral 4 evaluated in the direction to that pixel. One dot on the graph represents one healpy pixel, and the black straight line corresponds to the theoretical prediction of the equation 4. In the vicinity of zero, the dots follow a theoretical linear law, however, in the region of large values of the integral, the gain becomes much stronger.

In Fig. 4 the evolution of medium-scale anisotropies of the distribution of cosmic rays on a sphere is presented as a function of distance to the source for Larmor radius $R_{L}=10$ Mpc (for $B=1$ nG and $E=10$ EeV) and coherence length $\lambda_{C}=0.3$ Mpc. The top left panel shows an almost isotropic distribution for $D=1$ Mpc from the source. At 10 Mpc, top right panel, anisotropy is strongest and sharpest, then fades away by 100 Mpc, see bottom panels.

In Fig. 5 we show how the flux averaged over the five brightest spots on the sphere changes as a function of the distance to the source for changing values of two parameters: Larmor radius $R_{L}$ and coherence length $\lambda_{C}$. In the left panel, we fix the Larmor radius equal to $R_{L}=10$ Mpc and study the dependence of the gain on distance for several coherence lengths $\lambda_{C}=0.1,0.3,1.3$ Mpc. It can be seen that with an increase in the coherence length, the gain maximum is reached at a greater distance from the source, but varies within 0.8 and 1.5 $R_{L}$. For the right panel of Fig. 5 the coherence length $\lambda_{C}=0.3$ Mpc is fixed and the distance dependence of the gain is studied for several Larmor radii $R_{L}=3,10,30,100$ Mpc. In all cases, the flux enhancement as a function of distance rises to a maximum at a distance from the source of about 1-2 Larmor radius $R_{L}$, and then drops to an average value after about 10 $R_{L}$.

In Fig. 6 we show a sky map with UHECR density at a fixed radius of the sphere around the source. This figure shows an example of a UHECR with $E=10$ EeV propagating in a turbulent magnetic field with strength $B=1$ nG and coherence scale $\lambda_{C}=1$ Mpc. With such a field and energy of cosmic rays, we expect that the maximum anisotropy in the UHECR distribution will be reached at 10 correlation lengths, i.e. at a distance of $D=10$ Mpc. This distance to the source corresponds to the lower panel. The bar under the figure shows the color correspondence to the density of cosmic rays in the structures. Three types of structures can be identified: knots with the highest density, filaments and voids with the lowest density. Note that the typical size of a magnetic field domain with $\lambda_{C}=1$ Mpc has an angular scale of 6-10 degrees on this map and corresponds to small features. However, the most prominent and highly visible are the medium-scale structures, with a typical scale of about 60 degrees. To understand the angular scale of these structures, we have shown in Fig. 6 the UHECR density at $D=10$ Mpc, but when the magnetic field is set to zero outside the 1 and 3 Mpc spheres, see top and middle panels respectively. It can be seen that most of the global structure observed at a distance of 10 Mpc is formed during the passage of the first 3 correlation lengths from the source, its physical size corresponded to the correlation length at that time. Later on the structure mainly sharpens.

It is clear that the observer stands only at one given point on the sphere. By randomly changing its location, we calculate the cumulative probability of observing a flux with an intensity higher than a given value, which is shown in the Fig. 7. The x-axis is normalized to the mean density on the sphere. We plot the cumulative probability at a distance of 1, 10, and 100 Mpc from the source, from the top panel to the bottom, respectively for $\lambda_{C}=1$ Mpc, $B=1$ nG, and $E=10$ EeV. Both at small and large distances, the probability distribution is similar to Gaussian. However, in the middle panel it is very far from Gaussian. The probability of getting a flux below average is 80%. At this distance from the source, the initial flux is likely to decrease, e.g. by a factor of 10 with a probability of about 20%. The probability of an order of magnitude gain is only 2% .

Finally, we have verified that the appearance of caustic-like structures does not depend on specific parameters of the magnetic field spectrum. We first checked that our results are robust to changes in the ratio of maximum and minimum turbulence scales. We reduced the size of the smallest eddy by an order of magnitude and made sure that the results did not change. Secondly, we replaced the turbulent magnetic field with the Kolmogorov spectrum by a field consisting of cells with a size equal to the correlation length. In each cell, the field is uniform with strength $B$, but the direction changes randomly between cells. Repeating the simulation in such a magnetic field, we saw that caustics appear in this case as well.