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

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.


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 ultra-high-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;Ryu et al. 2003;Kang et al. 1997).Galaxy clusters are the largest virialized structures, having typical radii  cl = 1 − 2 Mpc and total masses  ≃ 10 14 − 10 15  ⊙ , including both baryonic and dark matter.Strong turbulent magnetic fields, with root mean square values  ≃ few 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 radiolobes 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, whose profile is determined by Bremsstrahlung emission in X rays (Sarazin 1986).
Upper limits on the flux of neutrinos and gamma rays at ultra-high energies rule out a dominant origin of UHECRs from exotic particles (Abreu et al. 2023(Abreu et al. , 2022a)), which should then originate from extragalactic astrophysical sources.An extragalactic 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.
UHECR propagation in a specific cluster (e.g. the Virgo cluster) has been already treated in different works (Dolag et al. 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 (Ding et al. 2021;Abbasi 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 parametrization 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 UHECR propagation 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.

INTRACLUSTER MEDIUM 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 matter overdensity overcomes the expansion of the universe.Massive clusters have typical total masses of the order of 10 15  ⊙ , mostly in the form of dark matter (70 − 80 % of the total mass), while baryonic matter is harbored by galaxies (few %) and composes the hot ( ∼ 10 8 K) and tenuous ( gas ≃ 10 −1 − 10 −4 cm −3 ) gas (15 − 20 %) that forms the intracluster medium (ICM, Voit 2005).To model UHECR propagation in this environment, we need estimates of the gas density profile, of the magnetic field profile and of 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 measurementssometimes 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.

Density profile
An interesting feature of galaxy clusters is that they are self-similar objects at first order, so that 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 with  500 ( 500 , ) the self-similar normalization (Nagai et al. 2007),  ( 500 , ) a small mass-dependence correction, and where  0 ,  500 ,  UPP ,  UPP ,  UPP are parameters that describe the shape of the profile as a function of the scaled radius  = / 500 . 1 Similarly, we use the UDP as measured by Pratt et al. (2022), which can be expressed as The quantity ( 500 , ) describes the normalization as a function of mass and redshift, and the parameters  0 ,   ,  UDP ,  UDP ,  UDP describe the shape (see also Ghirardini et al. 2019, for another calibration of the UDP).
The gas density and pressure profile are expected to be connected.This provides us with an alternative way to describe the ICM thermal density 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 where  = 1.19 and  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  gas = 0.61,   = 1.16,  p = 1.39 and  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 gives our reference model, i.e., the one obtained using the UPP profile combined with the polytropic relation.For further comparison, the UDP profile, 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 Thermal proton density (cm 3 ) Coma, gas density Briel et al. (1992) UPP + polytropic relation UDP (Ghirardini et al. 2019) UDP (Pratt et al. 2022) R500 Figure 1.Thermal proton density profiles for the Coma cluster: the blue line gives the best-fit model to the ROSAT data (Briel et al. 1992), the red line gives 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 give the UDP from Pratt et al. (2022) and Ghirardini et al. (2019), respectively.
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  ∈ [0.2 500 ,  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.

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 with  0 the vacuum permeability.(Bonafede et al. 2010), the red and orange lines give the model obtained when scaling the magnetic-field strength to the UPP and using  pl = 200 and  pl = 77, respectively, the orange line give our model when scaling the magnetic-field strength to the gas density, with the UDP from Pratt et al. (2022).
The normalization  ref , taken at the radius  ref is defined using the reference Coma cluster, for which detailed measurement 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 gives the profile estimated using Equation 4, with  pl = 200, combined with the UPP from Arnaud et al. (2010).The orange line is based instead on  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 strong assumptions involved in our modeling, the prediction follows relatively well the measurement.This is also the case in the inner region of the cluster, where the environment is expected to play a major role for UHECR propagation.
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  pl = 200.We compare the different assumptions in appendix for a set of clusters with different morphologies and discuss the impact of these assumptions on UHECR propagation in the following section.

Interactions and diffusion in a cluster
We compute the typical timescales for photo-hadronic and hadronic interactions of UHECRs in the cluster environment from a modified version of the Monte Carlo code SimProp (see Aloisio et al. 2012Aloisio et al. , 2015Aloisio et al. , 2016)).We account for interactions with photons of the cosmic microwave and infrared backgrounds (CMB, CIB), as well as for hadronic interactions within the ICM.
Under the assumption of a monochromatic photon field of number density   , the typical interaction rate between a relativistic atomic nucleus () and a low energy photon is approximately  −1  ≃     , where   represents the cross section of the process and  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 (Aloisio et al. 2013): where Γ is the Lorentz factor of the interacting nucleus.Note that primed symbols (e.g. ′ ) 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 negligible impact in the extragalactic medium, their role can be substantial in the ICM considering the effective time that relativistic particles spend in this environment.The timescale for the spallation process reads: where  ICM is the ICM gas density and  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, Sibyll 2.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 diffusion timescale reads:  D =  2 /, where  is the radius of the environment , and where  is the UHECR diffusion coefficient computed in the context of quasi-linear theory (Lee et al. 2017).The expression of the diffusion coefficient is: .Interaction and escaping lengths as a function of magnetic rigidity at the center of a prototypical galaxy cluster: photo-hadronic interaction times (dashed-dot 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: the particle Larmor radius,  c is the coherence length of the magnetic field and  is the slope of the turbulence power-spectrum, while  is the strength of the turbulent magnetic field.We assume  = 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   ≳   .In this energy range, the diffusion coefficient is estimated as  =  0 (  /  ) 2 , where  0 is the value of the diffusion coefficient computed at the energy  0 such that   ( 0 ) =   .At the highest energies, the particle propagates ballistically so that the diffusion time tends to /.
Figure 3 summarizes the typical length-scales for interactions and escape in the source environment for a prototype cluster (see 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), photo-interaction lengths are the shortest at high rigidities for the chosen parameters of the cluster (see caption).Clusters with larger magnetic fields also present higher target densities, which reduces the hadronic interaction length and makes hadronic interactions predominant at lower rigidities.

Implementation of ICM propagation in SimProp
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;Halim et al. 2023;Luce et al. 2022).SimProp implements different photo-disintegration cross sections and different models for the CIB.In this work, we adopt TALYS (Koning et al. 2005;Koning & Rochman 2012) for the photo-disintegration 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×  500 , a distance beyond which the ICM has negligible impact with respect to the extragalactic medium.
We also consider the impact of the magnetic field on UHECR propagation.
Charged particle moving through a uniform magnetic field undergo an angular deflection upon traversing a distance,   , of ≃ . A particle of energy  and charge  =  traversing a distance  suffers an overall angular deflection given by •  (Hooper et al. 2007), which depends on the properties of the environment (,  and   ) and of the particles (, ).Such deflections result in an increase in the effective propagation length,  eff , in the ICM given by (Armengaud et al. 2005): 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 in a given number of steps,  steps ≥ , with  atomic mass of the injected nuclei, 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 diffusion probability 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 not propagated anymore.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.

UHECR FLUX ESCAPING THE ICM
Once particles escape from the magnetized environment, it is possible to evaluate what is the impact of the ICM on the UHECR spectrum as a function of the injection point.We inject 10 4 particles logarithmically distributed in the energy range 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 ≈ +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 ( ≤ 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 ICM transmission are an artifact of the normalization procedure, in a regime where interactions of UHECRs with the CMB are important.
More massive clusters present 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.

Parametrization of the UHECR escaping flux
We provide a parametrization of the escaping fluxes as a function of the mass of the cluster , of the position of injection point  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, 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;Halim et al. 2023;Luce et al. 2022).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  ≤ 0 would have an escaping fluxes shaped by the propagation in the cluster environment, while the effect is weaker for host sources placed at  > 0 where the traversed magnetic field is milder.For this reason, we assume clusters to be transparent for accelerators placed at  > 0, while we provide a single parametrization of the transparency of the clusters for  ≤ 0. We define the transparency  () 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  by a broken power-law, with full transparency at the highest energies: We notice that, in our equation above, the break rigidity, , depends on the mass of the cluster, , following to first order: log  = log  0 +  log(/10 15  ⊙ ). (10) We parametrize the low-rigidity slope, Γ, of the transparency function so that it reaches a maximum value of 2 at high cluster masses and softens at lower masses: We find that log( free / ⊙ ) = 14.4 ± 0.5 is consistent with the transparency functions of both nuclei and protons.The parameter  governs the evolution of the index with cluster mass.Also in this case, we find a common value  = 0.25 ± 0.10 for both nuclei and protons.
The parameters are determined by fitting the model in Equation 9to the escaping fluxes for different position of the sources at  ≤ 0 and for different cluster masses, considering either protons only or nuclei only.We find best parameter values of log( 0 /V) = 20.0 ± 0.2 for protons and log( 0 /V) = 24.3 ± 0.3 for nuclei, while  = 0.6 ± 0.1 for protons and  = 1.7 ± 0.2 for nuclei.The 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.  0 and , 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  () = 1 happens at higher rigidities for nuclei.
Our simulations show that clusters of mass  = 10 14  ⊙ or  = 10 15  ⊙ , with central magnetic fields of 3 and 9 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 G.
The proposed parametrization 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  > 0 describes well the results of the simulation for weakly magnetized clusters.For clusters with  ≥ 3 G, the proposed parametrization for  > 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  ≥ 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  ≥ 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.

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 parametrization of the diffusion time, which is based in this work on the scaling laws expected from diffusion theory (see Sec. 3.1).Using instead the equation 20 of Harari et al. (2014), which is based on Monte-Carlo simulations, results in an even more opaque environment, with a transparency reduced by a factor ≃ 2 above the ankle for protons injected in the center of the cluster.
The coherence length of the magnetic field in the ICM also influences our results.In this work, all the clus-ters 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 Figure 1 and 2 (see also 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 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  ∝ √  500 ∝  1/3 500 , so that only an orderof-magnitude uncertainty has a strong impact on the UHECR transparency illustrated in Figure 6.
We only considered interactions with CMB and CIB photons, neglecting the contribution of stellar and dustgrain 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 looking at Figure 3, where the change of slope at around 10 19.6 eV in the curves labeled "photo-interaction" 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 few percents.For this reason, as in Harari et al. (2016); 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 are currently implemented in SimProp.This investigation should be discussed in future works.

DISCUSSION
In this work, we develop a detailed model to explore the extent to which galaxy clusters impact UHECR propagation.In particular, the modelling of the cluster environment and the use of a HIM for propagation in this environment represent novelties for UHECR propagation studies.
We work under the assumption of self-similarity.From this assumption, it is possible to derive the important quantities for UHECR propagation, namely the magnetic-field and the 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 a 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 which predicted that galaxy cluster are hostile environment 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  = 1  and coherence length of  c = 10 kpc, the environment is completely transparent for protons above the ankle energy (> 5 • 10 18 eV), while it affects slightly 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 parametrization of the diffusion coefficient, we can compare these results with those obtained for a structure of mass  = 10 13  ⊙ , which corresponds to a central magnetic field of ≃ 1 .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 a 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 UHECRs 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 intermediateheavy 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 ( ≃ 16 Mpc,  ≃ 1.2•10 14  ⊙ , 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 in Biteau et al. (2021).Assuming that the UHECR production rate follows the star-formation rate or stellar mass of nearly half million of galaxies, Biteau (2021) found that the computed skymaps should show some excess in the direction of the Virgo cluster, not present in the observed skymaps (Abreu et al. 2022b).Our work confirms that this tension is lifted by magnetic trapping of UHECRs in Virgo, as was already hypothesized in Biteau et al. (2021) 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 ( ≃ 74 Mpc,  ≃ 5.8•10 14  ⊙ , from Urban et al. 2014), in the direction of which the Telescope Array collaboration claims an indication of excess above 5.7 • 10 19 eV (Telescope Array Collaboration et al. 2021).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 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 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 multi-messenger 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 ( ≃ 100 Mpc,  ≃ 7 • 10 14 M ⊙ , from Planck Collaboration et al. 2013), as well as to determine the expected sensitivity of upcoming gammaray and neutrino facilities at higher energies.

ACKNOWLEDGMENTS 1
The authors would like to thank the reviewer for constructive suggestions, which helped to improve the quality of this work.AC and JB gratefully acknowledge funding from ANR via the grant MultI-messenger probe of Cosmic Ray Origins (MICRO), ANR-20-CE92-0052.AC acknowledges co-developers 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, A 194, Hydra A, A 2634 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); Walker et al. (2017) for Perseus, Briel et al. (1992) for Coma, Govoni et al. (2017) for A 194, Vogt & Enßlin (2003); Enßlin & Vogt (2006); Laing et al. (2008); Kuchar & Ensslin (2011) for Hydra A and Schindler & Prieto (1997) for A 2634.These are compared to experimental data from the Plank/ROSAT project (Eckert et al. 2012), the XCOP project (Eckert et al. 2017), and the ACCEPT data base (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, constant magnetic to thermal energy density for magnetic field.The other possible choices are detailed in Section 2.
In Figure 5, 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 oversimplistic.Magnetic-field estimates derive from the matching of cluster simulations to radio data and/or Faradayrotation 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 above-mentioned caveats in the measurement uncertainties and simplifying modeling assumptions.
Figure3.Interaction and escaping lengths as a function of magnetic rigidity at the center of a prototypical galaxy cluster: photo-hadronic interaction times (dashed-dot 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:  500 = 1 Mpc,  = 1 G,   = 10 kpc,  ICM = 1 • 10 −4 cm −3 .

Figure 4 .
Figure 4. Escaping proton spectra from a cluster of  = 10 14  ⊙ as a function of the injection point.Positive (negative)  valued 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.

Figure 5 .
Figure5.Escaping proton spectra from cluster of different magnetic field values at the center, taken at 1 kpc (see 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.

Figure 6 .
Figure6.Transparency as a function of energy for protons, helium and nitrogen nuclei for different cluster magnetic field (see legend), assuming an 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 parametrization.The vertical line shows the ankle energy.

Fig.
Fig. Set 7. Thermal-proton and magnetic-field densities of reference clusters.