Simulation of the Isotropic Ultra-High Energy Photon Flux in the Solar Magnetic Field

Abstract

Both the lack of observation of ultra-high energy (UHE) photons and the limitations of the state-of-the-art methodology being applied for their identification motivate studies on alternative approaches to the relevant simulations and the related observational strategies. One such new approach is proposed in this report and it concerns new observables allowing indirect identification of UHE photons through cosmic ray phenomena composed of many spatially correlated extensive air showers or primary cosmic rays observed at one time. The study is based on simulations of interactions of UHE photons with the magnetic field of the Sun using the PRESHOWER program with some essential modifications. One of the expected results of such interactions is a generation of cosmic ray ensembles (CREs) in the form of very thin and very elongated cascades of secondary photons of energies spanning the whole cosmic ray energy spectrum. Upon entering the Earth’s atmosphere, these cascades or their parts may generate uniquely characteristic walls of spatially correlated extensive air showers, and the effect is expected also in cases when primary UHE photons are not directed towards the Earth. Particle distributions in these multi-primary UHE photon footprints are expected to have thicknesses of the order of meters and elongations reaching even hundreds of millions of kilometers, making them potentially observable with a global, multi-experiment approach, including re-exploring of the historical data, with the expected event rate exceeding the capabilities of even very large cosmic ray observatories. In this report, we introduce for the first time the methods allowing for simulating the isotropic flux of UHE photons in the Sun’s vicinity. Presented methods were verified and optimised in such a way that they would successfully model the cumulative spatial distribution of secondary photons at the top of the atmosphere. The preliminary results of simulations for the UHE photon flux of energy 100 EeV demonstrate the possibility of simulating potentially observable quantities related to CRE induced by UHE photons: densities, energy spectra and geographical orientations of secondary particles at the top of the Earth’s atmosphere. A measurement of at least one of these quantities would be equivalent to a confirmation of the existence of UHE photons, which would give an insight into fundamental physics processes at unprecedentedly high energies, far beyond the reach of man-made accelerators. On the other hand, a lack of such an observation would allow for further constraining of these fundamental processes with the physically new upper limits on UHE photon fluxes after careful analysis of the technical observation ability. The novel advantage of such an approach would lay in the purely electrodynamical character of the underlying simulations which are fully independent on extrapolations of hadronic interaction models by many orders of magnitude. Such extrapolations are necessary in the UHE photon identification methods based on the analyses of properties of individual extensive air showers presently used to determine the UHE photon upper limits.

1. Introduction

The nature of the cosmic rays of the highest energies is one of the unsolved problems in physics. Their sources and mechanisms of acceleration remain unknown to us, as they are the subject of ongoing research (see, e.g., Ref. [1] and the references therein). Therefore, a detailed analysis of physical phenomena at the highest energy regime may improve our understanding of the universe and enable us to probe modern physical theories. As far as we understand the interactions of subatomic particles described by the Standard Model, there should be the cut-off of the cosmic ray energy spectrum, due to the interactions of cosmic radiation components with cosmic microwave background radiation. At cosmic ray energies of about 5–6 $\times {10}^{19}$ eV, the collisions of protons in the cosmic ray flux with cosmic microwave background photons may cause the production of the short-lived particle ${\Delta }^{+}$, which decays into a nucleon with energy lower than the initial proton energy. This phenomenon is called the Greisen–Zatsepin–Kuzmin (GZK) limit (see [2,3]). In turn, heavier nuclei of cosmic radiation at energies of about 10${}^{20}$ eV may undergo photodisintegration and in this way also lower their energy. Therefore, the maximum distance from the source for the possible observation of the highest energy cosmic rays should be limited.

Recent experimental data [4] indicate that the cosmic ray energy spectrum is suppressed, but it is impossible to determine whether the observations are consistent with the mentioned theory. On the other hand, the nature of detected cosmic radiation particles with energies exceeding 10${}^{20}$ eV is unknown, as there are no potential sources of such particles in the Earth’s vicinity. Thus, investigation of ultra-high energy cosmic rays is significant.

The problem of the highest energy cosmic rays gave rise to new theoretical models and predictions (see, e.g., Ref. [1] and the references therein) that can be justified or falsified only by detailed analysis of processes at the ultra-high energy regime. One of these predictions is the existence of super heavy dark matter, the decay or annihilation of which may lead to the production of ultra-high energy cosmic ray flux composed mainly of photons (see Refs. [5,6,7,8]). Hence, the study of ultra-high energy photons can be considered as an indirect search for heavy dark matter particles. Similarly, the rest of the new physical theories predict ultra-high energy (UHE) photon creation. However, the expected UHE photon flux is much higher than the observed one. A possible explanation for lowering the UHE photon flux is the cosmic ray ensembles (CREs) phenomenon. Providing that a UHE photon, namely the photon of energy greater than 10${}^{18}$ eV, is cascading already in space, then a group of secondary cosmic rays might reach the top of the atmosphere instead of an individual particle. Thus, instead of only one extensive air shower (EAS) created in the atmosphere, a group of EASs is generated. This correlation may involve arrival time relations and/or spatial distributions to provide a characteristic signature, although a detection of this phenomenon might require an extended detector array, even about the size of the Earth, as proposed by the Cosmic-Ray Extremely Distributed Observatory (CREDO) [1]. The project aims at analysing data from professional observatories, educational cosmic ray facilities, and also from the smartphones with the CREDO Detector [9,10] mobile application, where a CMOS camera is used as a particle detector. For proper collected data interpretation dedicated to identification of CRE, the process of UHE photon cascading, formation of correlated air showers, and detector response need to be correctly modelled. This is the reason why, among other investigations, the CREDO project probes the scenario of UHE photon cascading in the solar magnetic field [11]. The specific methods dedicated to utilising this phenomenon for planning novel indirect UHE photon searches are presented in this article.

Previous studies of the CRE effect [1,11] show that it is possible to observe some photons from cascades generated by the primary UHE photons also in cases where they are not aimed at the Earth. Here, we use the PRESHOWER program [12,13] with modifications to include such cases in the simulations for the first time, and deliver complete distributions of CRE photons at the top of the atmosphere which point to relevant observational strategies based on the properties of these distributions.

2. Model of Ultra-High Energy Photon Interaction

We modelled the Sun’s vicinity by a surrounding sphere with radius $6{\mathrm{R}}_{\odot }$. It was also assumed that UHE photon flux reaching the starting sphere is isotropic. However, in the case of the UHE photon direction, only photons that come from outside of the sphere are considered to avoid the case duplication. Limits on the UHE photons are applied in accordance with [14]. Then, the cascading of the UHE photon in the solar magnetic field is considered.

Physical processes that are responsible for cascading the UHE photon in space are magnetic pair production, the bending of the charged particle trajectory, and synchrotron photon emission. In the considered case, a primary UHE photon may convert to a positron and electron in a sufficiently strong magnetic field. Then, trajectories of the created pair bend due to the presence of a magnetic field. Hence, the positron and the electron emit synchrotron photons that may be observable as photons cascade on the top of the atmosphere. The following mathematical description of the UHE photon cascading is based on Ref. [12] and the references therein. The mathematical formulations of quantum mechanics of magnetic pair production used in this work were taken from [15,16]. On the basis of [15], the number of pairs created by $n_{photons}$ photons in transverse magnetic field B is described by the following expression

$$n_{pairs}=n_{photons}\left[1-exp\left(-\alpha \left(\chi \right)dl\right)\right],$$

(1)

where $dl$ is the photon path length and $\alpha \left(\chi \right)$ is the attenuation coefficient that depends on parameter

$$\chi =\frac{1}{2}\frac{h\nu }{m_e{c}^2}\frac{B}{B_{cr}},$$

where $h\nu$ is the photon energy, $m_e$ is the electron mass, and $B_{cr}=4.414\times {10}^{13}$ G is the critical magnetic field strength. The attenuation coefficient depends on the probability transition for magnetic pair production. Using the first-order perturbation method and assuming that $B\ll B_{cr}$, the approximated formula for $\alpha \left(\chi \right)$ may be expressed as

$$\alpha \left(\chi \right)=\frac{1}{2}\frac{{\alpha }_{QED}}{{¯\lambda }_c}\frac{B}{B_{cr}}T\left(\chi \right),$$

where ${\alpha }_{QED}$ is the fine-structure constant, ${¯\lambda }_c$ is the reduced Compton wavelength, and an auxiliary function $T\left(\chi \right)$ is approximated by

$$T\left(\chi \right)\approx \frac{0.16}{\chi }{K}_{1/3}^2\left(\frac{2}{3\chi }\right)$$

(2)

whereas the $K_{1/3}\left(\frac{2}{3\chi }\right)$ in the expression (2) is the modified Bessel function.

According to Equation (1), the probability of conversion within an infinitesimal distance $dl$ is

$$p_{conv}=1-exp\left(-\alpha \left(\chi \right)dl\right)\approx \alpha \left(\chi \right)dl$$

(3)

and for larger distance L it should be calculated as

$$P_{conv}=1-exp\left[-{\int }_0^L\alpha \left(\chi \right)dl\right].$$

The conversion rate is negligible unless the condition $\chi =\frac{1}{2}\frac{h\nu }{m_e{c}^2}\frac{B}{B_{cr}}\approx 0.1$ is satisfied. Thus, the magnetic pair production has not been observed so far.

The parameter $\chi$ also determines the shape of the pair-member energy distribution, which is given by the expression

$$\frac{dn}{dϵ}\approx \frac{{\alpha }_{QED}}{{¯\lambda }_c}\frac{B}{B_{cr}}\frac{\sqrt{3}}{9\pi \chi }\frac{[2+ϵ(1-ϵ\left)\right]}{ϵ(1-ϵ)}{K}_{2/3}\left[\frac{1}{3\chi ϵ(1-ϵ)}\right],$$

(4)

where $ϵ$ is a fraction of the primary photon energy carried by the pair-member. The derivation of this formula is presented in detail in [16]. The radial acceleration of electron and positron in the solar magnetic field is calculated using

$$\frac{d\widehat{\mathbf{v}}\left(t\right)}{dt}=\frac{q{c}^2}{E}\widehat{\mathbf{v}}\times \mathbf{B},$$

where $\widehat{\mathbf{v}}$ is a velocity versor, q is particle charge, E is particle kinetic energy and $\mathbf{B}$ is an external magnetic field. For a short time interval $\Delta t$, altered particle direction $\widehat{\mathbf{v}}(t+\Delta t)$ may be approximated with the Taylor series expansion

$$\widehat{\mathbf{v}}(t+\Delta t)\approx \widehat{\mathbf{v}}\left(t\right)+\frac{d\widehat{\mathbf{v}}\left(t\right)}{dt}\Delta t$$

and in this case

$$\widehat{\mathbf{v}}(t+\Delta t)\approx \widehat{\mathbf{v}}\left(t\right)+\frac{q{c}^2}{E}(\widehat{\mathbf{v}}\times \mathbf{B})\Delta t.$$

(5)

In the PRESHOWER simulator, an appropriate time interval $\Delta t$ is chosen according to the particle energy and the magnetic field variability. During the first half of the $\Delta t$ period, the particle propagates in the previous direction. Then the particle direction is recalculated in accordance with (5) and in the second half of the $\Delta t$ the particle is moving in an altered direction. In the end, the particle direction is recalculated once again.

The energy distribution of the emitted synchrotron radiation for an electron at an ultra-relativistic regime, according to [17], is expressed by

$$f\left(y\right)=\frac{9\sqrt{3}}{8\pi }\frac{y}{{(1+\xi y)}^3}\left[{\int }_y^{\infty }{K}_{5/3}\left(z\right)dz+\frac{{\left(\xi y\right)}^2}{1+\xi y}{K}_{2/3}\left(y\right)\right],$$

with the parametrisation $\xi =\frac{3}{2}\frac{B_{\perp }}{B_{cr}}\frac{E}{m_e{c}^2}$, where E is the electron/positron energy, and

$$y\left(h\nu \right)=\frac{h\nu }{\xi (E-h\nu )},$$

where $h\nu$ is the emitted photon energy. In turn, the energy emitted by the electron or the positron within the path length $dl$ is equal to

$$W_{dl}=P\frac{dl}{c}=\frac{2}{3}\frac{r_0^{2}c}{{\left(m_e{c}^2\right)}^2}{E}^2{\left|\mathbf{B}\right|}^2\frac{dl}{c},$$

where P is the radiation power and $r_0$ is the classical electron radius. Therefore, one can obtain that the probability of emission of a synchrotron photon within infinitesimal path length $dl$ is

$$P_{brem}(B_{\perp },E,h\nu ,dl)=dl{\int }_0^{E}I(B_{\perp },E,h\nu )\frac{d\left(h\nu \right)}{h\nu },$$

(6)

where

$$I(B_{\perp },E,h\nu )=\frac{h\nu dN}{d\left(h\nu \right)dl}=f\left(y\right)\frac{W_{dl}}{dl}\frac{dy}{d\left(h\nu \right)}$$

and $dN$ is the number of photons with energy $[h\nu ,h\nu +d(h\nu \left)\right]$, emitted within the path length $dl$. According to [18], at an ultra-relativistic regime, synchrotron photons are emitted in a half opening angle

$$\theta \approx \frac{1}{\gamma }$$

(7)

around charge particle direction, where $\gamma$ is the Lorentz factor. The azimuthal angle of emitted photon direction is random, consistent with uniform distribution. Therefore, at the ultra-relativistic regime, the direction of the emitted synchrotron photon is nearly tangential to the electron or positron trajectory. In this way, synchrotron photons may reproduce the projection of a pair-member trajectory onto the atmosphere plane. For that reason the spatial distribution of secondary photons on the top of the atmosphere is expected to be a very thin curve.

The efficiency of the simulation method presented in this article is independent of the Sun magnetic field model. Here, for simplicity, we use a dipole approximation with a magnetic dipole moment equal to $6.87\times {10}^{32}\mathrm{G}{\mathrm{cm}}^3$, but any other could also be used, e.g., DQCS [19], as was done in Ref. [11].

3. Materials and Methods

In order to model the UHE photon cascading in the solar magnetic field, the PRESHOWER program is used in this work. In this program, interactions are modelled by Monte Carlo simulations. However, the program had to be modified to enable the time cumulative distribution of secondary photons on the top of the atmosphere. The detailed description of the PRESHOWER program is presented in [13]. Unfortunately, the article refers to the simulator of the UHE photon interaction with the geomagnetic field, so the simulator version applied in this work is different. However, this section and the indicated article provide sufficient knowledge to understand research methods.

The PRESHOWER simulator generates sample data for determined initial conditions. Simulation parameters that could be changed were the primary photon energy, the impact parameter, and the heliocentric latitude of the primary photon trajectory. Therefore, only two coordinates of the simulation starting point could be adjusted. What is more, the primary photon direction was completely determined after the choice of parameters. For these reasons, it was necessary to update the PRESHOWER simulator, with the purpose of applying the more realistic model. The introduced modifications and program functioning are described in the following sections.

The coordinate system used in the simulator is shown in Figure 1. The applied coordinate system is similar to the heliocentric ecliptic coordinate system, but in contrast to this system, the Earth’s motion is not simulated. Therefore, the x axis is determined by the centre of the Earth and the centre of the Sun. The Earth’s centre is also the coordinate system origin. The $xy$ plane contains the ecliptic. In fact, the simulator neglects the Earth and its atmosphere curvature. The atmosphere flat is modelled by a plane with equation $x=112.83$ km. In the scheme, simulation parameters–the impact parameter and the heliocentric latitude of the primary photon trajectory—are indicated. The impact parameter in this model is the distance between the primary photon trajectory and the Sun’s centre, whereas the heliocentric latitude is defined as the angular distance between the ecliptic and a simulation starting point in the heliocentric coordinate system. In this way, the simulation starting point is completely determined because the x coordinate of the starting point is a constant value for all simulations equals to $156.9\times {10}^6$ km, which is implemented in the code and cannot be adjusted.

Moreover, the primary photon is always aimed at the centre of the atmosphere plane. Therefore, it is assigned after a simulation parameter choice. In this case, the heliocentric latitude remains unchangeable throughout the primary photon propagation.

For proper functioning of the program, many functions with different destinations need to be applied, among others, the CERN routine DBSKA for calculating the modified Bessel functions. The most important function of the simulator performs a loop on particles and simulates physical processes involving them and propagates them. The loop includes only the primary UHE photon and possibly the primary pair, because secondary photons’ directions are not altered and they do not need to be propagated step-wise. Thus, their position at the top of the atmosphere is calculated just after their generation and potential secondary photon cascading is neglected.

The function calculates a simulation starting point for given input parameters. Then, it propagates the UHE photon and simulates its conversion gradually. In each step, the probability for conversion is computed according to Formula (3).

For a stochastic sampling of the probability for the process occurrence, the long period pseudo-random number generator is used. Generated numbers are consistent with the uniform distribution. Therefore, thanks to imposing the condition that the process occurs if the generated number is less than the calculated probability, the simulation reflects the probability distribution for the considered process. When the electron-positron pair is created, the fractional energy is simulated similarly to the conversion probability, with respect to its probability distribution (4).

Then, the electron and positron are propagated step-wise according to Formula (5). The distance $\Delta s$ chosen in each step, which refers to the time interval $\Delta t$, depends on the magnetic field strength for the particle location in space and the current particle energy. Secondary photon emission is simulated by analogy with the primary photon conversion. The probability for this process is computed on the basis of Formula (6). The approximation (7) is adopted as the angle of synchrotron photon deflection from the electron or positron trajectory. In turn, the angle of photon direction in the plane horizontal to the emitting particle trajectory is sampled from the uniform distribution. The simulation ends when the primary photon or the created particles reach the top of the atmosphere. All obtained particle data are stored in a matrix, which has a fixed size. Therefore, the maximum number of particles generated in the simulation is limited.

On the basis of the simulation results, the spatial distribution of secondary photons at the top of the atmosphere may be analysed in detail in the case of the UHE photon aimed at the top of the atmosphere. The influence of simulation input parameters on the distribution of secondary photons at the top of the atmosphere was the subject of a previous study [11].

Previous results indicated that it is justified to consider a more realistic case, in which UHE photons reaching the Sun’s vicinity from every direction are modelled. Furthermore, it was assumed that the UHE photon flux is isotropic. Nevertheless, this approach required substantial changes in the program’s functioning.

The most crucial modification concerned the arbitrary direction of the UHE photon, which implies the need for new simulation parameters and a new method of calculating the simulation starting point. On the basis of assumptions made about the UHE photon flux, it was found that the most suitable for simulations is the heliocentric coordinate system. Figure 2 illustrates the modified coordinate system. It is worth mentioning that the inclination angle $\theta$ in the coordinate system is defined as an angle between the x axis and the position vector.

In the new version of PRESHOWER, simulation starting points should be isotropically distributed at the sphere surface. Therefore, the radius of this starting sphere was chosen as the simulation parameter, whereas the position on the starting sphere is generated in such a way that every solid angle is equally probable. This condition may be described as

$$\frac{dN}{d\Omega }=\mathrm{const}\equiv C,$$

(8)

where $dN$ is the number of UHE photons at the starting sphere in the solid angle $d\Omega$. Therefore, the probability of generating the primary photon at the position with the azimuthal angle from interval $[\phi ,\phi +d\phi ]$ should be equal to

$$P_{\phi }=\frac{1}{N}{\int }_{{\Omega }_0}\frac{dN}{d\Omega }d\Omega =\frac{1}{N}{\int }_0^{\pi }sin\theta d\theta {\int }_{\phi }^{\phi +d\phi }Cd\phi =\frac{C}{N}d\phi ,$$

where N is the total number of UHE photons reaching the starting sphere. The probability density function of the azimuthal angles is constant

$$f\left(\phi \right)=\frac{d{P}_{\phi }}{d\phi }=\frac{C}{N}=\mathrm{const}$$

(9)

and the variable is described by the uniform distribution. Similarly, one may obtain that

$$f(cos\theta )=\mathrm{const}$$

(10)

and $cos\theta$ is uniformly distributed. Therefore, in the applied function $\phi$ and $cos\theta$ are randomly generated from uniform distribution and then the simulation starting point is calculated using the standard transformation between spherical coordinates and Cartesian coordinates. In the transformation, the redefinitions of the $\theta$ angle and the Sun’s position are included. The isotropic primary photon direction is modelled analogically. Nevertheless, the transformations from $(\tilde{\theta },\tilde{\phi })$ coordinates to coordinates used in the simulation are more complex. Due to the fact that $(\tilde{\theta },\tilde{\phi })$ are defined in the local coordinate system related to the position on the starting sphere, generated angles are firstly transformed to the local Cartesian system

$$\begin{array}{c}{\widehat{e}}_{\tilde{x}}=-sin\theta cos\phi ,{\widehat{e}}_{\tilde{y}}=-sin\theta sin\phi ,{\widehat{e}}_{\tilde{z}}=-cos\theta .\end{array}$$

The minus sign results from considering primary photons coming from the outside of the starting sphere. The local versors may be identified as spherical versors

$${\widehat{e}}_{\tilde{z}}={\widehat{e}}_r,{\widehat{e}}_{\tilde{x}}={\widehat{e}}_{\theta },{\widehat{e}}_{\tilde{y}}={\widehat{e}}_{\varphi }$$

(11)

related to the position on the starting sphere $({\theta }_S,{\phi }_S)$. Therefore, they are expressed in coordinates used in the simulator (Figure 1) using formulas

$$\begin{array}{c}{\widehat{e}}_x={\widehat{e}}_{\tilde{z}}cos{\theta }_S-{\widehat{e}}_{\tilde{x}}sin{\theta }_S\\ {\widehat{e}}_y={\widehat{e}}_{\tilde{z}}sin{\theta }_{S}cos{\phi }_S-{\widehat{e}}_{\tilde{y}}sin{\phi }_S+{\widehat{e}}_{\tilde{x}}cos{\theta }_{S}cos{\phi }_S\\ {\widehat{e}}_z={\widehat{e}}_{\tilde{y}}cos{\phi }_S+{\widehat{e}}_{\tilde{x}}cos{\theta }_{S}sin{\phi }_S+{\widehat{e}}_{\tilde{z}}sin{\theta }_{S}sin{\phi }_S.\end{array}$$

Additionally, the drawn value $cos\tilde{\theta }$ belongs to the interval $[0,1]$ since we consider only the range of polar angles $[0,\pi /2]$, which responds to directions beneath the sphere surface.

The generalisation of the program for an arbitrary primary photon direction required also changes in the simulation end condition. The previous condition may be effectively applied only to primary photons aimed at the atmosphere plane. Thus, assuming the plane is infinite, the previous condition is placed on photons with the negative component of the direction versor. In other cases, the simulation end condition is the distance to the Sun’s centre larger than 7 $R_{\odot }$. What is more, the propagation of the primary photon, electron, or positron is aborted in case the particle lands on the Sun. In turn, secondary photons aimed at the Sun are not even included in the particle table. In the previous version of the program, these additional conditions were not needed, because of the specific determination of the primary photon direction. This approach enables considering also the potential special cases of the UHE photon cascading. For instance, it takes into account the primary photons not aimed at the Earth that may generate observable secondary photons after conversion thanks to the change of electron or positron direction in the magnetic field.

4. Results

Introduced amendments allowed simulation of the cumulative distribution of secondary photons at the top of the atmosphere. However, the new version of the simulator needed to be verified in detail beforehand. Therefore, both the validation of the simulator and the preliminary analysis of the results are described in the following sections.

4.1. Verification of the PRESHOWER Functioning

It was significant to analyse simulation results firstly with respect to the correctness. This verification was performed both for appended functions and for the simulation as a whole. In the case of added functions, we visualised their results separately and compared them with expected outcomes. For this purpose, we analysed the distribution of solid angles generated without performing complete simulations and then visualised calculated positions at the sphere and directions. On the basis of 100,000 generated $cos\theta$ and $\phi$ values, it was found that the distribution of solid angles is flat. Similarly, the visualisations of position and direction distributions were in line with our expectations. With the purpose of verifying the whole program and the interoperability of its functions, it was necessary to compare results from the present version of the PRESHOWER with those from the previous version. However, this comparison should concern an identical physical event. Therefore, the simulation starting point, the primary photon direction, and the energy must be the same in both used simulator versions. We considered the case when the primary photon with energy 100 EeV aimed at the Earth propagates in the Sun’s vicinity with the impact parameter 3 $R_{\odot }$ and the heliocentric latitude 0${}^{\circ }$. This responds to the following parameters in the current PRESHOWER version: the simulation starting point $(R,{\theta }_S,{\varphi }_S)=(10.44R_{\odot },0.29,0)$, where R is the starting sphere radius, and the primary photon direction $(\tilde{\theta },\tilde{\varphi })=(0.28,\pi )$. In this simple case, the primary photon, electron, and positron are propagating only in the $xy$ plane. Their movement modelled by both simulations is illustrated in Figure 3. The graph shows particle positions for each thousandth step of the simulation.

It was found that the propagation of particles is similar in both cases. Minor differences such as distances between the next steps are caused by a different conversion point. This, in turn, results from the stochastic nature of physical processes. Furthermore, the distribution of secondary photons at the top of the atmosphere is presented in Figure 4. The generated distribution has a line-like shape and a size $8.48\times {10}^4$ km, which is a reasonable outcome. In accordance with previous results [11], the distribution size for chosen input parameters should be approximately ${10}^5$ km.

Therefore, obtained results are correct and it was concluded that the simulator properly models UHE photon cascading.

4.2. Optimisation of the Simulation

The inclusion of more cases in the simulation, unfortunately, extended the computation time. It was found that the simulator needs to be optimised. Due to this fact, we reduced the cases included by the program and adjusted its parameters.

4.2.1. Reducing Considered Cases

The first simulations showed that only a few of them gave rise to observable effects on the Earth. Thus, most of the calculation time was used for cases that are not the subject of this study. With the purpose of effective use of computer resources, we refrained from modelling the special cases of the UHE photon cascading. Additional conditions were placed on the primary photon direction and its arrival point at the top of the atmosphere. Only primary photons moving towards the negative x axis are simulated. In turn, the primary photon arrival point is usually related to the centre of secondary photons distribution. In this way, limiting its distance from the x axis (Figure 2) increases the probability of the observation on the Earth. The limit on the distance between the point at the top of the Earth’s atmosphere where the UHE photon would have landed if no conversion had occurred, and the x axis connecting the Earth’s and the Sun’s centres is hereinafter referred to as the arrival point limit. Conditions are checked at the beginning of the program and further procedures are performed only for primary photons that fulfil limits. What is more, the number of secondary photons saved to the outfiles was limited to particles that are not further than 6500 km from the x axis at the top of the atmosphere.

4.2.2. Adjustment of Simulation Parameters

The implemented limit on the centre of the secondary photons spatial distribution and the starting sphere radius—the simulation input parameter—should also be adjusted for the optimisation of the simulation. The limit on the primary photon arrival point has an influence on the ratio of events observable on the Earth to all simulations that fulfilled the condition. The number of simulations giving observable effects increases with the decrease in the distance from the x axis to the centre of the secondary photons signature. On the other hand, the limit obviously affects the fraction of simulations that fulfil the condition and in this way the computation time. The strengthening of the limit on the primary photon arrival point results in a decline in the ratio of simulation fulfilling the condition for all running simulations and considerably prolongs the calculation time. These relations are shown in Figure 5 and Figure 6.

On the basis of the obtained data, we calculated the ratio of events observable on the Earth to all running simulations, which is illustrated in Figure 7. It was observed that the fraction of detectable events is greater for the relaxation of the condition. However, the statistics for this range are based on few entries, so it may be burdened with high risk.

Figure 5 shows that the number of simulations fulfilling the condition depends weakly on the starting sphere radius. Thus, the fraction of simulations observed on the Earth was calculated only for the starting sphere radius 6 $R_{\odot }$ (Figure 6). In order to consider the largest number of events observable on the Earth in a relatively short time, the arrival point limit was set at ${10}^4$ km. The number of CRE effects in all effects detectable on the Earth, in turn, is dependent on the starting sphere radius. Due to the fact that the starting sphere radius has an influence on impact parameters of primary photons, it also affects the probability of their conversion. On the other hand, the decrease in the starting sphere radius may also eliminate some cases that can give observable effects at the top of the atmosphere. Therefore the starting sphere radius $6R_{\odot }$ was chosen as an optimum value for the simulation parameter. The resulting preliminary cumulative distribution of secondary photons at the top of the atmosphere is shown in Figure 8. The results were obtained only for the primary photons of energy 100 EeV. However, the presented method may be applied to the whole energy spectrum of UHE photons.

We enlarged the considered observable area to take into account the cross section of the Earth and the surrounding atmosphere. The distribution includes 324 CRE effect signatures, which corresponds to one year of observation taking into account present limits on the UHE photon flux [14]. However, the presented distribution does not include cases reduced due to the calculation time and technical limitations of detectors. Hence, the number of predicted CRE effect signatures per year may be fraught with great uncertainty. Observed line-like signatures are parts of large cascades with centres beyond the Earth’s cross-section. Due to the fact that cascades are produced by UHE photons in the vicinity of the Sun, many secondary photons have directions contained in the solid angle determined by the Sun. Therefore, CRE effect cascades should distinguish themselves from the rest of the observations since secondary photons may reach higher energies than photons emitted by the Sun. The anomaly in the very-high-energy photon flux was already observed (see [20,21]). However, the UHE photon cascading was not taken into account as a possible source to date. Therefore, it is significant to analyse the cumulative distribution of secondary photons with regard to unique features that could help to classify them. The possible anisotropy may be observed in the distribution of the orientation of signatures at the top of the atmosphere. The horizontal or nearly horizontal orientation could be more frequent. It would also be theoretically motivated since signatures oriented in this way are generated by UHE photons that have a heliocentric latitude around $0^{\circ }$ or ${90}^{\circ }$ at the conversion point [1]. Simultaneously, the conversion occurs more often for the heliocentric latitude ${90}^{\circ }$ [11]. In order to describe the orientation of signatures, the rotation angle parameter was introduced. It is given by the equation

$$\alpha =arctan\left(\frac{z_1-z_2}{y_1-y_2}\right),$$

where points $(y_1,z_1)$ and $(y_2,z_2)$ are related to secondary photons that are the farthest, but still located in the observable area. The rotation angle is in the range of $(-\pi /2,\pi /2)$. The obtained distribution is shown in Figure 9.

5. Discussion

The previous study ([1,11]) showed that the UHE photon interaction with the solar magnetic field may result in the extremely extended, line-like shaped cascade of secondary photons reaching the top of the Earth’s atmosphere. Owing to the unique signature, the footprint of secondary photons at the top of the atmosphere should be possible to observe and classify. Due to the expected wide spatial span of these cascades, it should be possible to observe not only secondary photons produced by primary UHE photons aimed at the Earth, but also the tails of the secondary photons produced by primary UHE photons directed far from the Earth. In this article, we presented a successful method for modelling the interactions of the isotropic UHE photon flux with the solar magnetic field in order to obtain for the first time particle distributions containing secondary photons from both central and peripheral parts of the CRE cascades. The method is based on using the PRESHOWER program appropriately modified and optimised. We verified the functioning of the simulator and checked that obtained results are perfectly consistent with the previous analyses [11]. In order to investigate the cumulative distribution of secondary photons at the top of the atmosphere, we discussed the impact of simulation parameters on the modelled physical phenomena and found optimal values for these parameters. The results obtained in this way indicate that secondary photons produced by the UHE photons not aimed at the Earth contribute significantly to the cumulative distribution of secondary photons at the top of the atmosphere and to the distributions of CRE footprint orientations which might serve as a valuable observational hint. In relation to the methods described here, it is important to note that recent observations ([20,21]) demonstrate yet unexplained gamma-ray emission from the solar disk. In light of the simulations of CRE resulting from the interactions of UHE photons nearby the Sun, one might speculate that these observations may be a signature of the CRE effect, which motivates a dedicated analysis in this direction. We conclude that it is justified and indispensable to simulate the flux of UHE photons cascading in the Sun’s vicinity without neglecting the tails of the resultant CRE distributions, and that the presented method enables us to search for potentially observable signatures of the CRE effect, and hence offer a novel perspective for an indirect search for UHE photons. We proposed that one such signature might be seen in a distribution of the spatial orientations of the CRE cascades at the top of the atmosphere. What is more, the prepared method also allows us to estimate the flux of CRE cascades per year and quantise its increase, which could be directly compared to the anomalous gamma-ray emission from the solar disk. A quantitative evaluation of the appropriateness of proposed observables for the CRE and UHE photon searches would require a separate study, which has already been planned. Furthermore, a modified PRESHOWER program may be used along with standard simulators of the cosmic rays’ propagation to model the UHE photon conversion in the vicinity of other distant astrophysical objects. Therefore, it may also provide a new method for investigating the CRE cascades of galactic or extragalactic origin.