Propagation and Fluxes of Ultra High Energy Cosmic Rays in f(R) Gravity Theory

We study the effect of diffusion of ultra-high energy cosmic ray (UHECR) protons in the presence of turbulent magnetic fields in the light of the $f(R)$ theory of gravity. The $f(R)$ theory of gravity is a successful modified theory of gravity in explaining the various aspects of the observable Universe including its current state of expansion. Here, we consider the two most studied $f(R)$ gravity models, viz., the power-law model and the Starobinsky model. For these two models, the diffusive character of the propagation of UHECR protons is explored in terms of their density enhancement (a measure of how the density of UHECRs changes due to their diffusion and interactions in the intergalactic medium). Ankle, instep, and Greisen-Zatsepin-Kuzmin cutoff are spectral characteristics that extragalactic UHECRs acquire when they propagate through the cosmic microwave background. All these characteristics are analysed through the diffusive flux as well as its modification factor. We compare the UHECR protons spectra calculated for the considered $f(R)$ gravity models with the available data of the Telescope Array and Pierre Auger Observatory. Both models of $f(R)$ gravity predict energy spectra of UHECRs with all experimentally observed features, which lay well within the range of combined data of both experiments. However, this work is only to investigate the possible effects of $f(R)$ gravity theory on the UHECRs propagation, using pure proton composition as a simplified case study. At this stage, our results cannot be used to favour or disfavour the $f(R)$ cosmology over the $\Lambda$CDM cosmology.


I. INTRODUCTION
The discovery of cosmic rays (CRs) by V. F. Hess in 1912 [1] is one of the most significant milestones in the history of modern physics.CRs are charged ionizing particles, mostly consisting of protons, helium, carbon and other heavy ions up to iron emanating from outer space.Although the discovery of CRs occurred more than 110 years ago, the origin, acceleration, and propagation mechanisms of CRs are still not clearly known [2][3][4], especially in the higher energy range i.e. the energy range E ≥ 0.1 EeV (1 EeV = 10 18 eV).The sources of such usually referred ultra-high energy CRs (UHECRs) are not established yet [5][6][7][8].However, in the energy range E ≤ 0.1 EeV, it is assumed that the sources are of galactic origin and they are accelerated by supernova explosions [9], while those well above this range (∼ 1 EeV and above) are most probably extragalactic in origin and plausibly to accelerate in gamma-ray (γ-ray) bursts or in active galaxies [2].
The energy spectrum of CRs has an extraordinary range of energies.It extends over many orders of magnitude from GeV energies up to 100 EeV and exhibits a power-law spectrum.There is a small spectral break known as the knee at about 4 PeV (1 PeV = 10 15 eV) and a flattening at the ankle at about 5 EeV.In this spectrum, a strong cutoff near 50 EeV, which is called the Greisen-Zatsepin-Kuzmin (GZK) cutoff [10,11] appears due to the interaction with cosmic microwave background (CMB) photons.
The intergalactic medium (IGM) contains turbulent magnetic fields (TMFs), which impact significantly the propagation of extragalactic UHECRs.In the presence of any random magnetic field, the propagation of a charged particle depends on how much distance is traveled by that particle compared with the scattering length λ = 3D/c in the medium, where D denotes the diffusion coefficient and c is the speed of light in free space [12].If the traveled distance of the charged particle is much smaller than the scattering length, then the propagation is ballistic in nature while that is diffusive if the distance is much larger than the scattering length.Consideration of an extragalactic TMF and also taking into account the finite density of sources in the study of the propagation of UHECRs may result in a low-energy magnetic horizon effect, which may allow the observations to be consistent with a higher spectral index [9,13,14], closer to the values anticipated from diffusive shock acceleration.Other hypotheses rely on the assumption of acceleration of heavy nuclei by extragalactic sources, which then interact with the infrared radiation present in those environments to photodisintegration, producing a significant number of secondary nucleons that might explain the light composition seen below the ankle [15,16].In the presence of an intergalactic magnetic field, the propagation of UHECRs can be studied from the Boltzmann transport equation or by using some simulation methods.In Ref. [12], the author presents a system of partial differential equations to describe the propagation of UHCRs in the presence of a random magnetic field.In that paper, the author considered the Boltzmann transport equation and obtained the partial differential equations for the number density as well as for the flux of particles.A diffusive character of the propagation of CRs is also obtained in that paper.In Ref. [17] (see also Ref. [18]), an astrophysical simulation framework is proposed for studying the propagating extraterrestrial UHE particles.In their work, authors presented a new and upper version of publicly available code CRPropa 3. It is a code for the efficient development of astrophysical predictions for UHE particles.Ref. [19] presented an analytical solution of the diffusion equation for high-energy CRs in the expanding Universe.A fitting to the diffusion coefficient D(E) obtained from numerical integration was presented in Ref. [2] for both Kolmogorov and Kraichnan turbulence.Authors of Ref. [3] studied the effects of diffusion of CRs in the magnetic field of the local supercluster on the UHECRs from a nearby extragalactic source.In that study, the authors found that a strong enhancement at certain energy ranges of the flux can help to explain the features of the CR spectrum and the composition in detail.In Ref. [5], the authors demonstrated the energy spectra of UHECRs as observed by Fly's Eye [20], HiRes [21] and AKENO [22] from the idea of the UHE proton's interaction with CMB photons.A detailed analytical study of the propagation of UHE particles in extragalactic magnetic fields has been performed in Ref. [23] by solving the diffusion equation analytically with the energy losses that are to be taken into account.In another study [24], the authors obtained the ankle, instep and GZK cutoff in terms of the modification factor, which arises due to various energy losses suffered by CR particles while propagating through the complex galactic or intergalactic space [4].Similarly, in Ref. [25], authors obtained four features in the CR proton spectrum, viz. the ankle, instep, second ankle, and the GZK cutoff taking into consideration of extragalactic proton's interaction with CMB and assuming of resulting power-law spectrum.
General relativity (GR) developed by Albert Einstein in 1915 to describe the ubiquitous gravitational interaction is the most beautiful, well-tested, and successful theory in this regard.The discovery of gravitational waves (GWs) by LIGO detectors in 2015 [26] after almost a hundred years of their prediction by Einstein himself and the release of the first image of the supermassive black hole at the center of the elliptical supergiant galaxy M87 by the Event Horizon Telescope (EHT) in 2019 [27][28][29][30][31][32] are the robust supports amongst others in a century to GR.Even though the GR has been suffering from many drawbacks from the theoretical as well as observational fronts.For example, the complete quantum theory of gravity has remained elusive till now.The most important limitations of GR from the observational point of view are that it can not explain the observed current accelerated expansion [33][34][35][36] of the Universe, and the rotational dynamics of galaxies indicating the missing mass [37] in the Universe.Consequently, the modified theories of gravity (MTGs) have been developed as one of the ways to explain these observed cosmic phenomena, wherein these phenomena are looked at as some special geometrical manifestations of spacetime, which remain to be taken into account in GR.The most simplest but remarkable and widely used MTG is the f (R) [38] theory of gravity, where the Ricci scalar R in the Einstein-Hilbert (E-H) action is replaced by a function f (R) of R. Various models of f (R) gravity theory have been developed so far from different perspectives.Some of the viable as well as famous or popular models of f (R) gravity are the Starobinsky model [39,40], Hu-Sawicki model [41], Tsujikawa model [42], power-law model [43] etc.
Till now several authors have studied the propagation of CRs in the domain of GR [2-9, 12, 23, 24, 44].The enhancement of the flux of CRs is obtained in the framework of the ΛCDM model by a variety of authors [3,45].Besides these, differential flux as well as the modification factor have also been studied [4,5,12,[23][24][25]. Since MTGs have made significant contributions to the understanding of cosmological [46,47] and astrophysical [48] issues in recent times, it would be wise to apply the MTGs in the field of CRs to study the existing issues in this field.Keeping this point in mind, in this work, we study for the very first time the propagation of UHECRs and their consequent flux in the light of an MTG, the f (R) theory of gravity.For this purpose, we consider two f (R) gravity models, viz. the power-law model [43] and the Starobinsky model [40].Considering these two models, we calculate the expression for the number density of particles.From the number density, we calculate the enhancement factor as well as the differential flux and modification factor for the UHECRs.
The remaining part of the paper is arranged as follows.In Section II, we discuss the turbulent magnetic field and diffusive propagation mechanism.The basic cosmological equations that are used to calculate the cosmological parameters are introduced in Section III.In Section IV, we define f (R) gravity models of our interest and calculate the required cosmological parameters for those models.The fittings of predicted Hubble parameter values at different cosmological redshifts by those models to the observational Hubble parameter data are also shown in this section.In Section V, we calculate the number density of particles and hence the enhancement factor.Then the differential fluxes for both models along with the ΛCDM model were calculated and compared these results with the data of the Telescope Array (TA) experiment [49] and Pierre Auger Observatory (PAO) [50].We also compare the calculated modification factors for those two models with the observational data of TA and PAO.Finally, we compare the results of all three models including the ΛCDM model, and then conclude our paper with a fruitful discussion incorporating the Chi-square test in Section VI.

II. PROPAGATION OF COSMIC RAYS IN TURBULENT MAGNETIC FIELDS
It is a challenging task to build a model for the extragalactic magnetic fields since there are few observable constraints on them [51].Their exact amplitude values are unknown, and they probably change depending on the region of space being considered.In the cluster center regions, the large-scale magnetic fields have recorded amplitudes that vary from a few to tens of µG [52].
Smaller strengths are anticipated in the vacuum regions, with the typical boundaries in unclustered regions being 1 to 10 nG.This means that considerable large-scale magnetic fields should also be present in the filaments and sheets of cosmic structures.The evolution of primordial seeds impacted by the process of structure building may result in TMFs in the Universe [2].As a result, magnetic fields are often connected with the matter density and are therefore stronger in dense areas like superclusters and weaker in voids (≤ ∼ 10 −15 G).In the local supercluster region, a pragmatic estimation places the coherence lengths of magnetic fields in between 10 kpc and 1 Mpc, while their root mean square (RMS) strengths lie in the range of 1 to 100 nG [52][53][54].The regular component of the galactic magnetic field (GMF), which typically has a strength of only a few µG, may have an impact on the CRs' arrival directions, but due to its much lesser spatial extent, it is anticipated to have a subdominant impact on the CRs spectrum.
In the local supercluster region, the rotation measures of polarised background sources have suggested the presence of a strong magnetic field, with a potential strength of 0.3 to 2 µG [53].It is the magnetic field within the local supercluster that is most relevant since the impacts of the magnetic horizon become noticeable when the CRs from the closest sources reach the Earth.Thus we will not consider here the larger-scale inhomogeneities from filaments and voids.The propagation of CRs in an isotropic, homogenous, turbulent extragalactic magnetic field will then be simplified.The rms amplitude of magnetic fields B, and the coherence length l c which depicts the maximum distance between any two points up to which the magnetic fields correlate with each other, can be used to characterize such magnetic fields.The RMS strength of magnetic fields can be defined as B = ⟨B 2 (x)⟩, which can take values from 1 nG up to 100 nG and the strength of the coherence length l c can take the values from 0.01 Mpc to 1 Mpc.
An effective Larmor radius for charged particles of charge Ze moving with energy E through a TMF of strength B may be defined as A pertinent quantity in the study of diffusion of charged particles in magnetic fields is the critical energy of the particles.This energy can be defined as the energy at which the coherence length of a particle with charge Ze is equal to its Larmor radius i.e., r L (E c ) = l c and it is given by This energy distinguishes between the regime of resonant diffusion that occurs at low energies (< E c ) and the non-resonant regime at higher energies (> E c ).In the resonant diffusion regime, particles suffer large deflections due to the interaction with magnetic field B with scales that are comparable to l c , whereas in the latter scenario, deflections are small and can only take place across travel lengths that are greater than l c .Extensive numerical simulations of proton's propagation yielded a fit to the diffusion coefficient D as a function of energy [2], which is given by where m is the index parameter, a I and a L are two coefficients.For the case of TMF with Kolmogorov spectrum, m = 5/3 and the coefficients are a L ≈ 0.23 and a I ≈ 0.9, while that for Kraichnan spectrum one will have m = 3/2, a L ≈ 0.42 and a I ≈ 0.65.The diffusion length l D relates to the distance after which overall deflection of particles is nearly one radian and is given by l D = 3D/c.From Eq. (3), it is seen that for

III. BASIC COSMOLOGICAL EQUATIONS
On a large scale, the Universe appears to be isotropic and homogeneous everywhere.In light of this, the simplest model to be considered is a spatially flat Universe, which is described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric and is defined as where a(t) is the scale factor, δ ij is the Kronecker delta function with i, j = {1, 2, 3} and x µ = {x 0 , x 1 , x 2 , x 3 } are comoving coordinates with x 0 = t.Moreover, as a source of curvature, we consider the perfect fluid model of the Universe with energy density ρ and pressure p which is specified by the energy-momentum tensor T µ ν = diag(− ρ, p, p, p).At this stage, we are interested in the basic cosmological evolution equation to be used in our study and this equation is the Friedmann equation.
The Friedmann equation in f (R) gravity theory is derived by following the Palatini variational approach of the theory.In this approach both the metric g µν and the torsion-free connection Γ λ µν are considered as independent variables.In our present case the metric is g µν = diag(− 1, a 2 , a 2 , a 2 ) and the connection can be obtained from the f (R) gravity field equations in the Palatini formalism [38].Following the Palatini formalism the generalized Friedmann equation for our Universe in terms of redshift in f (R) gravity theory can be expressed as [55] H where In Eqs. ( 5) and ( 6), H 0 ≈ 67.4 km s −1 Mpc −1 [56] is the Hubble constant, Ω m0 ≈ 0.315 [56] is the present value of the matter density parameter and Ω r0 ≈ 5.373 × 10 −5 [57] is the present value of the radiation density parameter.f ′ (R) and f ′′ (R) are the first and second-order derivatives of the function f (R) with respect to R. It is seen that Eqs. ( 5) and ( 6) are f (R) gravity model dependent.
Secondly, in our study, it is important to know how the cosmological redshift is related to the cosmological time evolution.This can be studied from the connection between the redshift and cosmological time evolution, which is given by The expression of the Hubble parameter H(z) for different models of f (R) gravity will be derived using Eqs.( 5) and ( 6) in the next section IV.

IV. f (R) GRAVITY MODELS AND COSMOLOGICAL EVOLUTION
In this section, we will introduce the power-law model [43] and Starobinsky model [40] of f (R) theory of gravity, and then will derive the expressions for the Hubble parameter and evolution Eq. ( 7) for these two models.The least-square fits of the derived Hubble parameters for the models to the recent observational data will also be done here to constrain the parameters of the models.

A. Power-law model and cosmological equations
The general f (R) gravity power-law model is given by [58,59] where λ and n are two model parameters.Here the parameter n is apparently a constant quantity, but the parameter λ depends on the value of n as well as on the cosmological parameters H 0 , Ω m0 and R 0 as given by [59] This expression of the parameter λ implies that the power-law model has effectively only one unknown parameter, which is the n.For this model, the expression of the present value of the Ricci scalar R 0 can be obtained as [59] The expression of the Hubble parameter H(z) for the power-law model can be obtained from Eq. ( 5) together with Eq. ( 6) as [59] H In our study for the model parameter n, we use its value from the Ref. [59] where a detailed study has been made on this model in the cosmological perspective and the values of n = 1.25, 1.4 and 1.9 have been taken into account.Among these values the best-fitted value of n is 1.4 according to this Ref.[59].
The relation between the cosmological evolution time t and redshift z for the power-law model can be obtained by substituting Eq. ( 11) for H(z) in Eq. ( 7) as given by In Fig. 1, we plot the differential variation of cosmological time t with respect to redshift z i.e. the variation of dt/dz with the redshift z for different values for model parameter n along with that for the ΛCDM model.It is seen from difference of variation of dt/dz for the power-law model from the ΛCDM model is both redshift z and model parameter n dependent.The difference appears to be less significant for all values of n when values of z < 0.2, while for higher values of z, it has shown a notable deviation depending on the value of n.However, at around certain higher values z the power-law model predicts the same values of dt/dz as that of ΛCDM model depending on the parameter n.For example, at around z = 2.8 the power-law model with n = 1.4 and the ΛCDM model predict the same dt/dz.Beyond such values of z corresponding to n values the prediction of the power-law model deviates significantly from the ΛCDM model.It should be mentioned that although n = 1.4 is found as the most suitable value of the parameter of the power-law model as informed earlier, we use other two values of n in this plot to see how the model prediction varies from that of the ΛCDM model with different values of n.It is clear that the higher values of n obviously show more deviation from the ΛCDM model prediction for all appropriate values of z and hence the most favorable value n = 1.4 shows appreciable behavior in this regard.

B. Starobinsky Model and cosmological equations
The Starobinsky model of f (R) gravity considered here is of the form [40]: where α and β are two free model parameters to be constrained by using observational data associated with a particular problem of study.Similar to the previous case the expression of the Hubble parameter H(z) for the Starobinsky model can be obtained from Eq. ( 5) along with Eq. ( 6) as To use this expression of H(z) for further study we have to constrain the values of the model parameters α and β within their realistic values as the behaviour of H(z) depends significantly on these two model parameters.For this, we use the currently available observational Hubble parameter (H obs (z)) data set [60] as shown in Table I.Here we consider the combination of 43 observational Hubble parameter data against 43 distinct values of redshift z (as they are available in the references mentioned) to obtain the precise values of the aforementioned free model parameters, so that the predicted H(z) should be consistent with the ΛCDM model value at least around the current epoch i.e. at z ∼ 0. Using the least square fitting technique in ROOT software [61], we plot the best-fitted curve to this set of Hubble parameter data with respect to redshift as shown in Fig 2 .For this leastsquare fitting, we use an exponential function of the form: a exp(bz), where a and b are two constants whose values are found after the fitting as a = 69.750± 0.927 and b = 0.503 ± 0.012.Using this fitting we infer values of α and β as 1.07 and 0.00086 respectively by using the chi-square minimization method (as discussed in [59]).The value of χ 2 is 29.38 with the critical value (in 95 % confidence level) 56.94.Now, we are in a position to write the expression for dt/dz for this model and it can be expressed as In Fig. 3  model with the constrained set of parameters predicts the values of dt/dz which are almost comparable to the values of the same predicted by the ΛCDM model over the considered range of z (especially for z > 1).Whereas, except at z ∼ 2.8 there is a noticeable difference in the prediction of the power-law model from that of the ΛCDM model, although the difference is small at z < 0.2 as mentioned already.Power-law model predicts lower values of dt/dz than that for the other two models from z = 0.1 to z ∼ 2.8 and above this range the trend becomes reversed.Moreover, it is found that at z = 0, i.e. at the present epoch the ΛCDM model predicts the highest value and the power-law model predicts the lowest value of dt/dz.
In the next section, we will employ the results of this section to calculate the density and differential flux of CRs for the power-law and Starobinsky models.The first thing that piques our curiosity is how the density of CRs is being modulated at a certain distance from the originating source in a TMF.For this, it is necessary to calculate the density enhancement of CRs at a certain distance r s from the originating source while being surrounded by a TMF.Specifically, we wish to investigate the reliance of density enhancement on different CR parameters taking into account the diffusive propagation of CRs in the light of f (R) gravity theory.
In the diffusive regime, the diffusion equation for UHE particles propagating in an expanding Universe from a source which is located at a position x s can be expressed as [19] where H(t) = ȧ(t)/a(t) is the Hubble parameter as a function of cosmological time t, ȧ(t) is the time derivative of the scale factor a(t), x denotes the comoving coordinates, ρ is the density of particle at time t and position x, Q s (E) is the source function that depicts the number of emitted particles with energy E per unit time.Thus, at time t, which corresponds to redshift z, r s = x − x s .The energy losses of particles due to expansion of the Universe and interaction with CMB are described by Here H(t)E represents the adiabatic energy losses due to expansion and b int (E) denotes the interaction energy losses.The interaction energy losses with CMB include energy losses due to pair production and photopion production (for details see [2]).The general solution of Eq. ( 16) was obtained in Ref. [19] considering the particles as protons and it is given as where z i is the redshift at the initial time when a particle was just emitted by a source and E g is the generation energy at redshift z of a particle whose energy is E at z = 0, i.e. at present time.The source function with γ g as the spectral index of generation at the source.λ is the Syrovatsky variable [13,75] and is given by Here λ(E, z) refers to the usual distance that CRs travel from the location of their production at redshift z with energy E g , to the present time at which they are degraded to energy E. The expression of the rate of degradation of energy of particles at the source with respect to their energy at z = 0, i.e. dE g /dE is given by [19,25] The detailed derivation of this expression was nicely performed by Berezinsky et al. in Appendix B of Ref. [25].It is clear that using Eqs.( 12) and ( 15) in Eqs. ( 19) and ( 20) the density of UHE protons in the diffusive medium at any cosmological time t with energy E and at a distance r s from the source can be obtained for the f (R) gravity power-law model and the Starobinsky model respectively, as given by Eq. ( 18).So, in the following, we will implement the results of the power-law and Starobinsky models from Section IV to obtain the CR protons density enhancement factor, and subsequently their flux and energy spectrum as predicted by these two f (R) gravity models.

A. Projections of f(R) power-law model
To calculate the CR protons density from Eq. ( 18) and hence its enhancement factor in the TMF of extragalactic space projected by the power-law model of f (R) gravity, as a prerequisite we calculate first the Syrovatsky variable λ for this model from Eq. ( 19) using Eq.(12).In this calculation we use different values of the model parameters n taking the feasible values of field parameters as l c = 0.1 Mpc and B = 50 nG with the corresponding critical energy of protons as E c = 4.5 EeV.Then we study the behaviour of the variable λ for the both Kolmogorov spectrum and the Kraichnan spectrum.We also calculate this variable for the ΛCDM model for those two spectra for comparison.Here and rest of calculations we use the values of z = 0 − 5 keeping in view of possible source locations of CRs as well as the present and probable future cosmological observable range.The results of these calculations are shown in Fig. 4 with respect to energy E for the Kolmogorov spectrum (m = 5/3, a L ≈ 0.23 and a I ≈ 0.9) (left panel) and the Kraichnan spectrum (m = 3/2, a L ≈ 0.42 and a I ≈ 0.65) (middle panel).It is seen from the figure that the value of λ 2 increases substantially with increasing energy of particles.The power-law model predicts higher values of λ 2 for all values of n in comparison to that of the ΛCDM models for both spectra and this difference increases significantly with the increasing energy E. Similarly higher values of the parameter n give increasingly higher values λ 2 in comparison to the smaller values of n.No difference can be observed between the values of λ 2 obtained for the Kolmogorov spectrum and the Kraichnan spectrum from the respective plots.So, to quantify the difference of values of λ 2 for these two spectra to a visible one we calculate the percentage of per average bin difference between λ 2 values obtained for the Kolmogorov spectrum and the Kraichnan spectrum in each energy bin (∆λ 2 kk (%)) for the power-law model with n = 1.4,which is shown in the right panel of the figure.A peculiar behaviour of the variation of ∆λ 2 kk (%) with energy is seen from the plot.The ∆λ 2 kk (%) is energy dependent, it decreases rapidly with E up to ∼ 0.4 EeV, after which it shows oscillatory behaviour with the lowest minimum at ∼ 1.55 EeV.At energies above 0.1 EeV, the values of ∆λ 2 kk (%) are seen to be mostly below the 1%.Thus at these UHEs differences of λ 2 values for the Kolmogorov spectrum and the Kraichnan spectrum are not so significant.
In the diffusive regime, the density of particles has been enhanced by a factor depending on the energy, distance of the particles from the source and TMF properties.The density enhancement factor can be defined as the ratio of actual density to the density of particles that would lead to their rectilinear propagation, which is given by [3] ξ where L(E) is the spectral emissivity of the source, which has a power-law dependency on the energy of the particles.The results of the enhancement of the density for a proton source and for various parameters values obtained from Eq. ( 21) by numerically integrating Eq. ( 18) are displayed in Fig. 5.The distance to the source r s , the magnetic field amplitude B, and its coherence length l c are the major factors that determine the lower-energy suppression of the density enhancement factor.In the lower panels, the enhancement energy range is less as compared to the upper panels, which is lowest in the case of the lower right panel.As the distance from the source is far away, the enhancement of density is limited to a smaller range of energies but shifted towards the higher energy side.The final verdict from Fig. 5 is that as the distance from the source r s increases, the enhancement becomes gradually model independent.Also one can appreciate that the f (R) gravity power-law model has done a perfect job by enhancing density in a wider range of energies as compared to the ΛCDM model.For a given source distance of 25 Mpc and coherence length of 0.1 Mpc, we depict the enhancement factor ξ as a function of E/E c in Fig. 6 to better highlight the fact that for E/E c < 0.01 the Kolmogorov spectrum (left panel) and Kraichnan spectrum (right panel) have shown different behaviours while for E/E c > 0.01 both Kolmogorov and Kraichnan spectra have shown similar patterns.In this case, the f (R) power-law model is more suitable as it gives the enhancement in the higher as well as lower values of E/E c , while in the case of the ΛCDM the range it gives the enhancement is less wide than the power-law model.From this Fig. 6 it is seen that the Kolmogorov spectrum has given a better range of E/E c than the Kraichnan spectrum for both ΛCDM and f (R) power-law model.The diffusive character of the propagation of UHE protons is shown in Fig. 7.Here we plot the density enhancement factor ξ as a function of source distance r s .In these plots, we fix the coherence length l c = 0.1 Mpc, while energy E = 0.1 to 5 EeV have been taken into account.From these plots, we can say that the lower E/E c value results in a higher peak of the density enhancement with the peak position towards the smaller value of r s , and also the enhancement peak lies in the diffusive region for smaller E/E c value.Again, the ΛCDM model shows the highest peak in the CRs density enhancement, while the f (R) gravity power-law model depicts a better distribution of enhancement with the source distance.The power-law model with parameter values n = 1.25 and 1.4 results in a similar distribution, while for n = 1.9, it shows a larger distribution.In the lower right panel, we consider a larger value of E/E c which results in a very poor peak for both ΛCDM and f (R) power-law models.In this case, the enhancement peak is very far away from the diffusive regime.So from these results, we can finally say that for the suitable values of l c and E/E c , the ΛCDM model depicts a better peak, while f (R) power-law model depicts the enhancement in a much wider distribution.For a better illustration, we also draw contour plots of the density enhancement with the source distance r s (0 − 100 Mpc) and coherence length l c (0.05 − 0.5 Mpc) taking E = 0.1 EeV and E c = 4.5 EeV for the ΛCDM and f (R) power-law models as shown in Fig. 8.One can see that the density enhancement depends on the coherence length l c also.For the higher value of l c the density enhancement decreases, shifts away its maximum value from the source and takes place in the non-diffusive regime.
For reckoning the diffuse spectrum of UHE particles the separation between sources plays a crucial role.If the sources are distributed uniformly with separations, which are much smaller than the propagation and interaction lengths, then the diffuse spectrum of UHE particles has a universal form, regardless of the mode of propagation of such particles [23].To this end the explicit form of the source function Q(E, z) for the power-law generation of the particles can be written as [25] Q(E, z) = L 0 (1 + z) δ Kq gen (E g ), (22) where L 0 = L(E) dE is the total emissivity, (1 + z) δ represents the probable cosmological evolution of the sources with an index parameter δ, K is a normalisation constant with K = γ g − 2 for γ g > 2 and for γ g = 2, K = (ln E max /E min ) −1 , and q gen = E −γg g (see Appendix A for E g ).Utilizing the formalism of Ref. [19], it is possible to determine the spectrum of UHE protons in the model with a uniform source distribution and hence one can obtain the diffuse flux of UHE protons as Following Eq. ( 12) one can rewrite this diffuse flux Eq. ( 23) as The spectrum given by Eq. ( 23) is known as the universal spectrum as it is independent of the mode of propagation of particles which is the consequence of the small separation of sources as mentioned earlier.The shape of the universal spectrum may theoretically be changed by a variety of effects, which include fluctuations in interaction, discreteness in the source distribution, large-scale inhomogeneous source distribution and local source overdensity or deficit.However, the aforementioned effects only slightly change the form of the universal spectrum, except for energies below 1 EeV.Numerical simulations demonstrate that the energy spectrum is changed by the propagation of UHE protons in the strong magnetic fields depending on the separation of sources.For small separation of sources with their uniform distribution the spectrum becomes the universal one as mentioned already [76,77].In Fig. 9, we plot the diffusive flux with no cosmological evolution (δ = 0) [23,25].The emissivity L 0 is taken to fit the curve with the available observational data [25].The energy-rescaling data of the TA experiment and PAO have been taken from Ref. [78].It needs to be mentioned that the energy rescaling in the data of these two experiments is used to avoid the effect of the difference in the energy scales used by these two observatories.The uncertainty present in the energy scale contributes a significant impact on the uncertainty in the normalisation of the spectrum [78].The considered f (R) gravity power-law model has shown a very good agreement with the observational data in predicting the energy spectrum of UHECRs and has also predicted similar result with that of the ΛCDM one.However, only a slightly higher flux is obtained for the powerlaw model in comparison to the ΛCDM model above 4 EeV.In data, a dip (the ankle) is seen at the energy around 4.5 EeV, while at about 30 EeV a bump (the instep) is observed.The ankle predicted by both power-law and ΛCDM models is at slightly lower energy around 3.5 EeV, but the predicted position of the instep is the same as that of the data.These two signatures, the ankle and the instep are also observed in the modification factor of the energy spectrum plot as shown in Fig. 10.The modification factor of the energy spectrum is a convenient parameter for analysing the energy spectrum of UHECRs.This parameter corresponds to the enhancement factor of the density of UHECR particles discussed earlier.The modification factor of energy spectrum η(E) is calculated as the ratio of the universal spectrum J p (E) which accounts for all energy losses to the unmodified spectrum J unm p (E), in which only adiabatic energy losses due to the redshift are taken into consideration [25], i.e.Without any cosmological evolution, the unmodified spectrum can be written as The modification factor as a function of energy with the spectral index γ g = 2.7 is shown in Fig. 10 for the f (R) gravity powerlaw model and the ΛCDM model.At about 1 EeV, the ankle is seen in the spectrum as predicted by both models in agreement with the observation of the TA experiment [49] and PAO [50] as well as a good agreement for the instep in the spectrum is also seen.From Fig. 10, it can also be said that the modification factor of the energy spectrum is a weak model-dependent parameter.

B. Projections of Starobinsky f(R) gravity model
For this model of f (R) gravity also we will follow the same procedure as we have already done in the case of the power-law model.So here also we have to calculate the Syrovatsky variable λ 2 and for this purpose, we express λ 2 (E, z) from Eq. ( 19) 10 19  10 20 E (eV)  FIG. 10.Spectra of modification factor with γg = 2.7 for the f (R) gravity power-law model (n = 1.4) and the ΛCDM model in comparison with experimental data of the TA experiment [49] and PAO [50].using Eq. ( 15) for the Starobinsky model as In Fig. 11 we plot the variation of λ 2 with respect to energy for the f (R) gravity Starobinsky model and power-law model in comparison with the ΛCDM model.For this, we consider the source distance r s = 50 Mpc, the coherence length l c = 0.1 Mpc and the strength of the TMF, B = 50 nG, and use only the Kolmogorov spectrum of the diffusion coefficient.A noticeable variation with respect to the energy is observed in λ 2 values for all of the mentioned gravity models.Moreover, the f (R) gravity Starobinsky model gives the lowest value of λ 2 although its pattern of variation with respect to energy is similar for all three models.Similarly, using Eqs.( 18), ( 20) and ( 27) in Eq. ( 21) we calculate the density enhancement factor ξ(E, r s ) of UHE particles for the Starobinsky model, which can be written as the parameters we consider and for the different parameters we find a very distinct result in each of the cases.The distinction between enhancement factors for the Starobinsky model and ΛCDM model is clearly visible.The Starobinsky model gives a higher peak and wider range of the enhancement factor than that given by the ΛCDM model.Moreover, for smaller to medium values of r s the difference between the two models on the higher energy side is very small, while for higher values of r s , it is very small on the lower energy side of the enhancement factor plots.
Similar to the case of the f (R) gravity power-law model, here also we plot the density enhancement factor ξ with respect to the source distance r s by keeping fixed the coherence length l c = 0.1 Mpc for E = 0.1 EeV (upper left panel), E = 0.5 EeV (upper right panel), E = 1 EeV (lower left panel) and E = 5 EeV (lower right panel) in Fig. 13 to understand the propagation of UHECR protons in the light of the Starobinsky model in comparison with the ΛCDM model.From this figure, one can see that similar to the power-law model the peak of the enhancement is higher for smaller values of E/E c , whereas the peak of the distribution is higher for the Starobinsky model than that of the ΛCDM model.Also similar to the power-law model the ξ distribution becomes wider and the peak of it is shifted away from the source for higher E/E c values.As in the previous case for a clear understanding of the diffusive propagation, here also we draw contour plots in Fig. 14 for enhancement by keeping a range for the coherence length l c from 0.05 − 0.5 Mpc and that of the source distance r s from 0 − 100 Mpc for E = 0.1 EeV and E c = 4.5 EeV.We see that for increasing the value of l c the enhancement decreases and also shifts its maximum value from    [49] and PAO [50] along with the flux for the ΛCDM model.
The diffuse UHECR protons flux for the f (R) gravity Starobinsky model can be expressed as In Fig. 16, we plot this flux (29) as a function of energy by considering the Starobinsky model parameters as we discuss in Section IV.From the figure, we can see that the Starobinsky model predicts a spectrum that is also in good agreement with the TA experiment and PAO data.However, within 4 EeV to 8 EeV, the Starobinsky model's predicted spectrum remains slightly above the observational data range.It also gives noticeably higher flux in comparison to the ΛCDM model over the almost entire energy range considered and the difference increases with increasing energy.Moreover, it also predicts the ankle of the spectrum at a lower energy of around 3.5 EeV than data at the energy of around 4.5 EeV, but the position of the predicted instep remains as that of the data.A detailed comparison of the diffuse fluxes for all the models considered in this work will be discussed in the next section.Finally, for the calculation of the modification factor η of the energy spectrum, the unmodified flux of UHECR protons for the Starobinsky model is given by Fig. 17 shows the behaviour of the modification factor for the Starobinsky model along with that of the ΛCDM model, and is compared with experimental data as in the previous case.The observational data have given a good agreement with the calculated modification factor spectrum with the ankle as well as the instep for the Starobinsky model similar to the ΛCDM model.It is also clear that the modification factor is very weakly model dependent as seen in the case of the power-law model also.

VI. DISCUSSIONS AND CONCLUSIONS
The believable sources of UHECRs are extragalactic in origin [2,79].Accordingly, the propagation mechanisms of UHECRs through the extragalactic space have been one of the prime issues of study for the past several decades.It can be inferred that in the propagation of UHECRs across the extragalactic space, the TMFs that exist in such spaces and the current accelerated expansion of the Universe might play a crucial role.Thus this idea led us to study the propagation of UHECRs in the TMFs in the extragalactic space in the light of f (R) theory of gravity and to compare the outcomes with the experimental data of two world-class experiments on UHECRs.The f (R) theory of gravity is the simplest and one of the most successful MTGs FIG. 17. Spectrum of the modification factor for the f (R) gravity Starobinsky model along with that for the ΛCDM model with γg = 2.7, which is in comparison with the experimental data of TA experiment [49] and PAO [50].
that could explain the current accelerated expansion of the Universe.To this end, we consider two f (R) gravity models, viz., the power-law model and the Starobinsky model.The Starobinsky model of f (R) gravity is the most widely used and one of the most viable models of the theory [40,47,59].Similarly, the power-law model is also found to be suitable in various cosmological and astrophysical perspectives [59].The basic cosmological equations for these two f (R) gravity models, which are required for this study are taken from the Ref. [59].Independent parameters of the models are first constrained by using the recent observational Hubble data.The relation between the redshift z and the evolution time t is calculated for both models.The UHECRs density ρ(E, r s ) and hence the enhancement factor of the density ξ(E, r s ) are obtained and they are calculated numerically for both the models of f (R) gravity.A comparative analysis has been performed between the predictions of the power-law model and Starobinsky model of f (R) gravity along with the same of the ΛCDM model for the density enhancement factor ξ as a function of source distance r s in Fig. 18.In this analysis, we consider the coherence length l c = 0.1 Mpc and the fraction of energy and critical energy E/E c = 0.02.One can observe that at r s < 20 Mpc, the variation of ξ for the Starobinsky model and the ΛCDM is not very different but at the far distance from the source, the behaviour of these two models is quite different in terms of the peak position of the enhancement and the range of the source distance where the enhancement takes place.In the case of the f (R) power-law model, the enhancement is less than the Starobinsky model and the ΛCDM model, but it gives the density enhancement in a much wider range than the ΛCDM model.In fact, it gives the same range of source distance distribution in the enhancement and gives the peak of enhancement at the same distance as that of the Starobinsky model although the enhancement is comparatively low.Another comparative analysis has been done in Fig. 19  the best results as compared with the other two models.From the left panel of Fig. 19, we see that at lower energies i.e., below 1 Eev, the enhancement is different for different energy values including the peaks for all three models.But if we take a look at the higher values of energy, all three models depict almost similar results in the enhancement.One can say that the maximum value of enhancement for the power-law model and the ΛCDM model is approximately the same but the power-law model has covered a wider range of energy values than the ΛCDM model.While the Starobinsky model gives the highest enhancement value as well as the enhancement in a much wider range of energy values.The right panel of Fig. 19 is plotted to show the variation of density enhancement as a function of E/E c .In this panel, we consider the coherence length l c = 0.05 Mpc and source distance r s = 25 Mpc to demonstrate the behaviour of enhancement with the per unit increase of energy with respect to the critical energy.It is seen that at E/E c = 10 −4 , the values of enhancement for the Starobinsky model and ΛCDM are approximately the same, while the f (R) power-law model has shown a higher value of enhancement at this point.But as the fraction of energy is increased, the Starobinsky model has given a better result of enhancement as compared to the other two models.
We calculate the E 3 magnified flux numerically for the both f (R) gravity power-law and Starobinsky models and plot them along with that for the ΛCDM model in Fig. 20 (left panel).We compare our calculations with the available observational datasets of the TA experiment [49] and PAO [50] consisting of 15 and 18 numbers of data points respectively.All of these models have shown a very good agreement with the observational data of both the UHECRs experiments in predicting the signatures of UHECRs energy spectra.The Starobinsky model spectrum has shown a higher flux throughout the energy ranges considered.However, around 30 EeV it gives the flux very near to that of the power-law model.While the power-law model gives almost the same flux as that of the ΛCDM model below 4 EeV, above this energy the power-law model gives gradually higher flux than the ΛCDM model.The shaded regions have depicted the uncertainties in predicting the fluxes by the power-law and Starobinsky models.It is seen that the uncertainty regions in the plot are confined within the error bars of the observational data range.FIG.20.Calculated E 3 magnified spectra of UHECR protons for the ΛCDM model, the f (R) gravity power-law model (n = 1.4), and the Starobinsky model in comparison with the data of TA experiment [49] and PAO [50] with the uncertainty regions for the considered cosmological models (left panel).The modification factors of these spectra are shown in comparison with the TA experiment and PAO data in the right panel.To test the goodness of fit of our model's predictions to the experimental data, we implement the χ 2 test defined as where F i th is the ith theoretical value of flux that we obtained from a cosmological model and F i obs is the ith observational value of flux obtained from the TA experiment or PAO.σ is the standard deviation of the correspodng observed data.The values of χ 2 along with their critical values are shown in Table II.It is seen that the χ 2 value of the fit of predicted fluxes for each model to the data set is small in comparison to the corresponding critical value.Hence, this justifies the trustability of the model's predictions.It is to be noted that the critical value is calculated for the 95% confidence level of each dataset using the Python scipy library [80].
The analysis of the ankle and instep is more convenient with respect to the modification factor and we have compared it with the available data.Both the considered f (R) gravity models have shown a good agreement with the observational data.Thus it can be concluded that the f (R) gravity models considered here are found to be noteworthy with some limitations depending upon the range of energies in explaining the propagations of UHECRs and hence the observed data of their fluxes.However, we would like to clarify here that our aim was to study the possible effects of f (R) cosmology on UHECRs propagation by considering the pure proton composition of UHECRs as a conservative case study and at present our results may not be used to favour or disfavour whether it is the non-standard or standard cosmology, as we need to do more work to confirm our results and to rule out other possible explanations for our findings.Consequently, it is worth mentioning that by extending the work with these models, it would be interesting to study the localised low-scale anisotropies of CRs that arise at their highest energies.So, we keep this as one of the future prospects of study.

Fig. 1 9 FIG. 1 .
FIG. 1. Variation of dt/dz with the redshift z for different values of the power-law model parameter n along with the variation of the same for the ΛCDM model.
FIG.2.Least-square fitting to the observational Hubble data (OHD) as shown in TableIand the best-fitted curve for the Starobinsky model with parameters α = 1.07 and β = 0.00086.Also, a curve for the power-law model with the model parameter n = 1.4 is shown here along with the curve for the ΛCDM model.
FIG.3.Variation of dt/dz with respect to redshift z for both f (R) gravity models (power-law model and Starobinsky model) in comparison with the variation of the same for the ΛCDM model.Here the constrained parameter(s) is(are) used for the associated f (R) gravity model.
V. COSMIC RAYS DENSITY AND FLUX IN THE DOMAIN OF f (R) GRAVITY

FIG. 4 .
FIG. 4. Variations of λ 2 concerning energy E for the Kolmogorov spectrum (left panel) and the Kraichnan spectrum (middle panel) according to the f (R) gravity power-law model and the standard ΛCDM model.These plots are obtained by considering different values of the powerlaw model parameter n with lc = 0.1 Mpc, B = 50 nG and Ec = 4.5 EeV.The right panel shows the percentage of per average bin difference between λ 2 values for the Kolmogorov spectrum and the Kraichnan spectrum in each energy bin as per the power-law model with n = 1.4.Here and in the rest of the corresponding plots we use z = 0 − 5.

For r s = 25
Mpc, l c = 0.5 Mpc, B = 10 nG, and E c = 4.5 EeV (upper left panel), the enhancement has become noticeable for different gravity models in the energy range E < 1 EeV.For the energy range 0.01 < E < 10 EeV, r s = 50 Mpc, l c = 0.1 Mpc, B = 50 nG and E c = 4.5 EeV (upper right panel) are taken into account.In this case, below 1 EeV the variation of enhancement for different gravity models is more distinguished compared to E > 1 EeV.In the lower left panel, r s = 75 Mpc, l c = 0.05 Mpc, B = 40 nG and E c = 1.8 EeV are used to plot the enhancement factor for the ΛCDM and f (R) power-law models, while this is done for r s = 100 Mpc, l c = 0.05 Mpc, B = 80 nG and E c = 3.6 EeV in the lower right panel.

FIG. 6 .
FIG.6.Variation of density enhancement ξ with E/Ec.The left panel is for the Kolmogorov spectrum while the right panel is for the Kraichnan spectrum obtained by considering the ΛCDM and f (R) gravity power-law models with lc = 0.1 Mpc and rs = 25 Mpc.

FIG. 7 .
FIG. 7. Variation of ξ with source distance rs for the ΛCDM model and f (R) power-law model obtained by considering lc = 0.1 Mpc with E = 0.1 EeV (upper left panel), 0.5 EeV (upper right panel), 1 EeV (lower left panel) and 5 EeV (lower right panel).

FIG. 8 .
FIG. 8. Contour plots of variation of density enhancement factor ξ as a function of source distance rs and coherence length lc obtained by considering E = 0.1 EeV and Ec = 4.5 EeV for the f (R) gravity power-law model and the ΛCDM model.

FIG. 9 .
FIG.9.UHECR protons flux is shown for the f (R) gravity power-law model (n = 1.4) and the ΛCDM model in comparison with experimental data of the TA experiment[49] and PAO[50].

FIG. 14 .
FIG. 14. Contour plots of variation of density enhancement factor ξ with respect to source distance rs and coherence length lc obtained by considering E = 0.1 EeV and Ec = 4.5 EeV for the Starobinsky model and the ΛCDM model.

FIG. 15 .
FIG. 15.Variation of density enhancement ξ with E/Ec for the Starobinsky model in comparison with the ΛCDM model.The left panel is for the Kolmogorov spectrum while the right panel is for the Kraichnan spectrum obtained by considering different sets of coherence length lc and source distance rs.

FIG. 16 .
FIG.16.UHECR protons flux is shown for the f (R) gravity Starobinsky model and compared with the experimental data of TA experiment[49] and PAO[50] along with the flux for the ΛCDM model.

FIG. 18 .
FIG.18.Density enhancement factor ξ as a function of source distance rs is shown for the power-law (n = 1.4) and the Starobinsky models of f (R) gravity in comparison with that for the ΛCDM model by considering lc = 0.1 Mpc, E = 0.1EeV and Ec = 4.5 EeV.

FIG. 19 .
FIG. 19.Density enhancement factor ξ as a function of energy E of UHECR protons obtained by considering rs = 50 Mpc, lc = 0.1 Mpc, B = 50 nG and Ec = 4.5 EeV (left panel), and also as a function of E/Ec of the same particles at rs = 25 Mpc with lc = 0.05 Mpc (right panel).Both panels are shown for the power-law model (n = 1.4), the Starobinsky model, and the ΛCDM model.

TABLE II .
χ 2 values of the fits of the predicted UHECR protons energy spectra by the considered f (R) gravity models and the ΛCDM model to the data of TA experiment and PAO, and their associated critical values .