A FLUKA study towards predicting hadron-specific damage due to high-energy hadrons in inorganic crystals for calorimetry

Hadrons emerging from high-energy collisions, as it is the case for protons and pions at the CERN Large Hadron Collider, can produce a damage to inorganic crystals that is specific and cumulative. The mechanism is well understood as due to bulk damage from fragments caused by fission. In this paper, the existing experimental evidence for lead tungstate, LYSO and cerium fluoride is summarised, a study using FLUKA simulations is described and its results are discussed and compared to measurements. The outcome corroborates the confidence in the predictive power of such simulations applied to inorganic scintillators, which are relevant to their adoption as scintillators for calorimetry.


Introduction
The exposure of inorganic calorimeter crystals to large fluences of high-energy hadrons became an issue needing investigation with the construction of the Large Hadron Collider (LHC) and of its detectors at CERN. The choice of lead tungstate (PbWO 4 , or simply PWO) crystals for the CMS electromagnetic calorimeter (ECAL), where non-negligible exposures to high-energy hadrons were expected [1], raised questions on a possible, hadronspecific damage component, while studies of the ionising damage and its minimisation were already well advanced through crystal quality optimisation [2].
A series of pioneering systematic irradiation studies allowed to single out and understand the mechanism of a hadron-specific damage component. Through exposure to 20 -24 GeV protons, it was shown that lead tungstate crystals' light transmission is affected in a cumulative way, with no spontaneous recovery at room temperature [3]. Further, the scintillation mechanism was proven to be left unaffected [4]. Complementary irradiations with charged pions [5] allowed to understand the scaling behaviour between different particle types and the predictive power achievable within the same crystalline material through simulations using the FLUKA package [6]. In parallel, it became possible to draw predictions from data taken in situ for the CMS ECAL detector longevity [7] which led to the decision to replace its endcaps in view of the LHC upgrade for highluminosity running (HL-LHC) [8]. The predictions have recently been confirmed by observations of signal losses in the CMS ECAL that are by now quite significant, especially in the most exposed, high-η regions of the detector [9].
The need to find scintillating materials suitable for calorimetry in view of the HL-LHC upgrade led to study hadron damage in cerium-doped lutetium-yttrium orthosilicate (LYSO) [10] and in cerium fluoride (CeF 3 ) [11]. While LYSO exhibits a cumulative damage, although of smaller amplitude than lead tungstate, cerium fluoride spontaneously recovers from damage at room temperature, in a way which is consistent with a purely ionising damage mechanism [12].
Through all the studies above, the understanding was reached that a hadron-specific damage mechanism is at work in PWO and LYSO, an effect which was also observed in BGO [13,14]. It is due to the extremely high energy loss of fission fragments along their short path in the crystalline matter, leaving small regions of damage behind, that act as scatterers for the scintillation light. The effective reduction of light output is macroscopically observed as a loss in light transmission. Indirect evidence for it was the Rayleigh scattering λ −4 behaviour [15] of the damage amplitudes observed in [3], the visualisation of scattering centres by shining laser light through irradiated PWO and LYSO crystals and the fact that the scattered light is polarised [10,16]. All of these effects were not observed in cerium fluoride [11], as expected, since the material consists of light elements that do not undergo fission [17]. An ultimate proof was reached through the direct visualisation of fission tracks from proton irradiation of lead tungstate, following a method developed for geochronology and used for mineral dating [16].
The damage regions that cause hadron-specific signal losses by scattering the scintillation light have dimensions that are driven by the fission fragments' path lengths. In turn, these are linked to the fragments' stopping power. Fragment production and energy losses are processes that lie in the realm of FLUKA simulations. The possibility to verify how well FLUKA can reproduce measured quantities and thus to validate its predictive power for different materials was therefore found to be of interest for an informed design of future calorimeters and has thus inspired the work presented herein.

Fission track formation
The regions of damage left by fission fragments in minerals are commonly called "fission tracks" in literature. They have been studied foremost as a natural phenomenon occurring in Muscovite mica that is found near an uranium-containing ore and in meteorites that have been exposed to high-energy hadrons from outer space. The understanding of the fission track formation mechanism has been reached through a vast effort documented in literature [18]. The track formation is explained in [19] and references therein as a 3-stage process where a) a charged particle passage causes ionisation, b) ions are ejected from their lattice locations due to Coulomb repulsion and c) a region of atomic disorder is left behind. The main parameter involved is the primary ionisation density, i.e. the number of electrons displaced per unit length. Existing data are collected in figure 1 (left) from [18]. It is evident that track formation is observed above a material-specific ionisation density threshold which is common for all projectile hadrons in a given material, with inorganic materials exhibiting higher thresholds than organic ones. In inorganic crystals, tracks are observed to have been formed by projectiles whose mass number is larger than 30, while no tracks are formed for A < 20. The ionisation rate in figure 1 (left) is given in arbitrary units, however the energy loss per unit distance (dE/dx) threshold value for producing tracks is found by comparison with the plot in figure 1 (right), from [3]. For heavy nuclei, track formation occurs for dE/dx 1 × 10 5 MeV/cm.
Once tracks have been formed, light scatters against them. This phenomenon has been visualised for lead tungstate in [16], where tracks have been revealed through a technique used in geochronology: tracks crossing a surface are revealed through etching, since their material is in a disordered, "metamict" state [19,20], which is more easily soluble than the surrounding crystal lattice. An image of fission tracks for lead tungstate is presented in figure 2, which has been obtained as described in [16], and where typical etched track lengths of a few µm are evident.
One needs to make a distinction, however, between tracks that have been etched for visualisation purposes, and latent tracks, which are those affecting crystals' light transmission in calorimeters. Latent track studies tell us that 2. segment diameters amount to a few nanometers [23] and their length depends on the projectile; 3. light scatters against such segments; 4. the gaps in between segments contain only point defects and are etched more slowly [20].

Rayleigh Scattering
Having observed that fission tracks in scintillating crystals cause scattering of the scintillation light that exhibits all the features of Rayleigh scattering (RS), we summarize herein the relevant quantities of the RS approximation. We shall then attempt obtaining a theoretical input to the equations from FLUKA. The cross section for Rayleigh scattering, σ RS , is given by with d the dimension of the scatterers, λ 0 the vacuum wavelength of the scintillation light and n med the refractive index of the crystal), and finally m = n sc /n med the ratio of the refractive index of the scatterers to that of the crystal [24]. It should be recalled that the RS approximation is applicable for spherical particles that have radii small compared to the wavelength of the scattered light. With the so-called "size parameter" x defined [24] as the requirement is usually formulated as xm 1.
The fraction F RS of light that gets Rayleigh-scattered when going through a crystal of length L is given by where N is the density of scatterers. The intensity I of light transmitted through the crystal is to first order for µ < L, with I 0 the transmitted light intensity before the creation of scatterers through irradiation. By combining equations 3.4 and 3.5, one can express F RS as By combining eqs. 3.1, 3.4 and 3.6 one obtains for the ratio R µ of the induced absorption coefficients between PWO and LYSO: The emission spectrum [25] and the index of refraction as a function of wavelength for LYSO [26] and for PWO [27] have been used to obtain the average scintillation light wavelength in each material. While the scintillation light wavelength in vacuum is nearly identical for LYSO and PWO, in the material it amounts in average to λ LY SO = 248 nm and λ PW O = 198 nm because of the different index of refraction of the two crystals, with uncertainties discussed in section 5.4. Measurements of the ratio between the refraction index of a metamict state and the one of the crystalline lattice, m, for different crystals lie around 0.93 ± 0.01 [28,29]. We thus use m PW O = m LY SO = 0.93 and propagate the uncertainty on each one of them into the determination of the uncertainty on R µ (see section 5.4) that we call ∆R µ . With the values above, one obtains The density of scatterers and their dimension need to be determined in FLUKA, to allow comparing with the measured ratio, R µ = 4.5 ± 0.21 [10].
1uncertainty inferred from the publication

FLUKA simulation setup
For the simulations study, which involves particle interactions and transport in matter, the FLUKA package [6] has been adopted. FLUKA offers the advantage of focusing on the general particle interactions and their consequences, it is especially designed for the evaluation of radiation environments, it allows for averaged material compositions and simple combinatorial geometry and it enables accessing several quantities related to radiation damage. Even non-standard, additional information can be retrieved with user-written routines, as it has been the case for the study reported herein. The results have been obtained with FLUKA version 2011.2x.3 for a total of one million generated primaries reaching each crystal front face. The experimental setups of the irradiations and measurements [3,10,11] have been reproduced in FLUKA. The main parameters of the irradiations were the geometry and composition of the exposed crystals and of the surrounding elements of the irradiation facility IRRAD1 located in the CERN PS T7 beam line [30], the proton beam energy of 20 GeV, respectively 24 GeV, and the beam spot distribution. The crystals studied are: 1. a lead tungstate crystal in the shape of a truncated equal-sided pyramid, 23 [11].
The determination of quantities of interest inside the crystals have been performed according to the experimental measurement conditions. In particular, the light transmission, wherefrom the radiation-induced absorption coefficient is determined, was measured in the experiments with a spectrophotometer whose light beam was aligned on axis with the crystal. Its light beam spot size was 7 mm wide and 10 mm high [3].

Simulations at the infinitesimal level
As a first approach, the main quantities characterising the radiation field have been explored (density of inelastic interactions, dose, density of charged hadrons and of neutrons, displacement-per-atom) as they had been successful in reproducing the scaling of damage in lead tungstate between different particle types and energies [5]. However, when trying to scale between different materials, the ratios of the quantities above, which are of no fundamental nature, are insufficient to reproduce observations. A "microscopic" approach has thus been adopted in the study presented here, where particles are followed at every step of their propagation in matter, and where regions of damage are reconstructed following the experimental observations described in section 2.
For each beam proton impact on the crystal, every secondary particle (briefly "secondaries") has been generated and transported over infinitesimal steps ∆x i . At every step i, the initial particle kinetic energy E kin i has been retrieved and the energy loss has been calculated as dE dx i ∆E i ∆x i with ∆E i the energy lost over the step, and ∆x i has been set by FLUKA to allow the fraction of kinetic energy to be lost in a step to be at most 5%. are reached only in PWO and LYSO, and there just for heavy fragments with mass number A 30, in agreement with literature, as already discussed in section 2. For cerium fluoride, no fragments exceed local energy loss values dE/dx = 1 × 10 5 MeV/cm as would be needed for track formation, in agreement with experimental evidence [11] of no hadron-specific damage in this material. This is also illustrated in figure 4: the fractions of steps by secondaries are plotted, and it is clearly visible where the drop in population occurs for a given A value. For this reason, in the following, studies in cerium fluoride have not been pursued any further, due to the substantiated absence of fission track formation. A first milestone can be considered as reached, however: the agreement between the microscopic observations from simulations described in this section and the experimental results, of no hadron-specific damage in cerium fluoride, constitute a first validation of the predictive power of the FLUKA simulations.

Simulations of fission track segments
Fission tracks can be composed of several segments, separated by gaps where the energy loss falls below the required threshold for track formation (see section 2). In the simulations presented here, consecutive steps have been merged into one track segment if the requirement on energy loss in equation 5.1 is satisfied by both of them. Following this, the length of each segment, the number of segments in a track and the number of segments per incoming proton have been determined. The distribution of secondary particles kinetic energy values versus mass number is shown in figure 5, for a reconstruction requiring an energy loss dE/dx ≥ 1 × 10 5 MeV/cm for individual steps to be joined into segments. The range of mass numbers that give origin to segments in    simulation are observed to lie above A 40, in agreement with literature as summarised in section 2. However, there is an exception: a region of E kin 1 MeV where all mass number values are populated. These entries disappear though if E kin /A > 10 keV is required, as visible in figure 6. The requirement of a minimum threshold in kinetic energy per nucleon finds a justification if one estimates the infinitesimal displacement of such a secondary: the maximum displacement that can occur is ∆x max 10 keV/(1 × 10 5 MeV/cm) = 10 Å, which is of the order of the typical crystal lattice cell size. Such secondaries are therefore fragments that do not have sufficient kinetic energy to move away from their lattice location, and thus cannot possibly produce a track.

From segments to tracks
For the determination of track densities, towards a comparison with observed ones, the full spectrum of simulated track length distributions that are plotted in figure 7 needs to be considered. To allow appreciating more easily the features of the distributions, these have been histogrammed once with a vertical logarithmic scale, and once with a horizontal logarithmic one. One observes that the length of segments from simulations extends over three orders of magnitude, between a few nm and a few µm. A comparison can be made with figure 5 from [31], where two extreme cases are presented: the spectrum of track lengths observed in earth rocks containing radioactive isotopes is limited to short tracks (< 15µm), since these are due to fission recoils, while in meteorites a tail of longer tracks can be observed, due to the exposure to heavy, very high-energy cosmic rays, extending to 160 µm. The distributions in our test (figure 7, left) have a shape compatible with the observations in terrestrial rocks. This is understandably the case, since the high-energy primary protons in the irradiation tests produce a hadron shower in the crystalline matter, where heavy fission fragments just reach intermediate energies, compared to the two extremes quoted above.
The number of segments belonging to the same track are histogrammed in figure 8. While segments in latent tracks are scattering centres for light, in the process of track etching, which is performed for visualisation purposes, those segments get joined, as explained in section 2 and in [20]. This is evident in the images of etched tracks, as found in figure 2 (left) or in [16] e.g., where no pattern of tracks split into individual segments can be observed. To reach a further validation of the simulations, densities have thus to be compared at the track level.
From [16], where fission tracks from lead tungstate have been visualised, one obtains an observed track density crossing a surface per proton entering the crystal that amounts to φ = (1.76 ± 0.14) × 10 −6 tracks cm −2 /(p cm −2 ). In the FLUKA simulations for lead tungstate, the irradiation  has been performed with a proton fluence of 1 × 10 6 protons on the crystal front face, i.e. Φ p = 2.07 × 10 5 p cm −2 . From the number of entries in figure 7 and the lead tungstate crystal volume, the density of segments in the crystal results to be ρ seg = (1.01 ± 0.25) × 10 −2 segments cm −3 /(p cm −2 ). Here, the statistical uncertainty is negligible, the dominating contribution being due to the dE/dx threshold used. From figure 8 one reads off an average of 1.3 segments/track and thus a track density is obtained from FLUKA ρ trk = (7.8 ± 2.0) × 10 −3 tracks cm −3 /(p cm −2 ) for PWO. One can estimate the average track length considering that only tracks with a length of the order of L will be able to traverse a volume of thickness L. To get an estimate, we can thus write approximately wherefrom, in our case, we obtain L (2.3 ± 0.6)µm. The average length obtained can be compared to the average etched track length visible in figure 2, where it becomes evident that the calculated and the observed length agree in order of magnitude, yielding a second, independent validation of the FLUKA simulations.
It might be useful to point out that the track length was inferred from densities, rather than from a direct determination of the simulated track length: this is justified because we want to compare simulation findings with the length of observed etched tracks [16], which tend to be longer than latent ones. It should also be noted that the only tracks directly accessible to simulations are latent ones, and that their length in lead tungstate is found to be on average 2.2 µm, as can be inferred from figures 7 and 8.

Rayleigh Scattering ratio
After having validated the simulations in two independent ways, confidence was gained that the typical dimensions and densities of segments from simulations can be used to obtain the RS ratio between different crystalline materials. Recalling eq. 3.8, we shall call the two intervening fractions and To obtain an estimate of the Rayleigh scattering ratio, one has to take into account some caveats: the RS approximation assumes spheres, while dipole-shaped, randomly oriented tracks are present in calorimeter crystals after exposure to high-energy hadrons. However, for cylindrical scatterers, RS behavior has been observed [32], where either the radius of equivalent-surface spheres or the smallest particle dimension appears to be a good parameter for an estimate [33].
To understand the reliability of simulations, we estimate the effective scatterer dimension for each crystal type individually using the experimental results for the damage amplitudes and the already validated FLUKA segment densities.
From Eqs. 3.4 and 3.6, it follows that µ = N × σ RS . (5.6) This equation, combined with Eq. 3.1, allows to extract the relevant dimension for Rayleigh scattering. For PWO irradiated with a proton fluence Φ p = 2.07 × 10 5 p cm −2 , the density of segments from FLUKA amounts to N = 2100 segments cm −3 and the radiation-induced absorption coefficient is µ = 4.57 × 10 −8 m −1 , as of figure  15 in [3], so that the effective scatterer dimension obtained is d PW O = 32.0 nm (uncertainties are discussed below). For LYSO, analogously, the density of segments from FLUKA amounts to N = 1006 segments cm −3 for an irradiation with a proton fluence Φ p = 1.6 × 10 5 p cm −2 . The measured radiation-induced absorption coefficient is µ = 7.8 × 10 −9 m −1 , as obtained from [10] combined with figure 15 in [3]. These values lead to an effective scatterer dimension d LY SO = 31.7 nm. These sizes, well compatible between the two crystal types, are nearly an order of magnitude smaller than the scintillation emission wavelength, thus compatible with the requirements of the RS approximation in 3.2. Before considering the uncertainties affecting the estimated dimensions, we compare them to the equivalent-sphere radii that result from purely geometric considerations.
The latent segment lengths obtained from FLUKA are, for PWO, L = 1700 nm and, for LYSO, L = 920 nm ( figure 7). The latent track diameter found in literature is approximately 2.5 nm for Zircon [23]. Using this value yields an estimate for the surface-equivalent sphere radius R PW O = 46 nm for PWO and R LY SO = 34 nm for LYSO. While this estimate is strongly affected by the uncertainty on the track radius, which is not known for the two crystals considered, it offers a validation of the order of magnitude of the values determined above.
From the estimates based on experimentally validated quantities above, the assumption that R d 1 appears justified. The RS ratio is then obtained from the segment densities per proton hitting a unit surface derived from the total number of segments (see e.g. figure 7), which are, for PWO, ρ PW O seg = 1.01 × 10 −2 segments cm −3 /(p cm −2 ) and, for LYSO, ρ LY SO seg = 0.63 × 10 −2 segments cm −3 /(p cm −2 ), with uncertainties discussed below, resulting in a RS ratio R µ = 4.0.
Uncertainties on the prediction for R µ arise from various sources: (a) the uncertainties on the refractive indices for PWO (±0.01 [27]) and for LYSO (±0.005 [26]) contribute with ∆R µ,n = ±0.08; (b) uncertainties in the shape of the emission spectra of LYSO and lead tungstate affect the precision of the factor λ LY SO /λ PW O 4 determination and contribute with ∆R µ,em = ±0.2; (c) the uncertainties on the two ratios of refractive indices between metamict and crystalline state for PWO and for LYSO, ∆m PW O = ∆m LY SO = ±0.01, cause a significant contribution ∆R µ,m = ±1.60; (d) The uncertainty on the applicable dE/dx threshold affects the segment densities for PWO and LYSO in a correlated way, resulting in a contribution ∆R µ,dE dx1 = ±0.6; (e) A common dE/dx threshold has been used for both, PWO and LYSO. Allowing for a shift between the two threshold values causes an uncertainty ∆R µ,dE dx2 = ±0.8; (f) The uncertainties on d for PWO and for LYSO from segment densities contribute with ∆R µ,d = ±0.24; (g) The total statistical uncertainty of N PW O and N LY SO contributes with ∆R µ, N = ±0.04; The total uncertainty thus amounts to ∆R µ = ±1.9 and the FLUKA prediction is given by R µ = 4.0 ± 1.9. (5.7) This result has to be compared with the measured ratio R meas µ = 4.5±0.2: the Rayleigh scattering amplitude ratio between PWO and LYSO obtained from FLUKA simulations is compatible, within the uncertainties, with the measured one. However, differences in scatterer dimensions could lead to significantly different R µ values, due to the presence of a 6 th power dependence.
This last validation, of the RS ratio, along with those described in section 5.1 for the fulfilment of trackformation requirements and in section 5.2 for the segment lengths, corroborates the possibility of estimating, through simulations, the order of magnitude of relative damage amplitudes to be expected from hadrons in different inorganic crystals.

Conclusions
A simulation study has been performed of long, inorganic crystals irradiated by 24 GeV protons producing hadron showers in the materials, and several observable quantities have been compared with measured ones. The simulations confirm the qualitative understanding that has been reached of the hadron-specific damage component in the studied inorganic scintillating crystals. In particular, FLUKA simulations yield no heavy, highly ionising fragments in cerium fluoride, as would be needed for track creation, in agreement with the observed absence of a hadron-specific damage component [11]. Instead, the simulations yield heavy, highly ionising fragments in PWO and LYSO, as needed for track creation, in agreement with the observed hadron specific damage [3,10]. In lead tungstate, the simulations yield track densities that are compatible with experimentally observed ones [16]. Further, the simulated spectrum of segment lengths is in qualitative agreement with experimental data [19]. Finally, FLUKA simulations yield a Rayleigh Scattering amplitude ratio between PWO and LYSO that is consistent, within the uncertainties, with the measured one [10]. All these different points of validation confirm the reliability of FLUKA simulations, deployed at a microscopic level according to criteria that were established experimentally [18], in estimating the order of magnitude of damage densities and amplitudes that should be expected from hadrons in different inorganic crystals, enabling thus an informed choice of materials in calorimeter design.