# Magnetic Suppression of Cosmic Rays’ Flux in $\boldsymbol{f(R)}$ and
$\boldsymbol{f(Q)}$ Theories of Gravity

Swaraj Pratim Sarmah  0009-0003-2362-2080<EMAIL_ADDRESS>Umananda Dev Goswami  0000-0003-0012-7549<EMAIL_ADDRESS>Department of
Physics, Dibrugarh University, Dibrugarh 786004, Assam, India

###### Abstract

We investigate the effects of magnetic diffusion on the spectrum of ultra-high
energy cosmic rays (UHECRs) from a cosmological perspective. To this end, we
consider two modified theories of gravity (MTGs), namely, the $f(R)$ gravity
and a symmetric teleparallel gravity, also known as $f(Q)$ gravity. Utilizing
these two MTGs, we calculate the suppression in the flux of UHECRs for a
collection of sources. Non-evolution (NE) and cosmic star formation rate (SFR)
scenarios have been considered in our calculation of the suppression factor.
This study also includes a mixed composition scenario involving the nuclei
upto iron (Fe). Furthermore, we provide a parameterization of the suppression
factor for the proton and also for the mixed compositions within the $f(R)$
and $f(Q)$ theories, considering both NE and SFR scenarios. The influence of
the turbulent magnetic field on the suppression factor is also incorporated in
our work. Comparative analysis of all our results with the standard
$\Lambda$CDM model reveals significant effects of MTGs on the suppression
factor that the $f(R)$ power-law model predicts the lowest suppression factor,
while the $f(Q)$ model predicts the highest, and interestingly the results
from the standard model fall within the range predicted by these two
cosmological models.

Ultra-High Energy Cosmic Rays; Flux; Suppression; Modified Gravity Theory

##  I Introduction

In the annals of physics, the 1912 discovery of cosmic rays (CRs) by V. F.
Hess stands as a pivotal chapter [1]. CRs are charged particles, consisting of
protons, helium, and heavier ions up to iron. Despite more than a century-
spanning quest since their unveiling, the origin and propagation mechanisms of
CRs through the Universe largely remain shrouded in mystery [2, 3, 4], with
their narrative becoming most obscure beyond the energetic threshold of $E\geq
0.1$ EeV ($1$ EeV $=10^{18}$ eV). The provenance of these enigmatic ultra-high
energy CRs (UHECRs) eludes our grasp [5, 6, 7, 8]. Yet, for those within the
more modest energy confines of $E<0.1$ EeV, supernovae within the galactic
bounds are believed to be the celestial forges of acceleration [9, 10, 11, 12,
13]. Ascend higher, to realms around $1$ EeV and beyond, and one enters the
extragalactic domain, are supposed to be accelerated by gamma-ray bursts or
active galactic nuclei [2].

The CR energy spectrum is a vast continuum, stretching from the GeV scale to
the lofty heights of $100$ EeV, tracing a power-law spectrum. Along this
spectrum, a subtle inflection known as the ‘knee’ emerges at roughly $4\times
10^{15}$ eV, and a ‘second knee’ where the spectrum steepens more at around
$10^{17}$ eV [14, 15], followed by a flattening at the ‘ankle’ near $5\times
10^{18}$ eV. And then, near $5\times 10^{19}$ eV, one obtains the Greisen-
Zatsepin-Kuzmin (GZK) cutoff—a stark demarcation in the CR spectrum [16, 17,
18, 19]. CRs are the silent voyagers of space, carrying energies that dwarf
those produced by any human-made accelerator, a fact that has intrigued and
challenged physicists [20]. These particles, arriving at a rate of about one
per square kilometer per century with energies exceeding $6\times 10^{19}$ eV
[20], are the messengers from the unknown, holding secrets to questions that
have persisted for more than a century. The pursuit to understand these cosmic
wanderers has led to significant advancements in experimental physics,
revealing a surprising decrease in the flux of CRs above $4\times 10^{19}$ eV
[21], contrary to earlier predictions. This suppression, confirmed by multiple
studies [20, 22, 21, 23], suggests a transition from lighter to heavier
particles as energies increase, hinting at profound processes occurring over
cosmological distances. Such observations, including the setting of upper
limits on the presence of photons [20, 24, 25], neutrinos [20, 26, 27], and
neutrons [20, 28] among these UHECRs, have been crucial in our quest to
understand the universe. The observations of a shift in composition at
energies beyond $3\times 10^{18}$ eV [29, 30] and a suppression in the flux
above $4\times 10^{19}$ eV, which could be indicative of the GZK cutoff—a
theoretical prediction of energy loss over vast intergalactic travels, were
made nearly sixty years ago [18, 19]. However, despite these insights, the
complete story of CRs remains elusive, as current data is insufficient to
confirm if energy loss during propagation is the sole cause of this
suppression.

The vast intergalactic stage is set with turbulent magnetic fields (TMFs),
playing a pivotal role in the propagation of UHECRs. These charged particles,
propelled from distant galaxies, find their paths intricately influenced by
the TMFs they encounter. The path of a UHECR through space is a delicate
balance between its traversed distance and the scattering length
$l_{\text{D}}=3D/c$, where $D$ represents the diffusion coefficient and $c$ is
the speed of light [31]. When the cosmic path is shorter than $l_{\text{D}}$,
the particle’s journey is ballistic. On the other hand, if the distance
greatly exceeds $l_{\text{D}}$, the particle’s movement becomes diffusive. The
inclusion of extragalactic TMFs and the finite density of cosmic sources can
lead to a magnetic horizon effect due to which the CR flux is suppressed for
decreasing energies [32]. This effect may reconcile observations with a higher
spectral index [11, 33, 34], aligning with the theoretical values derived from
diffusive shock acceleration. Another intriguing possibility involves the
acceleration of heavy nuclei by extragalactic sources. These nuclei, upon
interacting with the ambient infrared radiation (IR), undergo
photodisintegration. This process releases a cascade of secondary nucleons,
potentially explaining the lighter composition observed below the cosmic ankle
[35, 36]. Intergalactic magnetic fields exert a profound influence on CRs, not
only altering their trajectories but also potentially shaping their energy
distribution. This latter effect, known as the magnetic horizon effect,
manifests when magnetic fields are sufficiently intense to bar lower-energy
CRs from reaching Earth, thereby reshaping the energy spectrum [37, 38, 39,
40]. The extent of this influence hinges on the magnetic field’s strength and
the average spacing between CR sources. The propagation theorem [41] posits
that magnetic fields do not impact the particle spectrum if source distances
are significantly less than both the diffusion and attenuation lengths. Yet,
in scenarios where magnetic fields are substantial and CR sources are sparsely
scattered, a notable suppression in the low energy end of the spectrum is
anticipated. This suppression is previously studied in Ref. [33] at the low
energy case. In Ref. [42], the authors have studied the suppression by both
analytic and numerical approaches for the primary and secondary nuclei.

Albert Einstein’s General Relativity (GR), formulated in 1915, stands as a
pinnacle of theoretical physics, elegantly describing gravitational phenomena.
A century later, the LIGO collaboration confirmed Einstein’s prediction of
Gravitational Waves (GWs) [43], and the Event Horizon Telescope captured the
first image of a black hole [44, 45, 46, 47, 48, 49], both are monumental
affirmations of GR. However, GR’s inability to reconcile with quantum
mechanics and account for the universe’s accelerated expansion [50, 51, 52,
53, 54] with dark energy [55, 56, 57, 58, 59] and dark matter in galactic
rotations [60, 61, 62, 63, 64] has led to the exploration of Modified Theories
of Gravity (MTGs) [65, 66]. Among these, the $f(R)$ gravity theory [67], which
replaces the Ricci scalar $R$ in the Einstein-Hilbert action by a function
$f(R)$, has gained traction. Models like the Starobinsky [68, 69], Hu-Sawicki
[70], Tsujikawa [71], and power-law [72, 73] models of $f(R)$ gravity have
been proposed to address these cosmic conundrums. Similarly, the symmetric
teleparallel gravity and its extension, the $f(Q)$ gravity are modifications
of the standard teleparallel gravity theory, where the non-metricity scalar
$Q$ and its arbitrary function respectively replace the torsion scalar [74,
75, 76, 77, 78]. The $f(Q)$ gravity has been used to investigate the
cosmological implications, including the behavior of dark energy, and has
shown the potential to address cosmological tensions. Given MTGs’ significant
role in recent cosmological [79, 80] and astrophysical research [81, 84, 85,
86, 82, 83], their application to UHECR studies, particularly in understanding
the UHECR flux suppression, is a promising frontier. Till now various research
groups have studied the anisotropy [87, 88, 89, 90, 91, 92, 93, 94, 95, 96,
97], propagation mechanism [99, 100, 41, 102, 103, 104, 101, 98], and
suppression [33, 42, 105, 106] of UHECRs from the standard cosmology as well
as the observatories [107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117,
118]. In our previous studies, we investigate the effects of $f(R)$ gravity on
the UHECRs’ propagation [119] and anisotropy [120] for a single source system.
In this work, we extend our previous studies of UHECRs for the diffusive
suppression of UHECR flux for an ensemble of sources and some scenarios of
mixed composition of nuclei including primary, secondary, and their
intermediate nuclei, focusing on the $f(R)$ and $f(Q)$ gravity models for the
very first time and compares their implications on UHECRs’ suppression with
the result from standard $\Lambda$CDM model.

The structure of this paper is organized as follows: In Section II, we delve
into the intricate process of CR diffusion in the context of a turbulent
magnetic field. Section III is neatly divided into two distinct parts. The
first part introduces the $f(R)$ power-law model, while the other part
presents a $f(Q)$ gravity model. In Section IV, we explore the suppression
factor of the CR flux, leaving room for calculations involving multiple
sources and a diverse composition of nuclei. The numerical results, along with
their corresponding analytical fittings, are discussed in Section V. This is
further divided into Subsection V.1 and Subsection V.2, which deal with $f(R)$
and $f(Q)$ respectively. Finally, we wrap up our paper with a conclusive
remark and a productive discussion in Section VI.

##  II Diffusive Cosmic Rays in Turbulent Magnetic Fields

Creating a model for extragalactic magnetic fields is difficult due to limited
observations [121]. The strength of these fields is not precisely known, which
also varies across different regions [122, 123]. Near the cluster centers,
field strengths range from several to tens of $\mu\text{G}$ [121]. In less
dense areas, they are weaker, typically between $1$ and $10$ nG, suggesting
that significant fields may exist along the cosmic structures like filaments.
Magnetic fields correlate with matter density, being stronger in dense regions
such as superclusters and weaker in voids, possibly less than $10^{-15}$ G.
Moreover, the magnetic fields corelate with each other up to a maximum
distance, known as the coherence length $l_{\text{c}}$. Estimates suggest that
in our Local Supercluster, magnetic fields have coherence lengths from $10$
kpc to $1$ Mpc and RMS strengths between $1$ and $100$ nG [100]. The galactic
magnetic field (GMF), with strengths of a few $\mu\text{G}$, influences CRs’
travel paths but has a limited effect on their spectrum due to its smaller
scale. In our Local Supercluster, the presence of strong magnetic fields are
indicated by the observed polarised rotations of background sources, with
estimated strengths of $0.3$ to $2$ $\mu\text{G}$ [3]. These fields are
significant for CRs arriving from nearby sources. We simplify our study by
focusing on CRs’ movement through a uniform, turbulent extragalactic magnetic
field. This field is characterized by its RMS strength $B$ and coherence
length $l_{\text{c}}$, with $B$ defined as $\sqrt{\langle B^{2}(x)\rangle}$
and ranging from $1$ nG to $100$ nG [124, 125, 126], while $l_{\text{c}}$
spans from $0.01$ Mpc to $1$ Mpc [127]. The magnetic field within the Local
Supercluster is crucial for understanding CRs’ arrival on Earth, as it is most
influential for CRs from nearby sources. Therefore, we exclude the effects of
larger structures like filaments and voids.

A critical energy $E_{\text{c}}$ can be defined as the energy where a charged
particle’s Larmor radius is equal to $l_{\text{c}}$, where the Larmor radius
is given by

$r_{\text{L}}=E/Z|e|B\simeq 1.1\frac{E/\text{EeV}}{ZB/\text{nG}}\,\text{Mpc,}$

and hence the critical energy $E_{\text{c}}$ is given by
$E_{\text{c}}=Z|e|B\,l_{\text{c}}=0.9Z\left({B/\text{nG}}\right)\left({l_{\text{c}}/\text{Mpc}}\right)\text{EeV}$.
$E_{\text{c}}$ distinguishes between two diffusion regimes: below
$E_{\text{c}}$, resonant diffusion occurs, and above it, particles experience
small deflections over $l_{\text{c}}$, leading to diffusion over much longer
distances. When $E<E_{\text{c}}$, the Larmor radius $r_{\text{L}}$ is smaller
than the coherence length $l_{\text{c}}$, leading to diffusion through
resonant scattering at wavelengths comparable to $r_{\text{L}}$. This is known
as the quasi-linear regime. The diffusion length $l_{\text{D}}$ is influenced
by the magnetic field’s strength at the Larmor radius scale, which is
dependent on the turbulent magnetic field spectrum. Generally,
$l_{\text{D}}\approx l_{\text{c}}\left(E/E_{\text{c}}\right)^{\alpha}$, where
$\alpha$ varies with the spectrum: $1/3$ for Kolmogorov, $1/2$ for Kraichnan,
and $1$ for Bohm diffusion, implying $l_{\text{D}}=r_{\text{L}}$.
$l_{\text{D}}$ is the distance over which a particle typically deflects by $1$
radian, related to the random walk’s basic step. The diffusion coefficient
$D(E)$ is $c\,l_{\text{D}}/3$, and the mean square displacement
$\langle(\Delta r)^{2}\rangle$ after time $\Delta t$ is $6D\Delta t$ [33].

For $E>E_{\text{c}}$ ($r_{\text{L}}>l_{\text{c}}$), scattering is non-
resonant, caused by several small deflections with each
$\sqrt{\delta\theta}\approx l_{\text{c}}/r_{\text{L}}$. After $N\approx
l_{\text{D}}/l_{\text{c}}$ steps, the total deflection $\Delta\theta\approx
N\delta\theta=1$ radian, yielding $l_{\text{D}}\approx
l_{\text{c}}\left(E/E_{\text{c}}\right)^{2}$, indicating a rapid increase of
diffusion length with energy. If $l_{\text{D}}\ll r_{\text{s}}$ (the source
distance), then spatial diffusion occurs; if $l_{\text{D}}>r_{\text{s}}$, the
particles propagate quasi-rectilinearly, causing minor angular diffusion. A
perfect fit to the diffusion coefficient was given by Ref. [2] as

$D(E)\simeq\frac{c\,l_{\text{c}}}{3}\left[4\left(\frac{E}{E_{\text{c}}}\right)^{2}+a_{\text{I}}\left(\frac{E}{E_{\text{c}}}\right)+a_{\text{L}}\left(\frac{E}{E_{\text{c}}}\right)^{\alpha}\right],$
(1)

and hence the diffusion length is given by

$l_{D}(E)\simeq
l_{\text{c}}\left[4\left(\frac{E}{E_{\text{c}}}\right)^{2}+a_{\text{I}}\left(\frac{E}{E_{\text{c}}}\right)+a_{\text{L}}\left(\frac{E}{E_{\text{c}}}\right)^{\alpha}\right],$
(2)

where $a_{\text{I}}$ and $a_{\text{L}}$ are two constants.
$a_{\text{I}}\approx 0.9$ and $a_{\text{L}}\approx 0.23$ according to the
Kolmogorov spectrum, whereas as per the Kraichnan spectrum,
$a_{\text{I}}\approx 0.65$ and $a_{\text{L}}\approx 0.42$ [2].

The average separation between the sources of UHECRs, denoted as
$d_{\text{s}}$, is intrinsically linked to their density $n_{\text{s}}$ and
can be expressed as $d_{\text{s}}\approx n_{\text{s}}^{-1/3}$ [33]. For a
value of $n_{\text{s}}=10^{-3}~{}\text{Mpc}^{-3}$, the separation is
approximately $d_{\text{s}}\approx 10~{}\text{Mpc}$. Consequently, for a lower
density of $n_{\text{s}}=10^{-6}~{}\text{Mpc}^{-3}$, the separation increases
to $d_{\text{s}}\approx 100~{}\text{Mpc}$.

The shaping of CR spectra for different elements is fundamentally influenced
by two main factors: their interactions with radiation backgrounds and the
impact of cosmological evolution. Energy losses in CRs are due to adiabatic
losses, represented as $dE/dt=-HE$, where $H=\dot{a}/{a}$ is the Hubble’s
parameter in terms of the scale factor $a$. These losses are prevalent across
all energy levels and also result from interaction losses from various
processes. Pair creation losses ($e^{+}e^{-}$) become significant when
interacting with Cosmic Microwave Background (CMB) photons of energy
$\epsilon\sim 10^{-3}$ eV. This requires CR Lorentz factor $\Gamma>10^{9}$,
which implies that $E>A$ EeV for nuclei with mass number $A$. For protons, the
energy loss length becomes less than the Hubble horizon
$R_{\text{H}}={c}/{H_{0}}$, when proton energies exceed $E_{\text{p}}>2$ EeV.
For heavier nuclei, due to an increased threshold and a larger pair production
cross section ($\sigma\propto Z^{2}$), the energy loss length becomes
comparable to $R_{\text{H}}$ when $E_{\text{A}}\approx A$ EeV [128]. The
process of photo-pion production becomes noticeable only at extremely high
energies, specifically when $E/A>50$ EeV. Photo-disintegration of nuclei comes
into play for $E>2A$ EeV when interacting with CMB photons, and even at lower
energies when interacting with higher energy background photons, such as those
in the IR range. However, when interacting with IR photons, the energy loss
length for photo-disintegration exceeds the Hubble horizon, except for heavier
nuclei like Iron (Fe). For these heavier nuclei, the energy loss length is
larger than $R_{\text{H}}$ only when $E<30$ EeV [128]. Generally, for energies
less than $E<Z$ EeV, the main energy losses are due to adiabatic effects,
while interaction losses are minor. Besides adiabatic losses, cosmological
influences also contribute to the increase in CMB density and temperature with
redshift. These influences may also affect the evolution of other radiation
backgrounds, the density and emissivity of CR sources, and possibly the
evolution of magnetic fields. The cosmological models also have some effects
on the propagation of UHECRs. In the next section, we will discuss the
cosmological models used for this work.

##  III Cosmological Models

In this section, we will discuss the cosmological models that will be used in
the calculations for the different parameters needed for this current work.
For this purpose, we consider one model each of two modified gravity theories:
$f(R)$ and $f(Q)$. We will also introduce the Hubble parameter $H(z)$ in this
section for the models of both modified gravity theories.

###  III.1 $\mathbf{f(R)}$ gravity model

We consider here the simplest model, the $f(R)$ gravity power-law model, whose
functional form is [73, 129, 72]

$f(R)=\beta R^{n},$ (3)

where $\beta$ and $n$ are model parameters. The parameter $\beta$ has the
dependence on the value of $n$ and other cosmological parameters such as
Hubble constant $H_{0}$, current matter density parameter $\Omega_{\text{m}0}$
and also radiation density parameter $\Omega_{\text{r}0}$ via Ricci scalar at
present time $R_{0}$ as given by [129]

$\beta=-\,\frac{3H_{0}^{2}\,\Omega_{\text{m}0}}{(n-2)R_{0}^{n}},$ (4)

where the expression for $R_{0}$ is given by [129]

$R_{0}=-\,\frac{3(3-n)^{2}H_{0}^{2}\,\Omega_{\text{m}0}}{2n\left[(n-3)\Omega_{\text{m}0}+2(n-2)\Omega_{\text{r}0}\right]}.$
(5)

By using the Palatini formalism, the Friedmann equation in $f(R)$ theory in
terms of redshift $z$ can be written as [130]

$\frac{H^{2}}{H_{0}^{2}}=\frac{3\,\Omega_{\text{m}0}(1+z)^{3}+6\,\Omega_{\text{r}0}(1+z)^{4}+\frac{f(R)}{H_{0}^{2}}}{6f^{\prime}(R)\zeta^{2}},$
(6)

where

$\zeta=1+\frac{9f^{\prime\prime}(R)}{2f^{\prime}(R)}\frac{H_{0}^{2}\,\Omega_{\text{m}0}(1+z)^{3}}{Rf^{\prime}(R)-f^{\prime}(R)}.$
(7)

Thus, from Eqs. (6) and (7), we can express the Hubble parameter as a function
of the redshift $z$ as

$H(z)=\left[-\,\frac{2nR_{0}}{3(3-n)^{2}\,\Omega_{\text{m}0}}\Bigl{\\{}(n-3)\Omega_{\text{m}0}(1+z)^{\frac{3}{n}}+2(n-2)\,\Omega_{\text{r}0}(1+z)^{\frac{n+3}{n}}\Bigl{\\}}\right]^{\frac{1}{2}}.$
(8)

In our current study, we take the value of the model parameter $n=1.4$, which
is the best-fitted value for the model as analyzed in Ref. [129]. The values
of the other cosmological parameters used in this study are taken as
$H_{0}\approx 67.4$ km s-1 Mpc-1 [131], $\Omega_{\text{m}0}\approx 0.315$
[131], and $\Omega_{\text{r}0}\approx 5.373\times 10^{-5}$ [132]. Thus the
relation between cosmological time evolution and redshift can be expressed as
[119]

$\bigg{|}\frac{dt}{dz}\bigg{|}_{f(R)}=\frac{1}{(1+z)\,H}=(1+z)^{-1}\left[-\,\frac{2nR_{0}}{3(3-n)^{2}\Omega_{\text{m}0}}\Bigl{\\{}(n-3)\Omega_{\text{m}0}(1+z)^{\frac{3}{n}}+2(n-2)\Omega_{\text{r}0}(1+z)^{\frac{n+3}{n}}\Bigl{\\}}\right]^{-\,\frac{1}{2}}\\!\\!.$
(9)

This relation will be used to calculate the CR flux for this model of $f(R)$
gravity in Section IV.

###  III.2 $\mathbf{f(Q)}$ gravity model

We consider the functional form of the $f(Q)$ gravity in this work as [133]

$f(Q)=\sigma Q,$ (10)

where $\sigma$ is the model parameter. For this $f(Q)$ gravity model, the
Hubble parameter in terms of the redshift $z$ can be written as [133]

$H(z)=H_{0}\left[(1+z)^{\frac{3\sigma+C_{1}+C_{2}}{2\sigma+C_{2}}}\left(1+\frac{C_{0}}{3\sigma+C_{1}+C_{2}}\right)-\frac{C_{0}}{3\sigma+C_{1}+C_{2}}\right],$
(11)

where $C_{0}$, $C_{1}$ and $C_{2}$ are three constant parameters, considered
as bulk viscous parameters. The best-fitted values of all these model
parameters are $\sigma=-1.03^{+0.52}_{-0.55}$, $C_{0}=1.54^{+0.83}_{-0.79}$,
$C_{1}=0.08^{+0.49}_{-0.49}$, and $C_{2}=0.66^{+0.82}_{-0.83}$ as reported in
the Ref. [133]. As in the previous case, the relation between cosmological
time evolution and redshift can be expressed for this $f(Q)$ gravity model as

$\bigg{|}\frac{dt}{dz}\bigg{|}_{f(Q)}=(H_{0}(1+z))^{-1}\left[(1+z)^{\frac{3\sigma+C_{1}+C_{2}}{2\sigma+C_{2}}}\left(1+\frac{C_{0}}{3\sigma+C_{1}+C_{2}}\right)-\frac{C_{0}}{3\sigma+C_{1}+C_{2}}\right]^{-1}.$
(12)

In Section IV we will use this relation to calculate the CR flux for this
model of $f(Q)$ gravity.

In the left panel of Fig. 1, we plot the Hubble parameters for the $f(R)$
gravity power-law model and the $f(Q)$ gravity model using Eqs. (8) and (11)
respectively considering the respective set of model parameters as mentioned
above along with the Hubble parameter for the standard $\Lambda$CDM model. We
compare the results with the observational Hubble data (OHD) obtained from
differential age (DA) and Baryon Acoustic oscillations (BAO) methods [133,
119, 129]. As one can see from this plot, both modified gravity models are
well-fitted with the observational data. However, it needs to be mentioned
that all three models predict slightly different values of Hubble constant
$H_{0}$ as $67.77$ km s-1 Mpc-1 ($\Lambda$CDM), $68.4$ km s-1 Mpc-1 ($f(R)$
power-law) and $69.0$ km s-1 Mpc-1 ($f(Q)$ model) respectively [133]. It is
clear that values for the $\Lambda$CDM and $f(R)$ power-law model are very
close to the observed value of $H_{0}$ by the Planck experiment [131]. For
completeness, we also plot the cosmological time evolution with respect to
redshift $z$ for all three models of gravity considered as above in the right
panel of Fig. 1. From this plot, one can see that at the present value of
redshift i.e. at $z=0$, there are no significant differences in the
predictions of considered cosmological models. However, as the value of $z$
increases, the differences become more pronounced, especially with the $f(R)$
gravity power-law model. Whereas at $z=2.5$, the $\Lambda$CDM and $f(R)$
power-law models tend to converge, while the $f(Q)$ model follows a flat
pattern.

Figure 1: Left: Variations of Hubble parameter $H(z)$ with redshift $z$ as
predicted by the $\Lambda$CDM model, $f(R)$ gravity power-law model, and the
$f(Q)$ gravity model in comparison with the observational Hubble data (OHD)
obtained from differential age (DA) and Baryon Acoustic oscillations (BAO)
methods [133, 119, 129]. Right: Evolution of cosmological time with $z$ as
predicted by the aforementioned models of gravity.

In the following section, we will apply these cosmological models to
understand CR flux suppression. In this context, another pertinent point to be
mentioned here is that variation of $H_{0}$ values for different models will
not affect the CR flux suppression as it is independent of $H_{0}$. This will
be clear in the following section.

##  IV The Suppression of CR Flux

In the past, the diffusion of CRs in the TMF was discussed by several
researchers [101, 139, 36, 127, 140, 134, 135, 136, 137, 138, 141, 142]. In
the expanding Universe, a detailed analysis has been performed generalising
the Syrovatskii solution [143] for a proton diffusion by Berezinsky and
Gazizov [38, 101]. The proton flux anticipated from a CR source at a distance
of $r_{\text{s}}$, significantly exceeding the diffusion length $l_{D}$, can
be determined by solving the diffusion equation within the expanding Universe
[101]. The resulting expression is as follows [42]:

$J(E)=\frac{c}{4\pi}\,\int_{0}^{z_{\text{max}}}\\!\\!dz\,\biggl{|}\frac{dt}{dz}\biggl{|}\,\mathcal{L}\left[E_{\text{g}}(E,z),z\right]\frac{\textrm{exp}\left[-r_{\text{s}}^{2}/4\lambda^{2}\right]}{(4\pi\lambda^{2})^{3/2}}\,\frac{dE_{\text{g}}}{dE},$
(13)

where $z_{\text{max}}$ is the maximum redshift of the source at which it
starts to emit CRs, $E_{\text{g}}(E,z)$ is the generation energy at redshift
$z$ whose value is $E$ at $z=0$, $\mathcal{L}$ is the source emissivity that
is obtained by summing up the contributions of the charge specific emissivity
$\mathcal{L}_{\text{Z}}$ from various charges. We will employ a power-law with
a rigidity cutoff $ZE_{\text{max}}$, represented as
$\mathcal{L}_{\text{Z}}(E,z)=\xi_{\text{Z}}f(z)E^{-\gamma}/\cosh(E/ZE_{\text{max}})$
[33]. Here $\xi_{Z}$ represents the relative contribution to the CRs flux of
charge $Z$ nuclei and $f(z)$ denotes the source emissivity evolution
concerning the redshift $z$. The Styrovatskii variable $\lambda^{2}$ is given
by

$\lambda^{2}(E,z)=\int_{0}^{z}dz\,\bigg{|}\frac{dt}{dz}\bigg{|}\,(1+z)^{2}D(E_{g},z).$
(14)

Even though Eq. (13) was originally derived for protons, it can also be
applied to nuclei by expressing it in terms of particle rigidities. When it
comes to photo-disintegration processes in nuclei, they tend to conserve the
Lorentz factor and rigidity of the main fragment. As a result, these processes
do not significantly impact the particle’s diffusion properties. However,
there is a potential complication in relation to the photo-disintegration
losses. As the source term $\mathcal{L}$ should refer to the primary nucleus
that led to the observed one, obtaining this information is challenging due to
the stochastic nature of the process. In the past, S. Mollerach et al.
discussed this scenario [33] and we will extend this work in a modified
gravity framework. Since in this work, we are interested in multiple sources
instead of a single source, so according to the propagation theorem [41], to
sum all the sources, one can use

$\int_{0}^{\infty}dr\,4\pi
r^{2}\frac{\textrm{exp}\left[-r^{2}/4\lambda^{2}\right]}{(4\pi\lambda^{2})^{3/2}}=1.$
(15)

To understand how the finite distance to the sources affects suppression, we
compute the sum using a specific set of distance distributions. These
distributions correspond to a uniform source density, and we consider the
source distances from the observer as
$r_{\text{i}}=(3/4\pi)^{1/3}d_{\text{s}}\Gamma(i+1/3)/(i-1)!$ [33, 42], where
$d_{\text{s}}$ is the distance between the sources and $i$ corresponds to the
$i^{th}$ sources from a average distance. Hence, for a discrete source
distribution, performing the sum over the sources, one can get a factor [33,
42]

$F\equiv\frac{1}{n_{s}}\sum_{i}\frac{\textrm{exp}\left[-r_{\text{i}}^{2}/4\lambda^{2}\right]}{(4\pi\lambda^{2})^{3/2}}$
(16)

instead of getting Eq. (15).

In Eq. (13), after summing all the sources, we can write the combined flux for
an ensemble of sources for the $f(R)$ gravity power-law model as

$\displaystyle
J_{\text{com}}(E)\Big{|}_{f(R)}\simeq\frac{R_{\text{H}}n_{\text{s}}}{4\pi}\int_{0}^{z_{\text{max}}}dz\,(1+z)^{-1}$
$\displaystyle\left[-\,\frac{2nR_{0}}{3(3-n)^{2}\Omega_{\text{m}0}}\Bigl{\\{}(n-3)\Omega_{\text{m}0}(1+z)^{\frac{3}{n}}+2(n-2)\Omega_{\text{r}0}(1+z)^{\frac{n+3}{n}}\Bigl{\\}}\right]^{-\,\frac{1}{2}}$
$\displaystyle\times\mathcal{L}\left[E_{\text{g}}(E,z),z\right]\,\frac{dE_{\text{g}}}{dE}\,F,$
(17)

where $R_{\text{H}}=c/H_{0}$ is the Hubble radius. Similarly, for the $f(Q)$
gravity model, the combined flux can be written as

$\displaystyle
J_{\text{com}}(E)\Big{|}_{f(Q)}\simeq\frac{R_{\text{H}}n_{\text{s}}}{4\pi}\int_{0}^{z_{\text{max}}}dz\,(1+z)^{-1}$
$\displaystyle\left[(1+z)^{\frac{3\sigma+C_{1}+C_{2}}{2\sigma+C_{2}}}\left(1+\frac{C_{0}}{3\sigma+C_{1}+C_{2}}\right)-\frac{C_{0}}{3\sigma+C_{1}+C_{2}}\right]^{-1}$
$\displaystyle\times\mathcal{L}\left[E_{\text{g}}(E,z),z\right]\,\frac{dE_{\text{g}}}{dE}\,F.$
(18)

The final suppression factor of CR flux that describes the magnetic horizon
effect can be expressed as

$G(E/E_{\text{c}})\equiv\frac{J_{\text{com}}(E)}{J_{\text{com}}(E)\big{|}_{d_{\text{s}}\rightarrow
0}},$ (19)

which is the ratio of the actual flux from discrete source distribution to
that of the continuous source distribution ($d_{\text{s}}\rightarrow 0$). The
continuous source distribution corresponds to the term $F=1$ in Eqs. (IV) and
(IV), and it means that the flux is independent of the modes of propagation of
CRs. Moreover, we can rewrite the Eq. (14) in terms of Hubble radius
$R_{\text{H}}$ and from Eq. (1) as

$\lambda^{2}(E,z)=H_{0}\frac{R_{\text{H}}l_{\text{c}}}{3}\int_{0}^{z}dz\,\bigg{|}\frac{dt}{dz}\bigg{|}\,(1+z)^{2}\left[4\left(\frac{(1+z)\,E}{E_{\text{c}}}\right)^{2}+a_{\text{I}}\left(\frac{(1+z)\,E}{E_{\text{c}}}\right)+a_{\text{L}}\left(\frac{(1+z)\,E}{E_{\text{c}}}\right)^{\alpha}\right].$
(20)

The suppression factor relies on the coherence length $l_{\text{c}}$ and the
distance between the sources $d_{\text{s}}$ by this following relation [33,
144]

$X_{\text{s}}=\frac{d_{\text{s}}}{\sqrt{R_{\text{H}}l_{\text{c}}}}.$ (21)

This relation $X_{\text{s}}$ is known as the finite density factor and it will
appear in the factor $F$ given in Eq. (16) after using Eq. (20) in it. Hence
$X_{\text{s}}$ will be used in Eq. (IV) and Eq. (IV) while running the
numerical calculations. In the next section, we will discuss the results
obtained from the computations of Eq. (19) for different scenarios.

##  V Numerical Analysis

This section is completely devoted to the numerical simulation and their
corresponding analytic fitting. It is divided into two subsections for a
systematic analysis: the first subsection for the $f(R)$ gravity model and the
latter one for the model of $f(Q)$ gravity. We extensively use the python
scipy library [145] for the calculation with some assistance from the CRPropa
3.2 [146] and SimProp v2r4 [147] codes. Unless we specify, we use the RMS
value of the magnetic field as $1$ nG, the coherence length of $1$ Mpc, and
spectral index $\gamma=2$ in all of the following plots.

###  V.1 Results from $\mathbf{f(R)}$ gravity model

In the left panel of Fig. 2, we plot the suppression factor $G$ vs
$E/E_{\text{c}}$ that is obtained by the magnetic horizon effect considering
the maximum redshift of $z_{\text{max}}=2$. Here, we consider the finite
density factor $X_{\text{s}}=0.5,1,2$ and $3$. As this factor is increasing,
the departure of the suppression factor is clearly seen. The filled dots
represent the numerical simulation of the $f(R)$ model, while the hollow dots
are for the $\Lambda$CDM model. The difference between the $f(R)$ model and
the standard $\Lambda$CDM model is visible in this plot. The flux is less
comparatively suppressed in the $f(R)$ gravity model. The solid lines
represent the corresponding fits for the $f(R)$ gravity model, while the
dotted lines for the $\Lambda$CDM model. We adopt the method of Ref. [42] for
the fitting formalism, but our parameterizations are a little bit different.
The fitting equation that is used in the whole work is given by [42]

Figure 2: Left: Suppression factors for different finite density factors of
sources $X_{\text{s}}$ as predicted by the $\Lambda$CDM model and the $f(R)$
gravity power-law model. The dotted and solid lines are for the analytic
fitting for these two cosmological models respectively. Right: Same as the
left panel but for a small scaling region. Both these plots are for the no
cosmological evolution (NE) scenario. In the plots the hollow symbol
represents the $\Lambda$CDM model and the filled symbol represents the $f(R)$
gravity power-law model. This representation will be followed in the rest of
similar figures.

$P(x)=\text{exp}\left[-\left(\frac{aX_{\text{s}}}{x+b(x/a)^{\eta}}\right)^{\vartheta}\right],$
(22)

where $a$, $b$, $\vartheta$ and $\eta$ are fitting parameters. The
parameterizations required for the fitting the expression throughout this work
are given in Table 1. In the right panel of Fig. 2, we plot the same results
as the left panel, but with a zoomed view in a small-scale region. It is seen
that although the suppressions seem the same between the $E/E_{\text{c}}$
ranges of $1$ to $100$ in the left panel, in the zoomed view of the same on
the right panel, they are not exactly the same.

Figure 3: Suppression factors for two different finite density factors of
sources $X_{\text{s}}=1$ (red points) and $X_{\text{s}}=2$ (black points) as
predicted by the $\Lambda$CDM model and the $f(R)$ gravity power-law model.
The hollow and filled points represent the results from the $\Lambda$CDM model
and the $f(R)$ gravity power-law model respectively. The circle and triangle
in the plot represent the NE and SFR scenarios respectively. The solid and
dotted lines represent the analytic fitting for the said cosmological models.

The suppression factor generally depends on the luminosity evolution of the
sources with respect to the redshift. In Fig. 2, we have adopted the constant
luminosity source evolution i.e. no cosmological evolution (NE). For the
cosmological evolution star formation rate (SFR) we adopt the formalism as
follows [148]:

$\displaystyle\mathcal{L}\propto\left\\{\begin{array}[]{ll}(1+z)^{3.44},&\text{if
}z\leq 0.97\\\ (1+z)^{-0.26},&\text{if }0.97<z\leq 4.48\\\
(1+z)^{-7.8},&\text{if }z>4.48\end{array}\right.$ (26)

For this SFR case, we simulate the suppression factor upto $z_{\text{max}}=4$.
In Fig. 3, we plot the SFR scenario (filled and hollow triangle) along with
the non-evolution (filled and hollow circle) case. Here, the hollow points
represent the results from the $\Lambda$CDM model and those of filled ones
represent the $f(R)$ gravity power-law model. For a certain energy, the cosmic
particles arrive from a higher redshift in the SFR case than that of the NE
case [42]. Here also the $f(R)$ gravity model is less suppressed in the lower
energy region as well as in the higher energy range. The middle region
suppression is similar (although not exactly the same as we have already shown
in the right panel of Fig. 2) for the both $f(R)$ model and the standard
$\Lambda$CDM model. The solid and the dashed lines represent the fitting of
Eq. (22) for the $f(R)$ and $\Lambda$CDM models respectively. The fitting
parameters for the $f(R)$ power-law model and $\Lambda$CDM model each with the
NE and SFR cases are given in Table 1.

Table 1: Parametrizations of Eq. (22) for the $\Lambda$CDM model and the $f(R)$ gravity power-law model with NE and SFR scenarios. NE |  |
---|---|---
|  |  $\Lambda$CDM |  |  |  |  $f(R)$ model |  |  |
| a | b | $\vartheta$ | $\eta$ | a | b | $\vartheta$ | $\eta$ |
| 0.224 | 0.162 | 1.541 | 0.140 | 0.305 | 0.262 | 1.55 | 0.210 |
SFR |  |
|  |  $\Lambda$CDM |  |  |  |  $f(R)$ model |  |  |
| a | b | $\vartheta$ | $\eta$ | a | b | $\vartheta$ | $\eta$ |
| 0.215 | 0.334 | 1.810 | 0.270 | 0.202 | 0.163 | 1.800 | 0.190 |

As we can see from Figs. 2 and 3 a change in suppression at higher energies is
also evident and this suppression is mainly due to the radiation backgrounds.
When the energy exceeds $60$ EeV, the attenuation length of protons decreases
rapidly due to the photopion production with the CMB, falling significantly
below $100$ Mpc for energies greater than $100$ EeV [42]. A similar phenomenon
occurs for heavier nuclei at energies exceeding approximately $5Z$ EeV, as a
result of the photodisintegrations with the CMB photons [42]. In simple words,
when the energy reaches a point where the attenuation length is about the same
or less than the distance to the nearest sources, we see a significant
decrease in the flux. This scenario is seen clearly in Fig. 4. For this
computation, we take the different primaries from $Z=1$ upto $Z=26$. We
consider $Z=7$ for the CNO group, $Z=14$ for the Si group, and $Z=26$ for the
Fe group. The left panel of Fig. 4 is for the NE scenario, while the right
panel is for the SFR by considering $X_{\text{s}}=1$. It is to be noted that
the cosmological model effect is more evident in SFR than in the NE case. Here
also the same results are obtained as in the earlier case that the $f(R)$
model exhibits a lower suppression in both high and low energy regimes than
the $\Lambda$CDM model.

Figure 4: Left: Variations of suppression factor with respect to the energy
$E$ for the different primaries as predicted by the $\Lambda$CDM model and the
$f(R)$ gravity model with NE and $X_{\text{s}}=1$. Right: Same as the left
pane but for the SFR case.

In Fig. 5, we also plot the suppression factor with respect to
$E/E_{\text{c}}$ for different nuclei. In the left panel, we plot the NE case,
while the SFR case is plotted in the right panel. We can see that all nuclei
show the same behaviour for their respective primaries. Again here the nuclei
in the $f(R)$ gravity model are less suppressed in comparison with the
$\Lambda$CDM model. In the higher energy range, a sharp drop in flux appears
which is mainly due to the photodisintegration.

Figure 5: Left: Variations of suppression factor with respect to
$E/E_{\text{c}}$ for the different primaries as predicted by the $\Lambda$CDM
model and the $f(R)$ gravity model with NE and $X_{\text{s}}=1$. Right: Same
as the left pane but for the SFR case.

In all the above scenarios, we consider the magnetic field strength $B=1$ nG.
Now, the effect of a variation of the magnetic field strength in the
suppression factor can be seen in Fig. 6. Here we consider the magnetic field
strength from $1$ nG to $15$ nG. In this figure, the left panel is for the NE
case while the right one is for the SFR case. The cosmological model effect is
the same as that of the earlier cases. At low energy levels, the predictions
of the $\Lambda$CDM model at $15$ nG closely align with those of the $f(R)$
gravity power-law model at $10$ nG. The variation of suppression factor
between magnetic fields of strengths $1$ nG to $5$ nG is quite large as
compared to that for the $5$ nG to $15$ nG. We also check (not included here)
for the higher magnetic field strength such as $50$ nG, but interestingly no a
significant variation can be seen as we see the variation for the $1$ nG to
$5$ nG.

Figure 6: Left: Variations of suppression factor with respect to energy $E$
for different strengths of the magnetic field as predicted by the $\Lambda$CDM
model and the $f(R)$ gravity power-law model with NE and $X_{\text{s}}=1$.
Right: Same as the left pane but for the SFR case.

Figure 7: Left: Variations of suppression factor with respect to the magnetic
field for different energy values as predicted by the $\Lambda$CDM and the
$f(R)$ gravity power-law model with NE and $X_{\text{s}}=1$. Right: Same as
the left pane but for the SFR case.

For a more clear picture of the variation of the suppression factor with
magnetic field strength, in Fig. 7, the flux suppression factor as a function
of the magnetic field is plotted. For this scene, we take the energy values
from $E=0.1$ EeV to $E=50$ EeV, and the NE and SFR cases are also shown in the
left and right panels of this figure respectively. The standard $\Lambda$CDM
model predicts higher suppression than the $f(R)$ gravity power-law model as
expected. The cosmological model’s effect is evident throughout the plot, but
as the TMF shifts to a higher strength, the model effect becomes more. It is
seen that the suppression decreases quickly with increasing magnetic field for
low energy particles and for substantially high energy particles it almost
remains constant up to very high magnetic field strength. One important result
to be noted is that the suppression predicted by the $f(R)$ power-law model at
$E=1$ EeV in the NE case is the same as the $\Lambda$CDM at $E=0.1$ EeV in the
SFR case.

###  V.2 Results from f(Q) gravity model

Following the same procedures as discussed above, here we discuss the
numerical results obtained from the $f(Q)$ gravity model. Fig. 8 shows the
variation of the suppression factor with respect to $E/E_{\text{c}}$ as
predicted by the $f(Q)$ gravity model in comparison with that of the
$\Lambda$CDM model. In the left panel of the figure, we have shown the NE
scenario while the SFR is shown in the right panel. In the plots, the empty
circles or triangles are for the $\Lambda$CDM model while the filled ones are
for the $f(Q)$ gravity model.

Figure 8: Left: Suppression factors with respect to $E/E_{\text{c}}$ for different density factors of sources $X_{\text{s}}$ as predicted by the $\Lambda$CDM model and the $f(R)$ gravity model. Right: The suppression factors for two different density factors of sources $X_{\text{s}}=1$ and $X_{\text{s}}=2$ as predicted by the $\Lambda$CDM model and the $f(R)$ gravity model. The round and triangle-shaped points represent the NE and SFR scenarios respectively. The hollow and filled points represent the results from the $\Lambda$CDM model and the $f(Q)$ gravity model respectively. The dotted and solid lines are for the analytic fitting for these two respective cosmological models. Table 2: Parametrizations of Eq. (22) for the $f(Q)$ gravity model with NE and SFR scenario, and with $\gamma=2$. NE |
---|---
| a | b |  $\vartheta$ |  $\eta$ |
| 0.194 | 0.172 | 1.448 | 0.160 |
SFR |
| a | b |  $\vartheta$ |  $\eta$ |
| 0.175 | 0.301 | 1.710 | 0.231 |

Unlike the previous case of MTG, here the $f(Q)$ gravity model depicts a
higher suppression than that of the $\Lambda$CDM model. In the right panel,
the red points are for the density factor of sources $X_{\text{s}}=1$ and the
black points are for $X_{\text{s}}=2$. The solid and dotted lines are for the
corresponding fitting of the numerical results of the $f(Q)$ gravity model and
the $\Lambda$CDM model. We adopt the same Eq. (22) for the analytic fitting.
The cosmological model’s effect is more evident in the SFR case as compared to
the NE case. The fitting parameters relevant to this calculation are given in
Table 2.

Figure 9: Variations of suppression factor $G$ with respect to the energy $E$
(top panels) and $E/E_{\text{c}}$ (bottom panels) for the different nuclei as
predicted by the $\Lambda$CDM model and the $f(Q)$ gravity model with NE and
$X_{\text{s}}=1$.

A mixed composition of nuclei upto iron is plotted in Fig. 9. The left panels
of Fig. 9 are for the NE scenario, while the right panels are for the SFR by
considering $X_{\text{s}}=1$. Here also the same results are obtained as in
the earlier case that the $f(Q)$ model exhibits a higher suppression in both
high and low energy regimes than the $\Lambda$CDM model. In the top panels of
Fig. 9, we plot the suppression factor with respect to the energy for
different nuclei. In the bottom panels of Fig. 9, we also plot the suppression
factor with respect to $E/E_{\text{c}}$ for the same nuclei. In the left
panel, we plot the NE case, while the SFR case is plotted in the right panel.
We can see that all nuclei show the same behaviour for their respective
primaries when plotting with respect to $E/E_{\text{c}}$. Again here the
nuclei in the $\Lambda$CDM model are less suppressed in comparison with the
$f(Q)$ model. A sharp drop in flux appears in higher energy cases due to the
photodisintegration as we have already discussed in the previous section.

Figure 10: Left: Variations of suppression factor with respect to the energy
$E$ for the different strengths of the magnetic field as predicted by the
$\Lambda$CDM and the $f(Q)$ gravity model with NE and $X_{\text{s}}=1$. Right:
Same as the left pane but for the SFR case.

A discernible variation in the suppression factor with respect to the magnetic
field strength is illustrated in Fig. 10. As in the previous case, the
magnetic field strength under consideration ranges from $1$ nG to $15$ nG. In
Fig. 10, the left panel corresponds to the NE case, while the right panel
represents the SFR case. The influence of the cosmological model is consistent
with previous cases. At lower energy levels, the $\Lambda$CDM model’s
predictions at $10$ nG align closely with those of the $f(Q)$ model at $15$
nG. We will discuss these types of scenarios in the next section. The
suppression factor exhibits a substantial variation between the magnetic field
of strengths $1$ nG and $5$ nG, compared to the range from $5$ nG to $15$ nG.
We also examined higher magnetic field strengths, such as $50$ nG (results not
included here). Interestingly, the variation observed in this range was not as
significant as the variation from $1$ nG to $5$ nG. In Fig. 11, we present the
flux suppression factor as a function of the strength of the magnetic field.
For this analysis, we consider energy values ranging from $E=0.1$ EeV to
$E=50$ EeV. The left and right panels of Fig. 11 depict the NE and SFR cases,
respectively. The $f(Q)$ gravity model predicts a higher degree of suppression
compared to the standard $\Lambda$CDM model. The influence of the cosmological
model is apparent throughout the plot. However, as the TMF shifts to a higher
strength, the model’s effect becomes more. The CR flux suppression decreases
quickly with increasing magnetic field strength for low energy particles and
it remains constant for substantially high energy particles up to very high
magnetic field strength as in the previous case.

Figure 11: Left: Variations of suppression factor with respect to the magnetic
field for the different energy values as predicted by the $\Lambda$CDM and the
$f(Q)$ gravity model with NE and $X_{\text{s}}=1$. Right: Same as the left
pane but for the SFR case.

##  VI Conclusions

We study in detail the CR flux suppression in both low and high energy
regions. We consider a model for each of two MTGs for this analysis viz. the
$f(R)$ gravity and the $f(Q)$ gravity. In the left panel of Fig. 1, we have
shown the difference of the cosmological models through the Hubble parameter
expression and in the right panel the cosmological time evolution with the
redshift. We have calculated the magnetic suppression factor for all the
cosmological models including the $\Lambda$CDM by considering the NE and SFR
scenarios for different density sources $X_{\text{s}}$ in Fig. 2 and 3. A
scenario of the nuclei flux suppression as a function of both $E$ and
$E/E_{\text{c}}$ is presented in detail in Fig. 4 and Fig. 5. We parameterize
all these numerical results through an analytical factor $P(E/E_{\text{c}})$
given in Eq. (22). We consider the effect of the magnetic field in the
suppression factor for all the cosmological models in Figs. 6 and 7. We
discuss all of the above cases for the $f(Q)$ gravity model in terms of Figs.
8 – 11 including NE and SFR scenarios. The effect of the cosmological model is
better distinguished in the SFR case as compared to the case of NE.

Figure 12: Suppression factors for the different scenarios (see text) for the
$\Lambda$CDM (dotted lines), $f(R)$ gravity power-law (dash-dot lines), and
$f(Q)$ gravity model (solid lines) in NE case.

In Fig. 12, we summarize the results of all three cosmological models for the
study of the suppression factor in the case of NE (all panels) and
$X_{\text{s}}=1$ (top right, bottom left and right panels). Since we have
observed the same kind of patterns in SFR cases also, we have not included
this scenario here. In the shaded regions, the upper boundary line denotes the
results from the $f(Q)$ gravity model (solid line), while the lower boundary
line denotes the $f(R)$ model (dash-dot line). The standard $\Lambda$CDM model
(dotted line) lies between these two. From the top left panel, we can see that
as the finite density of sources $X_{\text{s}}$ increases, the cosmological
model’s effect in the high energy regime becomes narrower. In the top right
panel, a scenario of the suppression factor is shown for the different primary
particles as predicted by the considered cosmological models. From the bottom
left panel, we can see that in between the magnetic field strength of $10$ nG
and $15$ nG, all three cosmological models exhibit a similar kind of results
in both low and high energy regimes. Again from the bottom right panel, one
can see that at the high energy regime, all considered cosmological models
behave as the same. Thus from this Fig. 12, we can conclude that the $f(R)$
gravity power-law model exhibits low suppression results, the standard
$\Lambda$CDM model exhibits a moderate suppression, and the $f(Q)$ gravity
model performs a higher suppression results in the UHECR flux.

The observations from the Pierre Auger Observatory indicate that as energy
increases, heavier elements become more dominant in the CRs composition. To
avoid the overlapping between different elements, strong suppression of heavy
elements is necessary at lower energies. While elementary spectra with a very
hard spectral index at the sources ($\gamma<1$) could explain this, it
contradicts the expectations from second-order Fermi acceleration [153, 42].
An alternative explanation discussed in this present study involves a magnetic
suppression effect, leading to a hardened spectrum for low rigidities. This
effect could account for the observed composition and spectrum patterns with
the spectral index at $\gamma=2$ [33, 154]. The MTGs play an effective role in
explaining the behaviours of CR flux suppression in the Universe with
accelerated expansion. The primary objective of this work is to investigate
the potential impact of MTGs on the suppression of CR flux. The various
predictions put forth by these models highlight the need for further
refinement in our comprehension of UHECR sources, their propagation, and the
underlying principles of gravity. Moving forward, this study can be expanded
by incorporating additional MTGs into both quantitative and more realistic
analyses including secondary particles along with the primaries, aiming to
gain deeper insights into the characteristics of CRs.

The extragalactic UHECRs that arrive from discrete sources differ from a
continuous distribution across the space. At higher energies, the attenuation
length of particles owing to interactions with the radiation backgrounds is
comparable to the separation distance between the sources, resulting in a
stronger suppression of flux. The background interactions can dramatically
reduce the UHECR flux from nearby sources. This can have a significant impact
on the CR flux recovery if the source spectrum extends beyond the attenuation
limit [42]. Due to the magnetic horizon effect, the low energy flux gets
suppressed from a discrete source distribution in the presence of a magnetic
field. When particles like protons and nuclei come from space, they can be
affected by the magnetic fields. If these particles have the same rigidity,
which is related to their charge and energy, they will follow the same
trajectory through these magnetic fields. Particles that keep their original
mass, undergo a strong effect from the magnetic field, similar to what happens
to protons. On the other hand, particles that arrive with a smaller mass than
they started with are less affected by the magnetic field. Interestingly, even
when particles lose some of their mass and charge during their journey as a
result of photodisintegration, they mostly keep the same rigidity. This is
because as their mass and charge decrease together, keeping their ratio the
same. So, these secondary particles follow a trajectory through the magnetic
field that is similar to that of the primary particles that exhibit without
photodisintegration. We will study these secondary particles generated for
different cosmological models in our next work.

## Acknowledgements

SPS is grateful to the authors of CRPropa and SimProp for making it available.
UDG is thankful to the Inter-University Centre for Astronomy and Astrophysics
(IUCAA), Pune, India for the Visiting Associateship of the institute.

## References

  * [1] V. F. Hess, _Uber Beobachtungen der durchdringenden Strahlung bei sieben Freiballonfahrten_ , Phys. Z. 13, 1084 (1912).

  * [2] D. Harari, S. Mollerach, E. Roulet, _Anisotropies of ultrahigh energy cosmic rays diffusing from extragalactic sources_ , Phys. Rev. D 89, 123001 (2014) [arXiv:1312.1366].
  * [3] S. Mollerach, E. Roulet, _Ultrahigh energy cosmic rays from a nearby extragalactic source in the diffusive regime_ , Phys. Rev. D 99, 103010 (2019) [arXiv:1903.05722].
  * [4] V. Berezinsky, A. Z. Gazizov, O. Kalashev, _Cascade photons as test of protons in UHECR_ , Astropart. Phys., 84, 52 (2016) [arXiv:1606.09293 ].
  * [5] V. Berezinsky, A. Z. Gazizov, S. I. Grigorieva, _Signatures of AGN model for UHECR_ [arXiv:astro-ph/0210095].
  * [6] M. Nagano, A. A. Watson, _Observations and implications of the ultrahigh-energy cosmic rays_ , Rev. Mod. Phys. 72, 689 (2000).
  * [7] P. Bhattacharjee, G. Sigl, _Origin and propagation of extremely high-energy cosmic rays_ , Phys. Rept. 327 (2000) [arXiv:astro-ph/9811011].
  * [8] A. V. Olinto, _Ultra High Energy Cosmic Rays: The theoretical challenge_ , Phys. Rept. 333 (2000) [arXiv:astro-ph/0002006].
  * [9] P. Blasi, _The Origin of Galactic Cosmic Rays_ , Astron. Astrophys. Rev. 21, 70 (2013), [arXiv:1311.7346].
  * [10] E. G. Berezhko, H. Volk, _Spectrum of cosmic rays, produced in supernova remnants_ , Astrophys. J. Lett. 661, L175 (2007) [arXiv:0704.1715].
  * [11] S. Mollerach, E. Roulet, _Anisotropies of ultrahigh-energy cosmic rays in a scenario with nearby sources_ , Phys. Rev. D 105, 063001 (2022) [arXiv:2111.00560].
  * [12] P. Abreu et al. (Pierre Auger Collaboration), _Constraints on the origin of cosmic rays above $10^{18}$ eV from large scale anisotropy searches in data of the Pierre Auger Observatory_, Astrophys. J. Lett. 762, (2013) L13 [ arXiv:1212.3083v1].
  * [13] R. U. Abbasi (Telescope Array Collaboration), _Search for EeV Protons of Galactic Origin_ , Astropart. Phys. 86, 21 (2017) [ arXiv:1608.06306v2].
  * [14] M. G. Aartsen et al. (IceCube Collaboration), _Cosmic ray spectrum and composition from PeV to EeV using 3 years of data from IceTop and IceCube_ , Phys. Rev. D 100, 082002 (2019).
  * [15] D. Ivanov et al. (Telescope Array Collaboration), _Energy Spectrum Measured by the Telescope Array Experiment_ , PoS ICRC 298, (2020).
  * [16] P. Abreu et al. (Pierre Auger Collaboration), _The energy spectrum of cosmic rays beyond the turn-down around $10^{17}$ eV as measured with the surface detector of the Pierre Auger Observatory_, Eur. Phys. J. C 81 (2021) 966 [arXiv:2109.13400v3].
  * [17] V. Novotný, the Pierre Auger Collaboration, _Energy spectrum of cosmic rays measured using the Pierre Auger Observatory_ , PoS ICRC 324 (2021).
  * [18] K. Greisen, _End to the Cosmic-Ray Spectrum?_ , Phys. Rev. Lett. 16, 748 (1966).
  * [19] G. T. Zatsepin, V. A. Kuzmin, _Upper Limit of the Spectrum of Cosmic Rays_ , JETP. Lett. 4, 78 (1966).
  * [20] D. Harari, _Ultra-high energy cosmic rays_ , Phys. Dark Universe 4, 23 (2014).
  * [21] J. Abraham et al. (The Pierre Auger Collaboration), _Observation of the Suppression of the Flux of Cosmic Rays above 4 $\times 10^{19}$ eV_, Phys. Rev. Lett. 101, 061101 (2008) [arXiv:0806.4302v1].
  * [22] R. U. Abbasi et al. (HiRes Collaboration), _First Observation of the Greisen-Zatsepin-Kuzmin Suppression_ , Phys. Rev. Lett. 100 (2008) 101101 [arXiv:astro-ph/0703099v2].
  * [23] T. Abu-Zayyad et al. (Telescope Array Collaboration), _The Cosmic-Ray Energy Spectrum Observed with the Surface Detector of the Telescope Array Experiment_ , Astrophys. J. 768, L1 (2013).
  * [24] J. Abraham et al., (Pierre Auger Collaboration), _Upper limit on the cosmic-ray photon fraction at EeV energies from the Pierre Auger Observatory_ , Astropart. Phys. 31, 399 (2009) [ arXiv:0903.1127v2].
  * [25] T. Abu-Zayyad et al. (Telescope Array Collaboration), _Upper limit on the flux of photons with energies above $10^{19}$ eV using the Telescope Array surface detector_, Phys. Rev. D 88, 112005 (2013) [arXiv:1304.5614].
  * [26] P. Abreu et al., (Pierre Auger Collaboration), _Search for Point-like Sources of Ultra-high Energy Neutrinos at the Pierre Auger Observatory and Improved Limit on the Diffuse Flux of Tau Neutrinos_ , Astrophys. J. 755, L4 (2012).
  * [27] M. G. Aartsen et al., (IceCube Collaboration), _Search for a diffuse flux of astrophysical muon neutrinos with the IceCube 59-string configuration_ , Phys. Rev. D 89, 062007 (2014) [arXiv:1311.7048].
  * [28] P. Abreu et al., (Pierre Auger Collaboration), _A Search for Point Sources of EeV Neutrons_ , Astrophys. J. 760, 148 (2012).
  * [29] J. Abraham et al., (Pierre Auger Collaboration), _Measurement of the Depth of Maximum of Extensive Air Showers above $10^{18}$ eV_ , Phys. Rev. Lett. 104, 091101 (2010) [arXiv:1002.0699v1].
  * [30] Pierre Auger Collaboration, _Interpretation of the depths of maximum of extensive air showers measured by the Pierre Auger Observatory_ , JCAP 02, 026 (2013).
  * [31] A. D. Supanitsky, _Cosmic ray propagation in the Universe in presence of a random magnetic field_ , JCAP 04, 046 (2021) [arXiv:2007.09063].
  * [32] A. A. Halim et al., (Pierre Auger Collaboration), _Impact of the Magnetic Horizon on the Interpretation of the Pierre Auger Observatory Spectrum and Composition Data_ , [arXiv:2404.03533].
  * [33] S. Mollerach, E. Roulet, _Magnetic diffusion effects on the Ultra-High Energy Cosmic Ray spectrum and composition_ , JCAP 10, 013 (2013) [arXiv:1305.6519].
  * [34] D. Wittkowski for The Pierre Auger Collaboration, _Reconstructed properties of the sources of UHECR and their dependence on the extragalactic magnetic field_ , PoS 563 (ICRC2017).
  * [35] M. Unger, G. R. Farrar, L. A. Anchordoqui, _Origin of the ankle in the ultra-high energy cosmic ray spectrum and of the extragalactic protons below it_ , Phys. Rev. D 92, 123001 (2015) [arXiv:1505.02153].
  * [36] N. Globus, D. Allard, E. Parizot, _A complete model of the CR spectrum and composition across the Galactic to Extragalactic transition_ , Phys. Rev. D 92, 021302 (2015) [arXiv:1505.01377].
  * [37] M. Lemoine, _Extra-galactic magnetic fields and the second knee in the cosmic-ray spectrum_ Phys. Rev. D 71, 083007 (2005) [arXiv:astro-ph/0411173]
  * [38] V. Berezinsky and A. Z. Gazizov, _Diffusion of Cosmic Rays in the Expanding Universe. II. Energy Spectra of Ultra-High Energy Cosmic Rays_ , Astrophys. J. 669, 684 (2007) [arXiv:astro-ph/0702102].
  * [39] V. Berezinsky, _Propagation and origin of ultra high-energy cosmic rays_ , Advances in Space Research 41, 2071 (2008).
  * [40] N. Globus, D. Allard, E. Parizot, _Propagation of high-energy cosmic rays in extragalactic turbulent magnetic fields: resulting energy spectrum and composition_ , A & A. 479, 97 (2008) [arXiv:0709.1541].
  * [41] R. Aloisio, V. Berezinsky, _Diffusive Propagation of Ultra-High-Energy Cosmic Rays and the Propagation Theorem_ , Astrophys. J. 612, 900 (2004) [arXiv:astro-ph/0403095].
  * [42] J. M. González, S. Mollerach, E. Roulet, _Magnetic diffusion and interaction effects on Ultrahigh Energy Cosmic Rays: protons and nuclei_ , Phys. Rev. D 104, 063005 (2021) [arXiv:2105.08138].
  * [43] B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), _Observation of Gravitational Waves from a Binary Black Hole Merger_ , Phys. Rev. Lett. 116, 061102 (2016) [arXiv:1602.03837].
  * [44] The Event Horizon Telescope Collaboration et al., _First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole_ , Astrophys. J. Lett. 871, L1 (2019).
  * [45] The Event Horizon Telescope Collaboration et al., _First M87 Event Horizon Telescope Results. II. Array and Instrumentation_ , Astrophys. J. Lett. 875, L2 (2019).
  * [46] The Event Horizon Telescope Collaboration et al., _First M87 Event Horizon Telescope Results. III. Data Processing and Calibration_ , Astrophys. J. Lett. 875, L3 (2019).
  * [47] The Event Horizon Telescope Collaboration et al., _First M87 Event Horizon Telescope Results. IV. Imaging the Central Supermassive Black Hole_ , Astrophys. J. Lett. 875, L4 (2019).
  * [48] The Event Horizon Telescope Collaboration et al., _First M87 Event Horizon Telescope Results. V. Physical Origin of the Asymmetric Ring_ , Astrophys. J. Lett. 875, L5 (2019).
  * [49] The Event Horizon Telescope Collaboration et al., _First M87 Event Horizon Telescope Results. VI. The Shadow and Mass of the Central Black Hole_ , Astrophys. J. Lett. 875, L6 (2019).
  * [50] A. G. Reiss et al., _Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant_ , Astron. J. 116, 1009 (1998) [arXiv:astro-ph/9805201].
  * [51] S. Perlmutter et al., _Measurements of $\Omega$ and $\Lambda$ from 42 High-Redshift Supernovae_, Astrophys. J. 517, 565 (1999) [arXiv:astro-ph/9812133].
  * [52] D. N. Spergel et. al., _Three-Year Wilkinson Microwave Anisotropy Probe ( WMAP ) Observations: Implications for Cosmology_ , Astrophys. J. Suppl. S 170, 377 (2007) [ arXiv:astro-ph/0603449].
  * [53] P. Astier et. al., _The Supernova Legacy Survey: Measurement of $\Omega_{M}$, $\Omega_{\Lambda}$ and $\omega$ from the First Year Data Set_, A & A 447, 31 (2006) [arXiv:astro-ph/0510447].
  * [54] G. Cognola et al., _A Class of viable modified f(R) gravities describing inflation and the onset of accelerated expansion_ , Phys. Rev. D 77, 046009 (2008) [arXiv:0712.4017].
  * [55] E. J. Copeland, M. Sami, S. Tsujikawa, _Dynamics of dark energy_ , IJMP D 15, 1753-1936 (2006) [arXiv:hep-th/0603057].
  * [56] U. D. Goswami, H. Nandan, M. Sami, _Formation of caustics in Dirac-Born-Infeld type scalar field systems_ , Phys. Rev. D 82, 103530 (2010).
  * [57] S. D. Odintsov and V. K. Oikonomou, _Geometric Inflation and Dark Energy with Axion $F(R)$ Gravity_, Phys. Rev. D 101, 044009 (2020) [arXiv:2001.06830]
  * [58] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, _Unifying Inflation with Early and Late-time Dark Energy in $F(R)$ Gravity_, Phys. Dark Universe 29, 100602 (2020) [arXiv:1912.13128].
  * [59] S. D. Odintsov and V. K. Oikonomou, _Unification of Inflation with Dark Energy in $f(R)$ Gravity and Axion Dark Matter_, Phys. Rev. D 99, 104070 (2019) [arXiv:1905.03496].
  * [60] J. H. Oort, _The force exerted by the stellar system in the direction perpendicular to the galactic plane and some related problems_ , Bull. Astron. Inst. the Netherlands. 6, 249 (1932).
  * [61] F. Zwicky, Helv. Phys. Acta. 6, 110-127 (1933); F. Zwicky, _Republication of: The redshift of extragalactic nebulae_ , Gen. Relativ. Gravit. 41, 207224 (2009).
  * [62] F. Zwicky, _On the Masses of Nebulae and of Clusters of Nebulae_ , Astrophys. J. 86, 217-246 (1937).
  * [63] K. Garrett and G. Duda, _Dark Matter: A Primer_ , Adv. Astron., 968283 (2011) [arXiv:1006.2483].
  * [64] N. Parbin, U. D. Goswami, _Galactic rotation dynamics in a new f(R) gravity model_ ,  Eur. Phys. J. C 83, 411 (2023) [arXiv:2208.06564].
  * [65] S. Nojiri and S. D. Odintsov, _Unified cosmic history in modified gravity: from F(R) theory to Lorentz non-invariant models_ , Phys. Rept. 505, 59 (2011) [arXiv:1011.0544].
  * [66] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, _Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution_ , Phys. Rept. 692, 1 (2017) [arXiv:1705.11098].
  * [67] P. Sotiriou, V. Faraoni, _$f(R)$ theories of gravity_, Rev. Mod. Phys. 82, 451 (2010) [arXiv:0805.1726].
  * [68] A. A. Starobinsky, _Disappearing cosmological constant in $f(R)$ gravity_, JETP. Lett. 86, 157 (2007) [arXiv:0706.2041].
  * [69] A. A. Starobinsky, _A New Type of Isotropic Cosmological Models without Singularity_ , Phys. Lett. B 91, 99 (1980).
  * [70] W. Hu and I. Sawicki, _Models of $f(R)$ cosmic acceleration that evade solar system tests_, Phys. Rev. D 76, 064004 (2007) [ arXiv:0705.1158v1].
  * [71] S. Tsujikawa, _Observational signatures of $f(R)$ dark energy models that satisfy cosmological and local gravity constraints_, Phys. Rev. D 77, 023507 (2008) [arXiv:0709.1391v2].
  * [72] D. J. Gogoi, U. D. Goswami, _Gravitational Waves in $f(R)$ Gravity Power Law Model_, Indian J. Phys. 96, 637 (2022) [arXiv:1901.11277v3].
  * [73] U. D. Goswami, K Deka, _Cosmological Dynamics of f(R) Gravity Scalar Degree of Freedom in Einstein Frame_ , IJMP D 22, 13 (2013) 1350083 [arXiv.1303.5868].
  * [74] J. B. Jiménez et al., _Coincident General Relativity_ , Phys. Rev. D 98, 044048 (2018) [arXiv:1710.03116].
  * [75] T. Harko et al., _Coupling matter in modified $Q$ gravity_, Phys. Rev. D 98, 084043 (2018) [arXiv:1806.10437].
  * [76] S. Mandal, P. K. Sahoo, J. R. L. Santos, _Energy conditions in $f(Q)$ gravity_, Phys. Rev. D 102, 024057 (2020) [arXiv:2008.01563].
  * [77] S. Mandal, D. Wang, P. K. Sahoo, _Cosmography in $f(Q)$ gravity_, Phys. Rev. D 102, 124029 (2020) [arXiv:2011.00420].
  * [78] N. Frusciante, _Signatures of $f(Q)$ gravity in cosmology_, Phys. Rev. D 103, 044021 (2021) [arXiv:2101.09242].
  * [79] P. Sarmah, A. De, U. D. Goswami, _Anisotropic LRS-BI Universe with $f(Q)$ gravity theory_, Phys. Dark Universe 40, 101209 (2023) [arXiv:2303.05905].
  * [80] D. J. Gogoi and U. D. Goswami, _A new f(R) Gravity Model and properties of Gravitational Waves in it_ , Eur. Phys. J. C 80, 1101 (2020) [arXiv:2006.04011].
  * [81] J. Bora, D. J. Gogoi, U. D. Goswami, _Strange stars in f(R) gravity Palatini formalism and gravitational wave echoes from them_ , JCAP 09, 057 (2022) [arXiv:2204.05473v2]
  * [82] N. Parbin et al., _Deflection angle, quasinormal modes and optical properties of a de Sitter black hole in f(T, B) gravity_ , Phys. Dark Universe 42, 101315 (2023) [arXiv:2211.02414].
  * [83] B. Hazarika, P. Phukon, _Thermodynamic topology of black holes in f(R) gravity_ , Prog. Theor. Exp. Phys. 2024, 043E01 (2024), [arXiv:2401.16756v3].
  * [84] N. Parbin, D. J. Gogoi, U. D. Goswami, _Weak gravitational lensing and shadow cast by rotating black holes in axionic Chern–Simons theory_ , Phys. Dark Universe 41, 101265 (2023) [arXiv:2305.09157]
  * [85] R. Karmakar, D. J. Gogoi, U. D. Goswami, _Thermodynamics and shadows of GUP-corrected black holes with topological defects in Bumblebee gravity_ , Phys. Dark Universe 41, 101249 (2023) [arXiv:2303.00297]
  * [86] R. Karmakar, U. D. Goswami, _Quasinormal modes, temperatures and greybody factors of black holes in a generalized Rastall gravity theory_ , Phys. Scr. 99, 055003 (2024) [arXiv:2310.18594].
  * [87] M. Chakraborty et al., (Grapes-3 Collaboration), _Large-scale cosmic ray anisotropy measured by the GRAPES-3 experiment_ , PoS ICRC 395 (2021).
  * [88] N. Globus et al., _Cosmic ray anisotropy from large-scale structure and the effect of magnetic horizons_ , MNRAS 484, 4167 (2019).
  * [89] H. Yoshiguchi et al., _Small-Scale Clustering in the Isotropic Arrival Distribution of Ultra-High-Energy Cosmic Rays and Implications for Their Source Candidates_ , Astrophys. J. 586, 1211 (2003).
  * [90] S. Mollerach, E. Roulet, _Anisotropies of ultrahigh-energy cosmic rays in a scenario with nearby sources_ , Phys. Rev. D 105, 063001 (2022) [arXiv:2111.00560v2].
  * [91] S. Mollerach , E. Roulet, O. Taborda, _Large-scale anisotropies of extragalactic cosmic rays below the ankle_ , JCAP 12, 021 (2022) [arXiv:2207.11540v2].
  * [92] A. Aab et al. (Pierre Auger Collaboration), _Observation of a Large-scale Anisotropy in the Arrival Directions of Cosmic Rays above 8 $\times 10^{18}$_ eV, Science 357, 1266 (2017) [ arXiv:1709.07321v1].
  * [93] A. U. Abeysekara et al., _All-sky Measurement of the Anisotropy of Cosmic Rays at 10 TeV and Mapping of the Local Interstellar Magnetic Field_ , Astrophys. J 871, 96 (2019).
  * [94] P. Mertsch, M. Ahlers, _Cosmic Ray Small-Scale Anisotropies in Quasi-Linear Theory_ , JCAP 11, 048 (2019) [arXiv:1909.09052v2].
  * [95] M. Ahlers, P. Mertsch, _Origin of Small–Scale Anisotropies in Galactic Cosmic Rays_ , Prog. Part. Nucl. Phys. 94, 184 (2017) [arXiv:1612.01873v1] .
  * [96] D. Harari, S. Mollerach, E. Roulet, _Cosmic ray anisotropies from transient extragalactic sources_ , Phys. Rev. D 103, 023012 (2021) [ arXiv:2010.10629v2].
  * [97] P. Erdogdu et al., _The dipole anisotropy of the 2 Micron All-Sky Redshift Survey_ , MNRAS 368, 1515 (2006).
  * [98] A. Y. Prosekin, S. R. Kelner, F. A. Aharonian, _Transition of propagation of relativistic particles from the ballistic to the diffusion regime_ , Phys. Rev. D 92, 083003 (2015) [arXiv:1506.06594].
  * [99] S. Mollerach, E. Roulet, _Ultrahigh energy cosmic rays from a nearby extragalactic source in the diffusive regime_ , Phys. Rev. D 99, 103010 (2019) [ arXiv:1903.05722v2 ].
  * [100] G. Sigl, M. Lemoine, and P. Biermann, _Ultra-high energy cosmic ray propagation in the local supercluster_ , Astropart. Phys.10, 141 (1999).
  * [101] V. Berezinsky, A. Z. Gazizov, _Diffusion of Cosmic Rays in the Expanding Universe. I._ , Astrophys. J. 643, 8 (2006) [arXiv:astro-ph/0512090].
  * [102] D. J. Bird et al., _Detection of a Cosmic Ray with Measured Energy Well beyond the Expected Spectral Cutoff due to Cosmic Microwave Radiation_ , Astrophys. J. 441 (1995) [arXiv:astro-ph/9410067].
  * [103] V. Berezinsky, A. Z. Gazizov, S. I. Grigorieva, _Dip in UHECR spectrum as signature of proton interaction with CMB_ , Phys. Lett. B 612, 147 (2005) [arXiv:astro-ph/0502550].
  * [104] V. Berezinsky, A. Z. Gazizov, S. I. Grigorieva, _On astrophysical solution to ultrahigh energy cosmic rays_ , Phys. Rev. D 74, 043005 (2006) [arXiv:hep-ph/0204357].
  * [105] V. A. Semenov, A. V. Kravtsov, D. Caprioli _Cosmic-Ray Diffusion Suppression in Star-forming Regions Inhibits Clump Formation in Gas-rich Galaxies_ , Astrophys. J 910, 126 (2021).
  * [106] A. L. Mesquita, D. Rodgers-Lee, A. A. Vidotto, R. D. Kavanagh, _The strong suppression of galactic cosmic rays reaching AU Mic b, c, and Prox Cen b_ , MNRAS 515, 1218 (2022).
  * [107] The High Resolution Fly’s Eye Collaboration, _Monocular Measurement of the Spectrum of UHE Cosmic Rays by the FADC Detector of the HiRes Experiment_ , Astropart. Phys. 23, 157 (2005) [arXiv:astro-ph/0208301].
  * [108] M. Takeda et al., _Energy determination in the Akeno Giant Air Shower Array experiment_ , Astropart. Phys. 19, 447 (2003) [arXiv:astro-ph/0209422].
  * [109] A. Aab et al. (Pierre Auger Collaboration), _Inferences on mass composition and tests of hadronic interactions from 0.3 to 100 EeV using the water-Cherenkov detectors of the Pierre Auger Observatory_ , Phys. Rev. D 96, 122003 (2017).
  * [110] A. Aab et al. (Pierre Auger Collaboration), _Large-scale Cosmic-Ray Anisotropies above 4 EeV Measured by the Pierre Auger Observatory_ , Astrophys. J. 868, 4 (2018).
  * [111] D. Ivanov (Telescope Array Collaboration), _Energy Spectrum Measured by the Telescope Array Experiment_ , Proc. Sci. ICRC 2019, 298 (2019).
  * [112] A. Aab et al. (Pierre Auger Collaboration) _Measrement of the cosmic-ray energy spectrum above 2.5 $\times 10^{18}$ eV using the Pierre Auger observatory_, Phys. Rev. D 102, 062005 (2020) [arXiv:2008.06486].
  * [113] M. Honda et al. (Akeno Collaboration), _Inelastic cross section for p-air collisions from air shower experiments and total cross section for p-p collisions up to $\sqrt{s}=24$ TeV_, Phys. Rev. Lett. 70, 525 (1993).
  * [114] R. U. Abbasi et al. (HiRes Collaboration), _A Study of the Composition of Ultra-High-Energy Cosmic Rays Using the High-Resolution Fly’s Eye_ , Astrophys. J. 622, 910 (2005).
  * [115] The Pierre Auger Collaboration, _First Estimate of the Primary Cosmic Ray Energy Spectrum above $3$ EeV from the Pierre Auger Observatory_, Proc. 29th ICRC, August 3-10, 2005, Pune, India [astro-ph/0507150].
  * [116] The Pierre Auger Collaboration et al., _Testing effects of Lorentz invariance violation in the propagation of astroparticles with the Pierre Auger Observatory_ , JCAP 01, 023 (2022) [arXiv:2112.06773].
  * [117] D. R. Bergman et al., (Pierre Auger and Telescope Array Collaborations), _The energy spectrum of ultra-high energy cosmic rays measured at the Pierre Auger Observatory and the Telescope Array_ , Eur. Phys. J Web Conf. 283, 02003 (2023).
  * [118] P. L. Biermann, V. Souza, _Centaurus A: The Extragalactic Source of Cosmic Rays with Energies above the Knee_ , Astrophys. J. 746, 72 (2012).
  * [119] S. P. Sarmah, U. D. Goswami, _Propagation and Fluxes of Ultra High Energy Cosmic Rays in $f(R)$ Gravity Theory_, Eur. Phys. J. C 84, 419 (2024) [arXiv:2303.16678].
  * [120] S. P. Sarmah, U. D. Goswami, _Anisotropies of Diffusive Ultra-high Energy Cosmic Rays in $f(R)$ Gravity Theory_ [arXiv:2309.14361].
  * [121] J. L. Han, _Observing Interstellar and Intergalactic Magnetic Fields_ , Annu. Rev. Astron 255, 111 (2017).
  * [122] Y. Hu et al., _Turbulent Magnetic Field Amplification by the Interaction of a Shock Wave and Inhomogeneous Medium_ , Astrophys. J. 941, 133 (2022).
  * [123] U. Chadayammuri, _Turbulent magnetic fields in merging clusters: a case study of Abell 2146_ , MNRAS 512, 2 (2022).
  * [124] L. Feretti et al., _Clusters of galaxies: observational properties of the diffuse radio emission_ , Astron. Astrophys. Rev. 20, 54 (2012).
  * [125] J. P. Vallée, _A Synthesis of Fundamental Parameters of Spiral Arms, Based on Recent Observations in the Milky Way_ , New Astro. Rev. 55, 91 (2011).
  * [126] F. Vazza et al., _Simulations of extragalactic magnetic fields and of their observables_ , Class. Quantum Grav. 34, 234001 (2017).
  * [127] G. Sigl, F. Miniati and T. A. Ensslin, _Ultrahigh energy cosmic ray probes of large scale structure and magnetic fields_ , Phys. Rev. D 70, 043007 (2004) [arXiv:astro-ph/0401084]
  * [128] D. Allard, _Extragalactic propagation of ultrahigh energy cosmic-rays_ , Astropart. Phys. 39-40, 33 (2012) [arXiv.1111:3290].
  * [129] D. J. Gogoi and U. D. Goswami, _Cosmology with a new $f(R)$ gravity model in Palatini formalism_, IJMP D 31, 2250048 (2022) [arXiv:2108.01409].
  * [130] B. Santos, M. Campista, J. Santos, J. S. Alcaniz, _Cosmology with Hu-Sawicki gravity in the Palatini formalism_ , A & A 548, A31 (2012) [arXiv:1207.2478].
  * [131] N. Aghanim et al., _Planck 2018 results_ (Planck Collaboration), A & A 641, A6 (2020) [arXiv:1807.06209].
  * [132] K. Nakamura and Particle Data Group, _Review of Particle Physics_ , J. Phys. G: Nucl. Part. Phys. 37, 075021(2010).
  * [133] R. Solanki et al., _Cosmic acceleration with bulk viscosity in modified $f(Q)$ gravity_, Phys. Dark Universe, 32, 100820 (2021), [arXiv:2105.00876].
  * [134] K. Kotera and M. Lemoine, _Optical depth of the Universe to ultrahigh energy cosmic ray scattering in the magnetized large scale structure_ , Phys. Rev. D 77, 023005 (2008) [arXiv:0801.1450].
  * [135] H. Yoshiguchi et al., _Small-Scale Clustering in the Isotropic Arrival Distribution of Ultra-High-Energy Cosmic Rays and Implications for Their Source Candidates_ , Astrophys. J 586, 1211 (2003).
  * [136] M. Lemoine, _Extra-galactic magnetic fields and the second knee in the cosmic-ray spectrum_ , Phys. Rev. D 71, 083007 (2005) [arXiv:astro-ph/0411173]
  * [137] D. Hooper, S. Sarkar, A. M. Taylor, _The intergalactic propagation of ultra-high energy cosmic ray nuclei_ , Astropart. Phys. 27, 199 (2007).
  * [138] _The Intergalactic Propagation of Ultra-High Energy Cosmic Ray Nuclei: An Analytic Approach_ , Phys. Rev. D 77, 103007 (2008) [arXiv:0802.1538].
  * [139] P. Blasi and A. V. Olinto, _Magnetized local supercluster and the origin of the highest energy cosmic rays_ ,  Phys. Rev. D 59, 023001 (1999).
  * [140] T. Stanev et al., _Propagation of ultrahigh energy protons in the nearby universe_ , Phys. Rev. D 62, 093005 (2000) .
  * [141] G. Sigl, _Nonuniversal spectra of ultrahigh energy cosmic ray primaries and secondaries in a structured universe_ , Phys. Rev. D 75, 103001 (2007) [arXiv:astro-ph/0703403].
  * [142] R. Aloisio _Acceleration and propagation of ultra high energy cosmic rays_ , Prog. Theor. Exp. Phys. 2017, 12A102 (2017) [arXiv:1707.08471]
  * [143] S. I. Syrovatskii, _The Distribution of Relativistic Electrons in the Galaxy and the Spectrum of Synchrotron Radio Emission_ , Soviet Astro. 3, 22 (1959).
  * [144] R. A. Batista, G. Sigl, _Diffusion of cosmic rays at EeV energies in inhomogeneous extragalactic magnetic fields_ , JCAP 11, 031 (2014) [arXiv:1407.6150].
  * [145] P. Virtanen et al., _SciPy 1.0: fundamental algorithms for scientific computing in Python_ , Nature Methods 17, 261 (2020).
  * [146] R. A. Batista et al., _CRPropa 3.2 – an advanced framework for high-energy particle propagation in extragalactic and galactic spaces_ , JCAP 09, 035 (2022) [arXiv:2208.00107].
  * [147] R. Aloisio et al., _SimProp v2r4: Monte Carlo simulation code for UHECR propagation_ , JCAP 11, 009 (2017) [arXiv:1705.03729].
  * [148] A. M. Hopkins, J. F. Beacom, _On the Normalization of the Cosmic Star Formation History_ , Astrophys. J. 651, 142 (2006) [arXiv:astro-ph/0601463].
  * [149] A. M. Taylor, M. Ahlers and F. A. Aharonian, _Need for a local source of ultrahigh-energy cosmic-ray nuclei_ ,  Phys. Rev. D 84, 105007 (2011).
  * [150] J. H. Matthews, A.R. Bell, K.M. Blundell and A.T. Araudo, _Fornax A, Centaurus A, and other radio galaxies as sources of ultrahigh energy cosmic rays_ , MNRAS Lett. 479, L76 (2018).
  * [151] R. Guedes Lang, A. M. Taylor, M. Ahlers and V. de Souza, _Revisiting the distance to the nearest ultrahigh energy cosmic ray source: Effects of extragalactic magnetic fields_ , Phys. Rev. D 102 063012 (2020).
  * [152] M. Kachelriess, P. Serpico, _The Compton–Getting effect on ultra-high energy cosmic rays of cosmological origin_ ,  Phys. Lett. B 640, 225 (2006)
  * [153] A. Aab et al. (Pierre Auger Collaboration), _Combined fit of spectrum and composition data as measured by the Pierre Auger Observatory_ , JCAP 04, 038 (2017) [arXiv:1612.07155].
  * [154] S. Mollerach, E. Roulet, _Extragalactic cosmic rays diffusing from two populations of sources_ , Phys. Rev. D 101, 103024 (2020) [arXiv:2004.04253].