Testing effects of Lorentz invariance violation in the propagation of astroparticles with the Pierre Auger Observatory

Lorentz invariance violation (LIV) is often described by dispersion relations of the form $E_i^2=m_i^2+p_i^2+\delta_{i,n} E^{2+n}$ with delta different based on particle type $i$, with energy $E$, momentum $p$ and rest mass $m$. Kinematics and energy thresholds of interactions are modified once the LIV terms become comparable to the squared masses of the particles involved. Thus, the strongest constraints on the LIV coefficients $\delta_{i,n}$ tend to come from the highest energies. At sufficiently high energies, photons produced by cosmic ray interactions as they propagate through the Universe could be subluminal and unattenuated over cosmological distances. Cosmic ray interactions can also be modified and lead to detectable fingerprints in the energy spectrum and mass composition observed on Earth. The data collected at the Pierre Auger Observatory are therefore possibly sensitive to both the electromagnetic and hadronic sectors of LIV. In this article, we explore these two sectors by comparing the energy spectrum and the composition of cosmic rays and the upper limits on the photon flux from the Pierre Auger Observatory with simulations including LIV. Constraints on LIV parameters depend strongly on the mass composition of cosmic rays at the highest energies. For the electromagnetic sector, while no constraints can be obtained in the absence of protons beyond $10^{19}$ eV, we obtain $\delta_{\gamma,0}>-10^{-21}$, $\delta_{\gamma,1}>-10^{-40}$ eV$^{-1}$ and $\delta_{\gamma,2}>-10^{-58}$ eV$^{-2}$ in the case of a subdominant proton component up to $10^{20}$ eV. For the hadronic sector, we study the best description of the data as a function of LIV coefficients and we derive constraints in the hadronic sector such as $\delta_{\mathrm{had},0}<10^{-19}$, $\delta_{\mathrm{had},1}<10^{-38}$ eV$^{-1}$ and $\delta_{\mathrm{had},2}<10^{-57}$ eV$^{-2}$ at 5$\sigma$ CL.


Introduction
According to the modern physics description of nature, Lorentz invariance (LI) is a fundamental symmetry. Several tests of LI have been performed and no evidence for symmetry breaking has been found. Nevertheless, the possibility of some level of Lorentz invariance violation (LIV) at very-high energies has been proposed in the context of theories beyond the Standard Model, quantum gravity and string theory [1].
Ultra-high energy cosmic rays (UHECR, E > 10 18 eV) are thus expected to be a sensitive probe of possible LIV [2]. Being the highest energy particles ever detected and dominantly coming from extra-galactic sources at energies E 10 19 eV [3], they represent a unique opportunity to test LIV beyond the energies available in man-made accelerators. UHE messengers interact with the background photon fields and the particles reaching Earth keep information about the main features of the interactions. Relativistic kinematics is changed under LIV assumptions and, therefore, the data measured on Earth might contain imprints of LIV.
In this work, we search for LIV signatures in the data collected at the Pierre Auger Observatory, while adopting realistic scenarios for UHECR characteristics. Two independent LIV sectors are tested. The pair production of UHE photons is simulated in the extragalactic space under LIV hypotheses and the flux of secondary photons is calculated accordingly. With that, our upper limits on the UHE photon flux are used to test LIV in the electromagnetic sector. Independently, the measured energy spectrum and distributions of depth of shower maximum, X max , are used together, for the first time, to set limits in the hadronic sector, as proposed in [4][5][6]. Dedicated simulations considering photopion production and photodisintegration of nuclei in the extragalactic propagation under LIV assumptions are compared to the data. In this work we do not take into account the possible effects of modifications of the development of the atmospheric showers due to LIV.
In section 2, the phenomenological approach of the LIV hypothesis used throughout this paper is explained. A description of the Pierre Auger Observatory and of the data used here are presented in section 3. The tested UHECR scenarios are reported in section 4. The results for the electromagnetic and hadronic sectors are presented independently in sections 5 and 6, respectively. Discussion and conclusions are finally addressed in section 7.

Lorentz invariance violation framework
A commonly used phenomenological approach for LIV was firstly developed by Coleman and Glashow [7,8]. It reduces the LIV effects to a change in the energy dispersion relation which can be expressed as: where i identifies the particle, E, m and p denote, respectively, the energy, mass and momentum of the particle, and s can be chosen to be either the energy, E, or the momentum, p, this choice being irrelevant for ultra-relativistic particles. We consider s = E in the following. The extra term in the energy dispersion relation thus modifies the propagation and kinematics of each particle type. The experimental fact that LI is preserved at low energies, or in other words that LIV could only be significant at high energies, can be expressed by taking as function of E and expanding it in a polynomial series, so that where n is the order of the perturbation. The parameters δ i,n are in principle independent for each particle species and define the energy scale associated with the violation. We note that because LIV is often associated with quantum gravity theories, an expansion in terms of the Planck scale is sometimes introduced as where M Pl ≈ 1.22 × 10 19 GeV/c 2 is the Planck mass. If we assume that only the lowest-order non-vanishing term has a non-negligible effect, it is possible to set individual limits on δ i,n .

The Pierre Auger Observatory and the datasets
The Pierre Auger Observatory [9] is the largest UHECR detector currently in operation. Located in the southern hemisphere in western Argentina, just northeast of the town of Malargüe (69 • W, 35 • S, 1400 m a.s.l.), it covers an area of 3000 km 2 with a Surface Detector array (SD) [10] overlooked by a Fluorescence Detector (FD) [11]. The SD consists of 1660 water-Cherenkov detectors arranged in a triangular grid operating with a nearly 100% duty cycle. Each SD station detects at ground level the secondary particles of the extensive air shower (EAS) produced by the primary UHECR interacting in the atmosphere. The FD consists of a set of telescopes that measure the UV fluorescence light from nitrogen molecules excited by the EAS particles along their path in the atmosphere. FD operations are limited to clear moonless nights, resulting in a duty cycle of ∼15% [12].
The combination of the information from both techniques results in a quasi-calorimetric determination of the energy scale, a geometric direction reconstruction and an estimator of the primary particle mass. Details of the reconstruction techniques and their efficiency can be found in [9]. The determination of the features of the energy spectrum measured by the Pierre Auger Observatory reached unprecedented precision [13,14], achieved with a continuous and stable operation of the detectors combined with the hybrid energy calibration. In this analysis, we use the energy spectrum measured using an accumulated exposure of 60,400 ± 1,800 km 2 sr yr, obtained by the combination of complementary data sets above 3 × 10 18 eV based on 215,030 events recorded with zenith angles below 60 • [15]. The events were measured by the SD between 1 January 2004 and 31 August 2018. A subset events with energy above 3 × 10 18 eV were collected between 1 January 2004 and 31 December 2017 simultaneously by the SD and FD detectors. This subset is used to calibrate the energy reconstruction of the events.
The atmospheric depth at which the air shower peaks in the rate of energy deposit, X max , is the most accurate available parameter with proven correlation with the mass of the primary particle. The FD measures the fluorescence light produced in the development of the shower in the atmosphere from which the X max can be reconstructed. In this paper, we use X max distributions of the events measured from 1 December 2004 to 31 December 2017 with energy above 10 17.8 eV and zenith angle below 65 • [16]. The 35425 events surviving the analysis and selection procedure [17] are sampled into the X max distribution for each of 19 energy bins [10 17.8 , 10 17.9 ) eV, . . ., [10 19.5 , 10 19.6 ) eV, and [10 19.6 eV,+∞).
No photon has yet been conclusively detected by the Pierre Auger Observatory, leading to the most restrictive upper limits on the photon flux at the highest energies [18][19][20]. Multivariate analyses based on SD and FD data was perfomed in the all events with zenith angle between 30 • and 60 • measured from 1 January 2004 to 30 June 2018. Since no photon-like event was found with an exposure of 40,000 km 2 sr yr, the null result was transformed into flux upper limits for energies above 10 18 eV

The UHECR scenario
In the present work, we test the sensitivity of the data collected by the Pierre Auger Observatory to violations of LI under specified UHECR scenarios. We consider isotropically distributed sources with (1 + z) m cosmological evolution emitting an energy spectrum given by where z is the redshift, Γ is the spectral index at the injection, R cut is the cutoff rigidity, f A is the fraction of nuclei with mass A, and J 0 is the normalization factor of the flux which enters into the computation of the total emissivity defined as dE is the number of nuclei with mass A injected per unit energy, volume and time.
UHECRs interact with photon backgrounds on the way to Earth. At these energies, the most important photon backgrounds are the cosmic microwave background (CMB) and the extra-galactic background light (EBL). Given the uncertainties in the models, we included two EBL distributions and two photo-nuclear cross section models among the UHECR propagation models used in this study, as used also in [21] and described in [22]. The comparison of events arriving at Earth with X max data is only possible if the arriving mass composition is transformed into an X max distribution with the use of EAS simulations in which the most relevant alternative is the hadronic interaction model. We used Gumbel functions [23] to parametrize the X max distribution. Two hadronic interaction models were considered.
Finally, a UHECR scenario is here defined by: a) the set of parameters from the source: Γ, R cut , J 0 , f A and m; b) the propagation model: the EBL and photo-nuclear cross sections and; c) the hadronic interaction model used to convert the arriving composition into an X max distribution. Table 1

LIV limits on the electromagnetic sector
In this section, we investigate LIV signatures in the electromagnetic sector. The interaction of UHECRs with the background radiation produces pions, as first studied by Greisen [32], Zatzepin and Kuzmin [33] (GZK). The neutral pions decay into UHE photons shortly after, which we refer to as GZK photons and are absorbed during propagation due to electronpositron pair production with background photon fields [34]. As explained below, we simulate GZK photons within a given UHECR scenario under LIV and LI assumptions and compare the predicted flux arriving at Earth with the upper limits on the photons flux measured by the Pierre Auger Observatory.

Simulation of the propagation of GZK photons under LIV assumption
If LIV is present in the photon sector, the kinematics of pair production is changed. Following the calculations from [35], we obtain the mean free path for this interaction considering subluminal LIV (δ γ < 0). Figure 1 shows the ratio between the LIV and LI mean free paths for different LIV coefficients corresponding to the LIV order n = 0. The main result is a significant increase of the mean free path above a critical energy that depends on δ γ,0 . Comparable results are found for n = 1 and n = 2. As a consequence, it is expected that if subluminal LIV is present, fewer UHE photons will be absorbed and, consequently, the arriving flux of these particles will be larger. As first proposed in [36,37] and later developed in [35], for some UHECR scenarios the predicted flux is larger than the upper limits imposed by the Pierre Auger Observatory, resulting in limits on the LIV coefficients. If superluminal LIV (δ γ > 0) were considered, the mean free path would decrease, leading to a weaker flux arriving on Earth. Therefore, the upper limits on flux are insensitive to superluminal LIV in the photon sector. The LIV-modified mean free path was implemented in the software packages CRPropa3/EleCa [38,39] in order to obtain the arriving GZK photon flux. We simulated UHECRs with energies ranging from 10 18.7 eV to 10 21 eV for five representative primaries (H, He, N, Si and Fe).
Sources were considered homogenously distributed (m = 0) going up to 9500 Mpc. The energy spectrum of UHECRs arriving on Earth was normalized at E = 10 18.75 eV. Different LIV coefficients with log 10 steps of 1 were considered for each order, n. The two limiting scenarios were also considered, LI (δ had = 0) and the case in which no interactions happen at any energy, here noted as δ had → −∞.

Results with the reference UHECR scenario
First, we consider the reference SPGE propagation model and the source parameters which best describes the energy spectrum and X max distributions measured by the Pierre Auger Observatory under LI assumption for that scenario: 08% and m = 0, as found in [21]. Using this UHECR scenario, we simulated the GZK photons arriving at Earth under LIV assumption. Figure 2 shows the simulated integral fluxes of GZK photons within the SPGE propagation model and LIV order n = 0 for several LIV coefficients (δ γ,0 ), including the two limiting scenarios. Even for the most extreme LIV case (δ γ,0 → −∞), the simulated GZK photon flux is smaller than the measured upper limits on the photon flux and, thus, this data is insensitive to subluminal LIV in the electromagnetic sector. The same conclusion is found for all combinations of propagation and hadronic interaction models in Table 1 and for LIV orders n = 1 and n = 2. Considering positive source evolution (m > 0) would slightly increase the simulated GZK photon flux but not enough to reach the measured upper limits on the photon flux.

Results with an alternative UHECR scenario
As shown in [35], under LIV hypothesis, the flux of GZK photons is expected to be more intense if a non-negligible proton component is present. In Ref. [40] it has been shown that the energy spectrum and composition data measured by the Pierre Auger Observatory can be described by models including a subdominant proton component for energies above 10 19 eV. Another possibility which could contribute to the increase of GZK photons, would be to consider an additional proton component in the energy range below the ankle, together with an additional Galactic or extragalactic component, as done in [41]. We consider only in this subsection an alternative scenario with the proton component above 10 19 eV, from Ref. [40], that best describes the energy spectrum and X max distributions measured by the Pierre Auger Observatory. The integral fraction of protons for E > 10 18.75 eV for this scenario is 31%. For more details regarding the ejected spectra, refer to Ref. [40]. The Gilmore EBL model and Talys cross-section were used. The EPOS-LHC hadronic interaction model was used to convert the arriving composition into X max distributions. Using this alternative scenario, we simulated the propagation of UHECR and GZK photons in order to obtain the flux of GZK photons arriving at Earth under a LIV assumption.

LIV limits on the hadronic sector
In this section, we investigate LIV signatures in the hadronic sector. During propagation, UHECRs interact with the background radiation and the kinematics of the dominant energy losses, photopion production and photodisintegration, are modified under LIV assumptions. We consider only positive LIV coefficients (δ had > 0). Less significant changes are found if negative coefficients (δ had < 0) are considered as the photopion production and the photodisintegration cross-sections prevent a significant decrease in the attenuation length, even if the energy threshold decreases significantly.
In this section, instead of assuming an a priori scenario, we fit the source parameters and the LIV coefficients to the energy spectrum and composition from the Pierre Auger Observatory. The latest works [42][43][44] considering effects of LIV in photopion production and photodisintegration with the aim of fitting UHECR data consider a simple pure-proton composition which is largely disfavored by the Pierre Auger Observatory data. Indeed, current data from the Pierre Auger Observatory exclude at 6.5σ a pure proton or a mixed proton and helium composition at the spectral ankle region (3−5) × 10 18 eV and is compatible with a near zero proton fraction just above 10 19 eV with a possible resurgence limited to a maximum of about 10% of the total flux above 10 19.3 eV [45]. As defined below, we allow the nuclei fractions, f A , to vary in the fitting procedure as well as the LIV coefficients. The different combinations of propagation and hadronic interaction models from Table 1 are used.
We define the LIV coefficient for hadrons in such way that δ had := δ p = δ π /2 = A n δ A , where δ p , δ π and δ A are the coefficients for protons, pions and nuclei of atomic mass A, respectively and n is the LIV order. As discussed in [43], the effects of LIV on photopion production depend on δ π − δ p , which justifies the use of the relation δ p = δ π /2 in order to have an effect of the same order in each of the interactions. The relation between δ p and δ A , on the other hand, is supported by the superposition model These dispersion relations will be used to implement the modifications to the kinematics of the interactions in simulations including LIV.

Simulation of the propagation of UHECR under LIV assumptions
The dominant energy loss for protons at the highest energies is the photopion production (p + γ → p + π) for which the attenuation length under LIV assumption was calculated based on [46]. Figure 4 shows the attenuation length as a function of energy for different values of δ had,0 . LI is represented by δ had,0 = 0. Above a critical energy, which depends on the value of δ had,0 , the attenuation length is significantly larger for LIV than for the LI assumption. Above this critical energy, if δ had,0 > 0 (LIV) the number of interactions during propagation is reduced and, consequently, the cosmic rays can travel farther than they would do under the LI hypothesis.
For nuclei, on the other hand, the leading energy loss is the photodisintegration ( A N + γ → (p or n)+ A−1 N ) with the emission of one nucleon. Multinucleon emission is less relevant and therefore not taken into account together with the LIV modifications in the propagation. Figure 5 shows the photon energy threshold in the nuclei reference frame for different values of δ had,0 and for two different nuclei (helium and iron), obtained following the calculations from [46]. Under the LIV assumption the energy threshold increases abruptly above an energy that depends on δ had,0 . Comparable results are found for the other nuclei and LIV orders, n. In a similar way to the photopion production, under the LIV assumption the UHECRs would interact less and, thus, travel farther than under the LI hypothesis.
Both the LIV modified attenuation length for pion production and the LIV modified energy threshold for photodisintegration were implemented in SimProp v2.4 [31].

Results of the combined fit considering LIV
The simulated energy spectrum and composition arriving at Earth are used to fit to the data of the Pierre Auger Observatory for energies above 10 18.7 eV, using Eq. 4.1 for weighing the UHECR spectrum at the source. The fit procedure follows the explanation in [21]. Within each UHECR scenario, the free parameters of the fit are: a) the nuclei fractions, f A , b) the index of the energy spectrum, Γ, c) the maximum rigidity, R cut , d) the normalization factor of the flux, J 0 and e) the LIV coefficient, δ had,0 . The cosmological evolution of the sources was fixed to m = 0. For each value of δ had,0 ranging from 10 −24 to 10 −18 in log 10 steps of 1 and for the LI case (δ had,0 = 0) a log-likelihood fit was done searching for the combination of the parameters which best describes the data. Figure 6 shows the evolution of the best fit parameters as a function of δ had,0 . Figure 7 shows the total deviance, D, and the deviance relative (D LIV − D LI ) to the LI case. Details on how the deviance is defined and calculated can be found in [21]. The deviance becomes too large for strong LIV (large values of δ had,0 ), and, thus, limits on δ had,0 can be imposed with a confidence level (CL) given by σ = √ D LIV − D LI . For all UHECR scenarios composed of options shown in Table 1, the data imply that δ had,0 < 10 −19 at 5σ confidence level.
Since the LIV effects are dominated by the most energetic particles, the higher order LIV coefficients can be estimated as with Z = 1 taken as the most conservative value. Considering δ had,0 < 10 −19 we have δ had,1 < 10 −38 eV −1 and δ had,2 < 10 −57 eV −2 at 5σ confidence level for all UHECR scenarios shown in Table 1. These limits are the first obtained with a fit of the spectrum and composition data including UHECR nuclei. Table 2 shows the value of δ had,0 for which we obtained the minimum deviance for each UHECR scenario. Taking this value (δ had,0 = 10 −21 ) for the STGE model, we show in figure 8 the energy spectrum and the first two moments of the X max distributions in comparison with these corresponding to the LI case. The comparison between the two panels with LI (left panel) and LIV (right panel) shows the correlation between mass composition and LIV in the fit. Lorentz violation suppresses the interactions during propagation which is compensated with a lighter composition at the sources in order to obtain the same composition on Earth.  Table 2. LIV coefficients δ had,0 at minimum deviance and the corresponding confidence level σ = √ D LIV − D LI with respect to the LI case. The model we fit to the data does not describe satisfactorily the measured energy spectrum and X max distributions, thus limiting the interpretation of the minimum deviance as an imprint of LIV in the data as discussed in Section 7.
The effect of the systematic uncertainty in the energy determination of showers measured by the Pierre Auger Observatory was considered. We chose the SPGE scenario for this study and defined the SPGE+14% and SPGE-14% scenarios representing the energy shift ±14% within the systematic uncertainties. The entire analysis procedure was repeated considering the energy shift. We verified that the systematic uncertainty of the energy scale has an impact of less than 10% on the CL of the limits we imposed on the LIV coefficient, δ had,0 . We have verified that a positive shift of the energy scale would worsen the quality of the fit. This is due to the fact that a positive shift would impose larger rigidity cutoffs or heavier masses at the source, while at the increase of the LIV parameter the interactions are suppressed and thus the preferred scenarios predict lighter masses at the sources. On the other side, corresponding to a negative shift, the same value of the LIV coefficient for which the best description of the data is found does not change.

Discussion and Conclusions
In this work, we have explored the use of the data from the Pierre Auger Observatory to test LIV effects. Two independent LIV sectors were explored. The propagation of GZK photons considering subluminal LIV was used to test LIV in the electromagnetic sector by comparison to the measured upper limits on the photon flux. The propagation of UHECR nuclei considering LIV in photopion production and photodesintegration was used to test LIV in the hadronic sector by fitting the measured energy spectrum and X max distributions.
For the electromagnetic sector, we simulated the flux of GZK photons considering LIV in the propagation within several reference and one alternative propagation model. For the reference UHECR scenarios, including the SPGE propagation model which was found to best describe the energy spectrum and X max distribution from the Pierre Auger Observatory under the LI assumption in [21], the upper limits on the photons flux imposed by the Pierre Auger Observatory are not constraining enough to set limits on the LIV coefficients. For an alternative scenario with subdominant proton component, limits of δ γ,0 > −10 −21 , δ γ,1 > −10 −40 eV −1 and δ γ,2 > −10 −58 eV −2 are imposed. The limits for such an alternative scenario are several orders of magnitude more constraining than the current most-restrictive limits imposed for subluminal LIV using TeV gamma-ray data [47] although the later limits are less dependent on the astrophysical scenario (see [48] for a compilation of current limits in the photon sector obtained with astrophysics).
In the hadronic sector, the measured energy spectrum and X max distribution were fitted under a LIV assumption for the first time. Limits at δ had,0 < 10 −19 , δ had,1 < 10 −38 eV −1 and δ had,2 < 10 −57 eV −2 were imposed at 5σ CL. These limits are valid for all the UHECR scenarios studied here. The model we fit to the data does not describe satisfactorily the measured energy spectrum and X max distributions. This can be seen by the large total deviance shown in Fig. 7 with values around 300 for 110 data points and 9 free parameters leading to a reduced deviance of about 300/101, thus limiting the interpretation of the minimum deviance between δ had,0 = 10 −22 and δ had,0 = 10 −20 (see Fig. 7) as an imprint of LIV in the data. The results of the fit showed that the data can be better described by a model with fewer interactions during propagation, but the range of possible UHECR scenarios was not explored as completely as it would be needed for a claim that the data is better described by LIV instead of LI.
In both electromagnetic and hadronic sectors, the predication of the LIV imprints is heavily dependent on the composition of the cosmic rays leaving the source. The upcoming upgrade of the Pierre Auger Observatory, AugerPrime [49,50], will improve our knowledge on the composition of the cosmic rays reaching Earth, and therefore future LIV searches will profit enormously from its operation.    Alternative UHECR scenario | n=2 Figure 3. Simulated integral flux of GZK photons as a function of the energy for an alternative scenario with subdominant proton component [40]. Continuous lines show the rejected LIV scenarios. The arrows show the flux determined by analysis of the Pierre Auger Observatory data [19,20].    Figure 6. Evolution of the fit parameters with respect to the LIV coefficient, δ had,0 . The panels show the spectral index, rigidity cutoff and integral fractions of each nuclei at the source, respectively. In the first two panels, each line represents a propagation model. In the last panel, each line is the percentage of each nucleus summed in the entire energy range for the SPGE model.  Figure 7. Evolution of the total and relative deviance with respect to the LIV coefficient, δ had,0 . Each line represents a propagation model. The top panel show the total deviance, while the bottom panel shows the difference in the deviation relative to the LI case. The model we fit to the data does not describe satisfactorily the measured energy spectrum and X max distributions, thus limiting the interpretation of the minimum deviance as an imprint of LIV in the data as discussed in Section 7.