<figure><img src="image_1.png" /></figure>

# Impact of Galaxy Clusters on the Propagation of Ultrahigh-energy Cosmic Rays

Antonio Condorelli<fnref n="1" />,*, Jonathan Biteau<fnref n="1" />,<fnref n="2" />,*, and Remi Adam<fnref n="3" />

Received 2023 June 19; revised 2023 September 8; accepted 2023 September 25; published 2023 November 2

## Abstract

Galaxy clusters are the largest objects in the Universe kept together by gravity. Most of their baryonic content is made of a magnetized diffuse plasma. We investigate the impact of such a magnetized environment on the propagation of ultrahigh-energy cosmic rays (UHECRs). The intracluster medium (ICM) is described according to the self-similar assumption, in which gas density and pressure profiles are fully determined by the cluster mass and redshift. The magnetic field is scaled to the thermal components of the ICM under different assumptions. We model the propagation of UHECRs in the ICM using a modified version of the Monte Carlo code *SimProp*, where hadronic processes and diffusion in the turbulent magnetic field are implemented. We provide a universal parameterization that approximates the UHECR fluxes escaping from the environment as a function of the most relevant quantities, such as the mass of the cluster, the position of the source with respect to the center of the cluster, and the nature of the accelerated particles. We show that galaxy clusters are an opaque environment, especially for UHECR nuclei. The role of the most massive nearby clusters in the context of the emerging UHECR astronomy is finally discussed.

*Unified Astronomy Thesaurus concepts: Cosmic ray astronomy (324)*

*Supporting material: figure set*

## 1. Introduction

Even though their existence has been known for more than a century, the nature and origin of cosmic rays at the highest energies remains elusive. Observations have allowed us to explore their spectral behavior and composition in terms of atomic mass on Earth (Coleman et al. 2023), but the sources of ultrahigh-energy cosmic rays (UHECRs), i.e., cosmic rays above $10^{18}$ eV, still remain unknown.

Rapid progress in computational high-energy astrophysics has dramatically advanced the study of acceleration mechanisms in systems ranging from the jets of stellar-sized objects such as gamma-ray bursts (GRBs; Sudilovsky et al. 2013) to the large-scale shocks surrounding galaxy clusters (Norman et al. 1995; Kang et al. 1997; Ryu et al. 2003). Galaxy clusters are the largest virialized structures, with typical radii of $R_{cl} = 1–2$ Mpc and total masses of $M \simeq 10^{14}–10^{15} M_\odot$, including both baryonic and dark matter. Strong turbulent magnetic fields, with rms values of $B \simeq$ few $\mu$G, are present inside clusters, having typical coherence lengths of 5–30 kpc (Bonafede et al. 2010; Donnert et al. 2018). This implies that cosmic rays accelerated in candidate sources inside the clusters, e.g., by hypernovae or GRBs in star-forming galaxies, or in the accretion shocks, jets, and radio lobes of active galactic nuclei (AGNs), can be confined for long times within clusters. They can undergo interactions with the enhanced baryonic content of the intracluster medium (ICM), whose profile is determined by Bremsstrahlung emission in X-rays (Sarazin 1986).

Upper limits on the flux of neutrinos and gamma rays at ultrahigh energies rule out a dominant origin of UHECRs from exotic particles (Abreu et al. 2023, 2022a), which should then

originate from extragalactic astrophysical sources. An extra-galactic origin is corroborated by the observation of a dipolar anisotropy above 8 EeV (Aab et al. 2018a) and an evidenced correlation of UHECRs above 40 EeV with extragalactic objects in the nearby Universe (Aab et al. 2018b; Abreu et al. 2022b). Some of these extragalactic sources could be hosted or shadowed by clusters.

Propagation of UHECRs in a specific cluster (e.g., the Virgo cluster) has been already treated in different works (Dolag 2009; Kotera et al. 2009; Harari et al. 2016; Fang & Murase 2018). Although some of these theoretical works suggested that galaxy clusters are efficient UHECR calorimeters, some authors recently claimed excesses of UHECRs from these structures (Abbasi et al. 2021; Ding et al. 2021). Revisiting the propagation of UHECRs in galaxy clusters is thus a timely topic. In the following, we evaluate whether UHECRs can escape from such environments and how clusters should be accounted for in UHECR astronomy. We provide in particular a single parameterization of the escaping flux, which depends on the mass of the cluster and on the UHECR features, such as energy and atomic mass.

The paper is organized as follows: we introduce the relevant properties of galaxy clusters and detail the way we compute the most important macroscopic quantities for our study in Section 2; the microphysics of the propagation of UHECRs in such environments is detailed in Section 3; we present our results and discuss the impact of our assumptions in Section 4. We finally draw our conclusions in Section 5.

## 2. ICM Modeling

Clusters of galaxies and the filaments that connect them are the largest structures in the present Universe in which the gravitational force due to the over-density of matter overcomes the expansion of the Universe. Massive clusters have typical

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Université Paris-Saclay, CNRS/IN2P3, IJCLab, F-91405 Orsay, France</footnote>

<footnote n="2">²Institut Universitaire de France (IUF), Paris, France</footnote>

<footnote n="3">³Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France</footnote>
total masses of the order of $10^{15} M_{\odot}$, mostly in the form of dark matter (70%–80% of the total mass), while the baryonic matter is harbored by galaxies (few percent) and composes the hot ($T \sim 10^8$ K) and tenuous ($n_{gas} \simeq 10^{-1}-10^{-4} \text{cm}^{-3}$) gas (15%–20%) that forms the ICM (Voit 2005). To model the propagation of UHECRs in this environment, we need estimates of the gas density profile, the magnetic field profile, and the coherence length. While the gas density is well understood and routinely derived from X-ray observations, this is not the case for the magnetic field, for which only a handful of measurements—sometimes model dependent—are available in the literature (Vacca et al. 2018). From theoretical arguments, however, the magnetic field strength is often assumed to scale with the ICM thermal density or pressure.

## 2.1. Density Profile

An interesting feature of galaxy clusters is that they are self-similar objects at first order, so their physical properties can be fully described given their mass and redshift (Kaiser 1986). For instance, their universal pressure profiles (UPP) and universal density profiles (UDP) are now well constrained from observations (e.g., Arnaud et al. 2010; Pratt et al. 2022). Following Arnaud et al. (2010), we use the UPP expressed as

$$P(x) = \frac{P_0 \times P_{500}(M_{500}, z) \times f(M_{500}, z)}{(c_{500}x)^{\gamma_{\text{UPP}}} \cdot (1 + (c_{500}x)^{\alpha_{\text{UPP}}})^{\frac{\beta_{\text{UPP}}-\gamma_{\text{UPP}}}{\alpha_{\text{UPP}}}}}, \quad (1)$$

with $P_{500}(M_{500}, z)$ the self-similar normalization (Nagai et al. 2007), $f(M_{500}, z)$ a small mass-dependence correction, and where $P_0$, $c_{500}$, $\alpha_{\text{UPP}}$, $\beta_{\text{UPP}}$, and $\gamma_{\text{UPP}}$ are parameters that describe the shape of the profile as a function of the scaled radius $x = r/R_{500}$.<fnref n="4" />

Similarly, we use the UDP as measured by Pratt et al. (2022), which can be expressed as

$$n(x) = \frac{A(M_{500}, z) \times f_0}{(x/x_s)(1 + (x/x_s)^{\gamma_{\text{UDP}}})^{\frac{3\beta_{\text{UDP}}-\gamma_{\text{UDP}}}{\gamma_{\text{UDP}}}}} \quad (2)$$

The quantity $A(M_{500}, z)$ describes the normalization as a function of mass and redshift, and the parameters $f_0$, $x_s$, $\alpha_{\text{UDP}}$, $\beta_{\text{UDP}}$, and $\gamma_{\text{UDP}}$ describe the shape (see also Ghirardini et al. (2019), for another calibration of the UDP).

The gas density and pressure profiles are expected to be connected. This provides us with an alternative way to describe the thermal density of the ICM given the pressure profile. Assuming a polytropic relation between gas density and pressure, using a sample of massive nearby clusters, Ghirardini et al. (2019) measured

$$P(x) = C \times n(x)^k, \quad (3)$$

where $k = 1.19$ and $C$ is a normalization constant.

With the gas density in hand, we can derive the electron-, proton-, and helium-density profiles by scaling through the mean molecular weights $\mu_{\text{gas}} = 0.61$, $\mu_e = 1.16$, $\mu_p = 1.39$, and $\mu_{\text{He}} = 14.6$ (see, e.g., Adam et al. 2020). The proton density profile of the Coma cluster, as obtained from the best-fit model describing the ROSAT data (Briel et al. 1992), is shown in Figure 1. It is compared to the model derived from our methodology, using the mass and redshift from the MCXC catalog (Piffaretti et al. 2011). The red line denotes our

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Thermal-proton density profiles for the Coma cluster: the blue line indicates the best-fit model to the ROSAT data (Briel et al. 1992), the red line indicates the density obtained from the Ghirardini et al. (2019) polytropic relation combined with the Arnaud et al. (2010) UPP profile, the green and orange lines indicate the UDP from Pratt et al. (2022) and Ghirardini et al. (2019), respectively.</figcaption></figure>

reference model, i.e., the one obtained using the UPP profile combined with the polytropic relation. For further comparison, the UDP profiles, as calibrated by Pratt et al. (2022) and Ghirardini et al. (2019), are given in green and orange, respectively. We can observe that the main differences between the data and the models, and among the models themselves, arise in the central part of the cluster. This reflects the increased intrinsic scatter among the cluster population relative to the self-similar approximation in the cluster cores, while the consistency significantly improves at $r \in [0.2R_{500}, R_{500}]$ (see, e.g., Ghirardini et al. 2019, for details). More specifically in Figure 1, the Coma cluster is a merging system with a very flat core, thus presenting a smaller central density than that given in our mean model (we also refer to the Appendix for further examples). The impact of the choice of the reference density model on our final results is discussed in Section 4.2.

## 2.2. Magnetic Field Profile

The profile of magnetic field strength can be scaled to the thermal gas density under several assumptions. Assuming the magnetic energy density to be proportional to the thermal energy, we have

$$\langle B^2(r) \rangle = 2\mu_0 P(r)/\beta_{\text{pl}}, \quad (4)$$

with $\mu_0$ the vacuum permeability. For the plasma, we set $\beta_{\text{pl}} = 200$ following the results by Walker et al. (2017) on the Perseus cluster. By combining the central magnetic field of the Coma cluster measured by Bonafede et al. (2010) and the central pressure obtained by Planck Collaboration et al. (2013), we would instead estimate $\beta_{\text{pl}} = 77$. Alternatively, assuming that the magnetic field is frozen into the plasma and amplifies under pure adiabatic compression with magnetic flux conservation, we have

$$\langle B^2(r) \rangle = B_{\text{ref}}^2 \left( \frac{n_{\text{gas}}(r)}{n_{\text{gas}}(r_{\text{ref}})} \right)^{4/3} \quad (5)$$

<footnote n="4">⁴ The mass $M_{500}$ is defined within $R_{500}$, the radius within which the cluster density is 500 times the critical density of the Universe at the cluster redshift.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Magnetic field strength profiles for Coma cluster: the blue line (and shaded region) indicates the best-fit model (and the constrained range) as obtained from the Faraday-rotation measure (Bonafede et al. 2010), the red and orange lines indicate the model obtained when scaling the magnetic field strength to the UPP and using $\beta_{pl} = 200$ and $\beta_{pl} = 77$, respectively, the orange line indicates our model when scaling the magnetic field strength to the gas density, with the UDP from Pratt et al. (2022).</figcaption></figure>

The normalization $B_{ref}$, taken at the radius $r_{ref}$ is defined using the reference Coma cluster, for which detailed measurements are available in Bonafede et al. (2010).

In Figure 2, we compare the magnetic field profile of the Coma cluster estimated from Faraday-rotation measures (Bonafede et al. 2010) to our models. The red line indicates the profile estimated using Equation (4), with $\beta_{pl} = 200$, combined with the UPP from Arnaud et al. (2010). The orange line is based instead on $\beta_{pl} = 77$. The green line uses Equation (5) with the density estimated from the UDP calibrated by Pratt et al. (2009). We observe that despite the strong assumptions involved in our modeling, the prediction follows the measurement relatively well. This is also the case in the inner region of the cluster, where the environment is expected to play a major role in the propagation of UHECRs.

In the following work, we use as a reference the UPP in Equation (1) to derive the density through the polytropic relation in Equation (3). The reference magnetic field is derived assuming constant magnetic to thermal energy density (Equation (4)) with $\beta_{pl} = 200$. We compare the different assumptions in the Appendix for a set of clusters with different morphologies and discuss the impact of these assumptions on the propagation of UHECRs in the following section.

# 3. Propagation of UHECRs in Galaxy Clusters

## 3.1. Interactions and Diffusion in a Cluster

We compute the typical timescales of photohadronic and hadronic interactions of UHECRs in the cluster environment from a modified version of the Monte Carlo code *SimProp* (see Aloisio et al. 2012, 2015, 2016). We account for interactions with photons of the cosmic microwave background (CMB) and cosmic infrared background (CIB), as well as for hadronic interactions within the ICM.

Under the assumption of a monochromatic photon field of number density $n_\gamma$, the typical interaction rate between a relativistic atomic nucleus (A) and a low energy photon is approximately $\tau_{A\gamma}^{-1} \approx c\sigma_{A\gamma}n_\gamma$, where $\sigma_{A\gamma}$ represents the cross section of the process and c is the speed of light in vacuum. If a

more realistic spectral energy distribution for the photon field is considered and the dependence of the cross section on the energy is taken into account, the interaction rate reads as (Aloisio et al. 2013)

$$ \frac{dN_{\text{int}}}{dt} = \frac{c}{2\Gamma} \int_{\epsilon}^{\infty} \sigma_{A\gamma}(\epsilon') \epsilon' \int_{\epsilon'/2\Gamma}^{\infty} \frac{n_{\gamma}(\epsilon)}{\epsilon^2} d\epsilon \, d\epsilon', \quad (6) $$

where $\Gamma$ is the Lorentz factor of the interacting nucleus. Note that primed symbols (e.g., $\epsilon'$) refer to quantities in the nucleus rest frame, whereas unmarked symbols refer to quantities in the laboratory frame.

Though spallation processes between UHECRs and gas have a negligible impact on the extragalactic medium, their role can be substantial in the ICM considering the effective time that relativistic particles spend in this environment. The timescale of the spallation process reads as

$$ \tau_{\text{spal}} = n_{\text{ICM}} \sigma_{\text{sp}} c^{-1}, \quad (7) $$

where $n_{ICM}$ is the ICM gas density and $\sigma_{sp}$ is the cross section for proton-proton or proton-nucleus interactions. This process has been implemented in *SimProp* making use of the most recent hadronic model, Sibyll2.3d (Riehn et al. 2020), a hadronic event generator. Details on the interface between the hadronic interaction model (HIM) and the in-source version of *SimProp* can be found in Condorelli et al. (2023).

In addition to interactions, diffusion in the magnetic field has to be taken into account. In fact, charged particles populating an astrophysical environment can be confined for a long time before escaping. The timescale of the diffusion reads as $t_D = R^2/D$, where R is the radius of the environment, and D is the UHECR diffusion coefficient computed in the context of quasi-linear theory (Lee et al. 2017). The expression of the diffusion coefficient is $D \approx cr_L^{2-\delta} l_c^{\delta-1}/3$, where $r_L = E/qB$ is the particle Larmor radius, $l_c$ is the coherence length of the magnetic field, and $\delta$ is the slope of the turbulence power spectrum, while B is the strength of the turbulent magnetic field. We assume $\delta = 5/3$ as prescribed for a Kolmogorov turbulence cascade. Following Subedi (2017) and Reichherzer et al. (2022), we additionally consider the transition in the diffusion regime taking place when $r_L \gtrsim l_c$. In this energy range, the diffusion coefficient is estimated as $D = D_0(r_L/l_c)^2$, where $D_0$ is the value of the diffusion coefficient computed at the energy $E_0$ such that $r_L(E_0) = l_c$. At the highest energies, the particle propagates ballistically so that the diffusion time tends to $R/c$.

Figure 3 summarizes the typical length scales for interactions and escape in the source environment for a prototype cluster (see the caption). The interplay between length scales governs the shape of the UHECR flux as well as the nuclear composition at the escape from the cluster. The shortest length scale for protons is always dictated by diffusion; this means that some protons can escape from the environment. For nuclei (e.g., nitrogen in Figure 3), photointeraction lengths are the shortest at high rigidities for the chosen parameters of the cluster (see the caption). Clusters with larger magnetic fields also present higher target densities, which reduces the length of the hadronic interaction and makes hadronic interactions predominant at lower rigidities.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Interaction and escape lengths as a function of magnetic rigidity at the center of a prototypical galaxy cluster: photohadronic interaction times (dashed-dotted lines), spallation times (dashed lines), and diffusion times (solid lines) for protons (red) and nitrogen nuclei (green). The Hubble radius (corresponding to the age of the Universe) is shown as a long-dashed line. Length scales have been calculated assuming the following parameters: $R_{500} = 1$ Mpc, $B = 1$ µG, $l_c = 10$ kpc, $n_{ICM} = 1 \cdot 10^{-4}$ cm$^{-3}$.</figcaption></figure>

## 3.2. Implementation of the Propagation in the ICM

In order to model the UHECR transport in clusters, we have developed an extension of *SimProp*. This software has been used so far in the context of the extragalactic propagation of UHECRs (see, for instance, Aab 2017; Luce et al. 2022; Halim et al. 2023). *SimProp* implements different photodisintegration cross sections and different models for the CIB. In this work, we adopt TALYS (Koning et al. 2005; Koning & Rochman 2012) for the photodisintegration cross sections and the CIB model of Gilmore et al. (2012), which are both representative of the state of the art. *SimProp* is a monodimensional propagator. Assuming spherical symmetry, all the particles are propagated along an axis of the cluster until they reach $3 \times R_{500}$, a distance beyond which the ICM has a negligible impact with respect to the extragalactic medium.

We also consider the impact of the magnetic field on the propagation of UHECRs. Charged particles moving through a uniform magnetic field undergo an angular deflection upon traversing a distance, $l_c$, of $\approx \frac{l_c}{r_L}$. A particle of energy $E$ and charge $q = Ze$ traversing a distance $L$ suffers an overall angular deflection given by $\theta(E, Z) \approx \left(\frac{L}{l_c}\right) \cdot \theta$ (Hooper et al. 2007), which depends on the properties of the environment ($B$, $L$, and $l_c$) and of the particles ($E$, $Z$). Such deflections result in an increase in the effective propagation length, $L_{eff}$, in the ICM given by Armengaud et al. (2005):

$$ \frac{L_{\text{eff}}}{L} \approx 65 \left( \frac{E/Z}{10^{20} \text{ eV}/26} \right)^{-2} \left( \frac{L}{1 \text{ Mpc}} \right) \left( \frac{l_c}{10 \text{ kpc}} \right) \left( \frac{B}{1 \text{ µG}} \right)^2 \quad (8) $$

Knowing the properties of the cluster, it is possible to compute the effective length and therefore the effective time that a particle spends in the environment.

The propagation inside the cluster environment is determined according to the following methodology: (1) the propagation axis is divided into a given number of steps, $n_{steps} \ge A$, with A atomic mass of the injected nuclei,

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Escaping proton spectra from a cluster of $M = 10^{14} M_{\odot}$ as a function of the injection point. Positive (negative) $y$ values correspond to positions closer to the nearest (furthest) edges of the cluster. The spectra are normalized to that expected without interactions in the ICM. The vertical line shows the ankle energy.</figcaption></figure>

sufficiently large to sample the interactions; (2) UHECRs are injected at a given point in the cluster and the propagation is performed only along the chosen axis; (3) the typical length scales are dependent on the position, according to the magnetic field and gas density profiles. The probability of interaction or escape changes as a function of the radius; (4) particles are moved to the following step if the interaction probability is smaller than the escape one, otherwise they lose energy, and their byproducts are accounted for in the following steps of the propagation: (5) once a particle has reached the border, if the probability of diffusion is larger than the interaction one, this particle escapes from the cluster environment and is propagated through the extragalactic medium; (6) particles that spend a time greater than the age of the Universe in the environment are considered trapped and are no longer propagated. This is a conservative assumption: the dominant time is the minimum between the age of the cluster and the age of the oldest accelerator inside it, both smaller than the age of the Universe.

## 4. UHECR Flux Escaping the ICM

Once particles escape from the magnetized environment, it is possible to evaluate what the impact of the ICM is on the UHECR spectrum as a function of the injection point. We inject $10^4$ particles logarithmically distributed in the energy range of $10^{17}-10^{21}$ eV. The results are shown in Figure 4, where the escaping fluxes are represented as a function of rigidity. The spectra are normalized to the spectrum expected if interactions and diffusion in the ICM were neglected. One can notice how the closer the injection is to the nearest edge of the environment (at $\approx +1$ Mpc), the more the escaping flux coincides with the injection spectrum: the gas and magnetic field densities are low and the propagating particles are less affected. If instead, the UHECRs cross the center of the cluster ($y \le 0$), the flux is reduced at low energies due to the trapping by the magnetic field (the so-called magnetic horizon for extragalactic propagation, e.g., Lemoine 2005; González et al. 2021). At the highest energies, fluctuations in the transmission of ICM are an artifact of the normalization procedure, in a regime where interactions of UHECRs with the CMB are important.
<figure><img src="image_7.png" /><figcaption>Figure 5. Escaping proton spectra from a cluster of different magnetic field values at the center, taken at 1 kpc (see the legend), assuming injection at the center of the cluster. The spectra are normalized to that expected without interactions in the ICM. The vertical line shows the ankle energy.</figcaption></figure>

More massive clusters present a more intense magnetic field at the center of the cluster, which shortens the magnetic horizon of UHECRs. The impact of the cluster magnetic field on the propagation of UHECRs is illustrated in Figure 5, where the escaping fluxes are shown as a function of rigidity, in this case assuming only protons at the injection.

### 4.1. Parameterization of the UHECR Escaping Flux

We provide a parameterization of the escaping fluxes as a function of the mass of the cluster M, of the position of injection point y, and of the nature of the accelerated particles (protons or nuclei) in order to describe the escaping flux above the ankle. Four representative nuclear masses are studied: $^1$H, $^4$He, $^{14}$N, and $^{28}$Si. The contribution from iron nuclei is neglected, as few, if any, are expected from simple cosmological models that describe data from the Pierre Auger Observatory (Aab 2017; Luce et al. 2022; Halim et al. 2023).

We notice in Figure 4 that a cluster mostly affects the escaping spectrum when UHECRs cross its center. In fact, it is the place where the magnetic field is most intense and where the target density is the highest. For this reason, sources placed at $y \le 0$ would have escaping fluxes shaped by the propagation in the cluster environment, while the effect is weaker for host sources placed at $y > 0$ where the traversed magnetic field is milder. For this reason, we assume clusters to be transparent for accelerators placed at $y > 0$, while we provide a single parameterization of the transparency of the clusters for $y \le 0$. We define the transparency $f(R)$ of a given cluster as the escaping flux divided by the one expected without interactions in the ICM. We approximate the transparency as a function of rigidity R by a broken power law, with full transparency at the highest energies of

$$ \log f(R) = \begin{cases} \Gamma \log(R/\rho) & R \le \rho, \\ 0 & R \ge \rho. \end{cases} \qquad (9) $$

We note that in our equation above, the break rigidity, $\rho$, depends on the mass of the cluster, M, following to first order:

$$ \log \rho = \log \rho_0 + \xi \log(M/10^{15} M_{\odot}). \qquad (10) $$

We parameterize the low-rigidity slope, $\Gamma$, of the transparency function so that it reaches a maximum value of 2 at high

<figure><img src="image_8.png" /><figcaption>Figure 6. Transparency as a function of energy for protons, helium, and nitrogen nuclei for different cluster magnetic fields (see the legend), assuming an injection point at the center of the environment. The points show the results obtained from the simulations with errors resulting from the number of injected particles. The solid lines display the proposed parameterization. The vertical line shows the ankle energy.</figcaption></figure>

cluster masses and softens at lower masses:

$$ \Gamma = \frac{2}{1 + \left(\frac{M}{M_{\text{free}}}\right)^{-\sigma}}. \qquad (11) $$

We find that $\log(M_{\text{free}}/M_{\odot}) = 14.4 \pm 0.5$ is consistent with the transparency functions of both nuclei and protons. The parameter $\sigma$ governs the evolution of the index with cluster mass. Also, in this case, we find a common value $\sigma = 0.25 \pm 0.10$ for both nuclei and protons.

The parameters are determined by fitting the model in Equation (9) to the escaping fluxes for different positions of the sources at $y \le 0$ and for different cluster masses, considering either protons only or nuclei only. We find the best parameter values of $\log(\rho_0/V) = 20.0 \pm 0.2$ for protons and $\log(\rho_0/V) = 24.3 \pm 0.3$ for nuclei, while $\xi = 0.6 \pm 0.1$ for protons and $\xi = 1.7 \pm 0.2$ for nuclei. A comparison to simulated data is shown in Figure 6, for an injection at the center of the cluster.

The two parameters that influence the rigidity at which the transition happens, i.e., $\rho_0$ and $\xi$, are larger for nuclei than for protons. This is due to the fact that the nuclei interact more than protons in the ICM, as shown in Figure 3; for this reason, the transition to $f(R) = 1$ happens at higher rigidities for nuclei.

Our simulations show that clusters of mass $M = 10^{14} M_{\odot}$ or $M = 10^{15} M_{\odot}$, with central magnetic fields of 3 and 9 $\mu$G, respectively, are able to trap nearly all protons up to the ankle. For lower magnetic fields, the effect of the ICM on protons is quite negligible above the ankle. Similar conclusions can be drawn for nuclei; nonetheless, it is important to stress that they are fully disintegrated up to at least the ankle for clusters with central magnetic fields larger than 1 $\mu$G.

The proposed parameterization describes well the impact of the galaxy cluster on the escaping flux above the ankle. The approximation of considering the environment as transparent for sources at $y > 0$ describes well the results of the simulation for weakly magnetized clusters. For clusters with $B \ge 3$ $\mu$G, the proposed parameterization for $y > 0$ overestimates the escaping fluxes on average by 0.4 dex.
Overall, we can conclude that only a few percent of nuclei can escape from clusters with $B \ge 1$ µG up to energies of $10^{19}$ eV for He and $10^{19.5}$ eV for N. Protons are strongly suppressed as well in the most massive clusters: only 40% escape at $10^{19}$ eV for a central magnetic field of 3 µG while practically none escape at this energy for $B \ge 9$ µG. Galaxy clusters are thus hostile environments for UHECRs. The filtering is more intense for nuclei, which are fully disintegrated in the most massive clusters even in the outer regions of the environment.

## 4.2. Impact of Our Assumptions

In this investigation, many assumptions have been made. This section aims to discuss their impact.

The most impacting assumption is the parameterization of the diffusion time, which is based in this work on the scaling laws expected from diffusion theory (see Section 3.1). Using instead Equation (20) of Harari et al. (2014), which is based on Monte Carlo simulations, results in an even more opaque environment, with transparency reduced by a factor $\simeq 2$ above the ankle for protons injected in the center of the cluster.

The length of the coherence of the magnetic field in the ICM also influences our results. In this work, all the clusters are assumed to have the same coherence length, 7 kpc, based on observations of Coma. More detailed constraints on this quantity would be instrumental in determining the UHECR transparency of clusters on a case-by-case basis.

Another important topic to be discussed is the assumption of the magnetic field and gas density profiles. Instead of the reference model detailed in Section 2, we performed the same analysis as in Section 4.1 using the best-fit models shown as blue lines in Figures 1 and 2 (see also the Appendix) for three different clusters: Virgo (Planck Collaboration et al. 2016), Coma (Adam et al. 2021), and Perseus (Churazov et al. 2004). In the three examined clusters, the differences in transparency are at the maximum of the order of 1%, irrespective of their morphology; therefore, it is possible to affirm that the assumption of a UPP does not influence the main results of this work. The self-similar framework is largely driven by the cluster mass, which can be a difficult quantity to measure. The accuracy of cluster-mass estimates is thus expected to be a primary source of uncertainty. One should nonetheless note that, under our reference approach, the magnetic field strength scales to first order as $B \propto \sqrt{P_{500}} \propto M_{500}^{1/3}$, so that only an order-of-magnitude uncertainty has a strong impact on the transparency of the UHECRs illustrated in Figure 6.

We only considered interactions with CMB and CIB photons, neglecting the contribution of stellar and dust-grain emission in the cluster. Harari et al. (2016) estimated such galactic contributions to be comparable to the CIB. Neglecting the galactic emission does not affect our results. This can be understood by looking at Figure 3, where the change in the slope at around $10^{19.6}$ eV in the curves labeled *photointeraction* corresponds to the transition from lower-energy interactions with CIB to higher-energy interactions with CMB. For protons, the CMB is the only relevant field within the Hubble time. For nuclei, we investigated the impact of doubling the CIB density to model the galactic emission. The transparency changes only by a few percent. For this reason, as in Harari et al. (2016) and Hussain et al. (2021), the galactic emission is neglected.

In this work, we adopt Sibyll2.3d (Riehn et al. 2020) as HIM. The systematics related to the use of this specific HIM

cannot be explored currently because no other HIM is currently implemented in *SimProp*. This investigation should be discussed in future works.

## 5. Discussion

In this work, we develop a detailed model to explore the extent to which galaxy clusters impact the propagation of UHECRs. In particular, the modeling of the cluster environment and the use of an HIM for propagation in this environment represent novelties for studies of the propagation of UHECRs.

We work under the assumption of self-similarity. From this assumption, it is possible to derive the important quantities for the propagation of UHECRs, namely, the magnetic field and gas density profiles given the mass and the redshift of the clusters. We find that the cluster environment acts as a high-pass filter, allowing a fraction of UHE protons to escape while the UHE nuclei interact with the gas and photons present in the ICM.

This work presents some advances with respect to the previous literature. The use of software dedicated to the treatment of the cluster environment is new in UHECR physics; the conclusions of this work are in line with other works that predicted that galaxy clusters are hostile environments for UHE nuclei (Kotera et al. 2009; Harari et al. 2016; Fang & Murase 2018), while they have a weaker although non-negligible effect on the propagation of UHE protons. For example, Harari et al. (2016) suggest in their Figure 6 that, for a cluster with a central magnetic field of $B = 1$ µG and coherence length of $l_c = 10$ kpc, the environment is completely transparent for protons above the ankle energy ($>5 \times 10^{18}$ eV), while it slightly affects the escaping flux of intermediate-mass nuclei (70% transparency for carbon nuclei at 10 EeV) and heavy nuclei (9% for iron nuclei at 10 EeV). In our case, neglecting the differences in the profiles and parameterization of the diffusion coefficient, we can compare these results with those obtained for a structure of mass $M = 10^{13} M_\odot$, which corresponds to a central magnetic field of $\simeq 1$ µG. The two results are comparable for protons, while our model predicts a significantly larger depletion of intermediate nuclei: 10% transparency in our work instead of 70% in Harari et al. (2016). We confirmed through simulations that this difference arises from the different treatment of the hadronic interactions; the use of an HIM, instead of a simple analytical model, can strongly influence the propagation of the cascade, affecting both the fraction of energy lost and the fragmentation of heavy nuclei. No direct comparison can be performed with the work of Fang & Murase (2018), which shows the cumulative spectra of UHECRs at Earth that escaped from a population of sources, nor with that of Kotera et al. (2009), where the UHECR flux escaping from a single cluster is arbitrarily normalized.

The present work leads to important conclusions for the emerging field of UHECR astronomy. Two different trends can be observed in the mass composition of UHECRs measured with the Pierre Auger Observatory (Yushkov 2020): a transition from heavy to light mass composition is observed up to $10^{18.3}$ eV, while data at higher energies suggest a transition to intermediate-heavy masses. Based on our simulations, we should not observe UHE nuclei coming from the inner regions of massive galaxy clusters above the ankle energy. This includes, in particular, the Virgo cluster, the closest galaxy cluster to us ($d \simeq 16$ Mpc, $M \simeq 1.2 \times 10^{14} M_\odot$,
from Planck Collaboration et al. 2016). The Pierre Auger Collaboration indeed does not see any indication of excess in this direction, which could have an important implication as pointed out in Biteau et al. (2022). Assuming that the UHECR production rate follows the star formation rate or stellar mass of nearly half a million galaxies, Biteau (2021) found that the computed sky maps should show some excess in the direction of the Virgo cluster, not present in the observed sky maps (Abreu et al. 2022b). Our work confirms that this tension is lifted by the magnetic trapping of UHECRs in Virgo, as was already hypothesized in Biteau et al. (2022) through a more naive argument (confinement time greater than the ballistic one). The result of our work reduces the discrepancies between the arrival direction model and the data, justifying the lack of UHE nuclei in the directions of the galaxy clusters and thus suggesting interesting pathways to investigate composition anisotropies.

Another application of our work is related to the dipole observed by the Pierre Auger Observatory above 8 EeV, whose direction is qualitatively explained from the distribution of local extragalactic matter and UHECR deflections in the Galactic magnetic field (Aab et al. 2017). The strong contribution to the dipole from the Virgo cluster inferred, e.g., by Ding et al. (2021) assuming that the UHECR production rate follows the distribution of matter should be significantly lowered when accounting for magnetic trapping and shadowing in the ICM. This statement is true also for the Perseus cluster ($d \simeq 74$ Mpc, $M \simeq 5.8 \times 10^{14} M_{\odot}$, from Urban et al. 2014), in the direction of which the Telescope Array Collaboration (2021) claims an indication of excess above $5.7 \times 10^{19}$ eV. From the analysis, it cannot be excluded that the UHECRs come from the vicinity or outer shocked region of the cluster. This work tends to exclude the possibility that the Telescope Array and the Pierre Auger Observatory see UHE nuclei accelerated by a host source close to the center of the Perseus or Virgo cluster; either they have to come from the cluster outskirts or they have to be UHE protons, both primaries or secondaries due to the fragmentation of heavy nuclei in the environment surrounding the accelerator.

An interesting result of this work concerns the role of filaments of the cosmic web (Kotera & Lemoine 2008). In these regions (Carretti et al. 2022), the turbulent component (inferred to be at the ≃nG level) is weaker than the regular magnetic field (≈30 nG), both much weaker than the central magnetic fields of clusters. This means that while UHECRs are trapped in the central regions of galaxy clusters, they can escape from filaments as stated in Kim et al. (2019). If, as suggested by the authors, UHECRs are correlated with filaments connected to the Virgo cluster, they should escape from galaxies in the filaments.

A possible critical aspect beyond the scope of this work concerns the secondary production. In fact, interactions of UHECRs lead to an excess of secondaries, namely, secondary cosmic rays, neutrinos, and photons, which can escape from the environment and could be in tension with the current measurements. It should be noted that secondary protons produced by the fragmentation of heavy nuclei would remain trapped in the environment, so that they would not show up at lower energies. A natural step forward in this analysis would

concern the multimessenger connection, by taking into account the emission and propagation of photons and neutrinos in the environment. In this way, it would be possible to compare the escaping gamma rays with the possible excess observed by Fermi-LAT in the direction of the Coma cluster ($d \simeq 100$ Mpc, $M \simeq 7 \times 10^{14} M_{\odot}$, from Planck Collaboration et al. 2013), as well as to determine the expected sensitivity of upcoming gamma-ray and neutrino facilities at higher energies.

## Acknowledgments

The authors would like to thank the reviewer for constructive suggestions, which helped to improve the quality of this work. A.C. and J.B. gratefully acknowledge funding from ANR via the grant Multi-messenger probe of Cosmic Ray Origins (MICRO), ANR-20-CE92-0052. A.C. acknowledges the codevelopers of SimProp for useful feedback in the development of the software.

## Appendix

In this Appendix, we show how the models of the gas density and magnetic field profiles discussed in Section 2 match the available measurements. We use reference clusters for which magnetic field profile estimates are available in the literature (Vacca et al. 2018), using the mass and redshift taken from the MCXC catalog (Piffaretti et al. 2011). For a sufficiently large sampling of mass and dynamical state, we select the Perseus, Coma, A194, Hydra A, and A2634 clusters. For each of them, a best-fit model (or a model that was matched to the data, in the case of the magnetic field) is reported: Churazov et al. (2004), Taylor et al. (2006), and Walker et al. (2017) for Perseus, Briel et al. (1992) for Coma, Govoni et al. (2017) for A194, Vogt & Enßlin (2003), Enßlin & Vogt (2006), Laing et al. (2008), Kuchar & Ensslin (2011) for Hydra A and Schindler & Prieto (1997) for A2634. These are compared to experimental data from the Planck/ROSAT project (Eckert et al. 2012), the XCOP project (Eckert et al. 2017), and the ACCEPT database (Cavagnolo et al. 2009) whenever available and to the models outlined in Section 2. The red lines correspond to the reference model of this work: UPP and polytropic assumption for determining the thermal gas density and the constant magnetic-to-thermal energy density for the magnetic field. The other possible choices are detailed in Section 2.

In Figure 7, the best-fit models of the thermal gas density profiles are matched to X-ray data. These models are beta models (or the sum of beta models), which present a flat core by construction and thus may be overly simplistic. Magnetic field estimates are derived from the matching of cluster simulations to radio data and/or Faraday-rotation measure, both generally assuming a scaling between the magnetic field strength and the thermal gas density. The uncertainties in these estimations are often not well defined. Some of the magnetic field profiles from the literature may not be fit to the data. We conclude that our reference model is in acceptable agreement with measurements from the literature in the explored mass and redshift range, given the abovementioned caveats in the measurement uncertainties and simplifying modeling assumptions.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Perseus cluster ($M \simeq 5.8 \times 10^{14} M_{\odot}$): typical relaxed cluster with a dense core. The complete figure set (11 images) is available in the online journal.
(The complete figure set (11 images) is available.)</figcaption></figure>

## ORCID iDs

Antonio Condorelli https://orcid.org/0000-0001-5681-0086

Jonathan Biteau https://orcid.org/0000-0002-4202-8939

## References

* Aab, A., Abreu, P., Aglietta, M., et al. 2017, *Sci*, 357, 1266

* Aab, A., Abreu, P., Aglietta, M., et al. 2018a, *ApJL*, 853, L29

* Aab, A., Abreu, P., Aglietta, M., et al. 2018b, *ApJ*, 868, 4

* Aab, A. E. A. 2017, *JCAP*, 2017, 038

* Abbasi, R. U., Abu-Zayyad, T., Allen, M., et al. 2021, arXiv:2110.14827

* Abreu, P., Aglietta, M., Albury, J. M., et al. 2022a, *ApJ*, 933, 125

* Abreu, P., Aglietta, M., Albury, J. M., et al. 2022b, *ApJ*, 935, 170

* Abreu, P., Aglietta, M., Albury, J. M., et al. 2023, *PhRvL*, 130, 061001

* Adam, R., Goksu, H., Brown, S., Rudnick, L., & Ferrari, C. 2021, *A&A*, 648, A60

* Adam, R., Goksu, H., Leingartner-Goth, A., et al. 2020, *A&A*, 644, A70

* Aloisio, R., Berezinsky, V., & Grigorieva, S. 2013, *ApH*, 41, 73

* Aloisio, R., Boncioli, D., di Matteo, A., et al. 2015, arXiv:1505.01347

* Aloisio, R., Boncioli, D., di Matteo, A., et al. 2016, arXiv:1602.01239

* Aloisio, R., Boncioli, D., Grillo, A., Petrera, S., & Salamida, F. 2012, *JCAP*, 2012, 007

* Armengaud, E., Sigl, G., & Miniati, F. 2005, *PhRvD*, 72, 043009

* Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, *A&A*, 517, A92

* Biteau, J. 2021, *ApJS*, 256, 15

* Biteau, J., Marafico, S., Kerfis, Y., & Deligny, O. 2022, *ICRC (Berlin)*, 37, 1012

* Bonafede, A., Feretti, L., Murgia, M., et al. 2010, *A&A*, 513, A30

* Briel, U. G., Henry, J. P., & Boehringer, H. 1992, *A&A*, 259, L31

* Carretti, E., Vacca, V., O'Sullivan, S. P., et al. 2022, *MNRAS*, 512, 945

* Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2009, *ApJS*, 182, 12

* Churazov, E., Forman, W., Jones, C., Sunyaev, R., & Böhringer, H. 2004, *MNRAS*, 347, 29

* Coleman, A., Eser, J., Mayotte, E., et al. 2023, *Aph*, 147, 102794

* Condorelli, A., Boncioli, D., Peretti, E., & Petrera, S. 2023, *PhRvD*, 107, 083009

* Ding, C., Globus, N., & Farrar, G. R. 2021, *ApJL*, 913, L13

* Dolag, K., Kachelrieb, M., & Semikoz, D. V. 2009, *JCAP*, 01, 033

* Donnert, J., Vazza, F., Brüggen, M., & ZuHone, J. 2018, *SSRv*, 214, 214

* Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, *AN*, 338, 293

* Eckert, D., Vazza, F., Ettori, S., et al. 2012, *A&A*, 541, A57

* Enßlin, T. A., & Vogt, C. 2006, *A&A*, 453, 447

* Fang, K., & Murase, K. 2018, *NatPh*, 14, 396

* Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, *A&A*, 621, A41

* Ghirardini, V., Ettori, S., Eckert, D., & Molendi, S. 2019, *A&A*, 627, A19

* Gilmore, R. C., Somerville, R. S., Primack, J. R., & Domínguez, A. 2012, *MNRAS*, 422, 3189

* González, J. M., Mollerach, S., & Roulet, E. 2021, *PhRvD*, 104, 104

* Govoni, F., Murgia, M., Vacca, V., et al. 2017, *A&A*, 603, A122

* Halim, A. A., Abreu, P., Aglietta, M., et al. 2023, *JCAP*, 2023, 024

* Harari, D., Mollerach, S., & Roulet, E. 2014, *PhRvD*, 89, 123001

* Harari, D., Mollerach, S., & Roulet, E. 2016, *JCAP*, 2016, 010

* Hooper, D., Sarkar, S., & Taylor, A. M. 2007, *Aph*, 27, 199

* Hussain, S., Alves Batista, R., de Gouveia Dal Pino, E. M., & Dolag, K. 2021, *MNRAS*, 507, 1762

* Kaiser, N. 1986, *MNRAS*, 222, 323

* Kang, H., Rachen, J. P., & Biermann, P. L. 1997, *MNRAS*, 286, 257

* Kim, J., Ryu, D., Kang, H., Kim, S., & Rey, S.-C. 2019, *SciA*, 5, eaau8227

* Koning, A., & Rochman, D. 2012, *NDS*, 113, 2841

* Koning, A. J., Hilaire, S., & Duijvestijn, M. C. 2005, in *AIP Conf. Ser.* 769, *TALYS: Comprehensive Nuclear Reaction Modeling* (Melville, NY: AIP), 1154

* Kotera, K., Allard, D., Murase, K., et al. 2009, *ApJ*, 707, 370

* Kotera, K., & Lemoine, M. 2008, *PhRvD*, 77, 123003

* Kuchar, P., & Ensslin, T. A. 2011, *A&A*, 529, A13

* Laing, R. A., Bridle, A. H., Parma, P., & Murgia, M. 2008, *MNRAS*, 391, 521

* Lee, J., Wright, J., Bertelli, N., et al. 2017, *PhPl*, 24, 052502

* Lemoine, M. 2005, *PhRvD*, 71, 083007

* Luce, Q., Marafico, S., Biteau, J., Condorelli, A., & Deligny, O. 2022, *ApJ*, 936, 62

* Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, *ApJ*, 668, 1

* Norman, C. A., Melrose, D. B., & Achterberg, A. 1995, *ApJ*, 454, 60

* Piffaretti, R., Arnaud, M., Pratt, G. W., Pointecouteau, E., & Melin, J. B. 2011, *A&A*, 534, A109

* Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, *A&A*, 554, A140

* Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, *A&A*, 596, A101

* Pratt, G. W., Arnaud, M., Maughan, B. J., & Melin, J. B. 2022, *A&A*, 665, A24

* Pratt, G. W., Croston, J. H., Arnaud, M., & Böhringer, H. 2009, *A&A*, 498, 361

* Reichherzer, P., Becker-Tjus, J., Zweibel, E. G., Merten, L., & Pueschel, M. J. 2022, *MNRAS*, 514, 2658

* Riehn, F., Engel, R., Fedynitch, A., Gaisser, T. K., & Stanev, T. 2020, *PhRvD*, 102, 102

* Ryu, D., Kang, H., Hallman, E., & Jones, T. W. 2003, *ApJ*, 593, 599

* Sarazin, C. L. 1986, *RvMP*, 58, 1

* Schindler, S., & Prieto, M. A. 1997, *A&A*, 327, 37

* Subedi, P. 2017, *ApJ*, 837, 140

* Sudilovsky, V., Greiner, J., Rau, A., et al. 2013, *A&A*, 552, A143

* Taylor, G. B., Gugliucci, N. E., Fabian, A. C., et al. 2006, *MNRAS*, 368, 1500

* Telescope Array Collaboration, Abbasi, R. U., Abu-Zayyad, T., et al. 2021, arXiv:2110.14827

* Urban, O., Simionescu, A., Werner, N., et al. 2014, *MNRAS*, 437, 3939

* Vacca, V., Murgia, M., Govoni, F., et al. 2018, *Galax*, 6, 142

* Vogt, C., & Enßlin, T. A. 2003, *A&A*, 412, 373

* Voit, G. M. 2005, *RvMP*, 77, 207

* Walker, S. A., Hlavacek-Larrondo, J., Gendron-Marsolais, M., et al. 2017, *MNRAS*, 468, 2506

* Yushkov, A. 2020, *ICRC (Madison, WI)*, 36, 482