Limits on the Contribution of Endogenic Radiolysis to the Presence of Molecular Oxygen in Comet 67P/Churyumov–Gerasimenko

One of the most surprising findings of the Rosetta mission (Glassmeier et al. 2007) is the high abundance of molecular oxygen (∼3.8% ± 0.85% relative to water) detected by the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) Double Focusing Mass Spectrometer (DFMS) in the coma of comet 67P/Churyumov–Gerasimenko (67P/C-G) (Bieler et al. 2015). The ROSINA observations show that the water and O₂ signals are strongly correlated, suggesting a homogeneous mixing of water ice with inclusions of oxygen molecules rather than the presence of pure O₂ ice (Mousis et al. 2016). The presence of O₂ was independently confirmed by the ALICE UV spectrometer aboard the Rosetta spacecraft (Keeney et al. 2017), but the O₂/H₂O ratios derived from these remote observations are considerably higher (median at 25%, up to 68%). The ALICE observations also suggest the emission of O₂ in bursts (Feldman et al. 2016), implying the presence of O₂ ice. In this investigation, we adopt the ROSINA values because the absorption cross-sections of several species (HS, S₂, CH₄O), which are poorly constrained, may have enhanced the signal attributed to O₂ in the measured UV spectra (Keeney et al. 2017).

Mousis et al. (2016) investigated exogenic sources of radiolysis, such as Galactic Cosmic Rays (GCRs), as a possible origin of cometary O₂. They find that GCR irradiation of the fully formed comet is not a satisfactory explanation, since GCRs would reach only within a few tens of meters from the surface (Cooper et al. 1998), which is the depth of the material lost via ablation at each orbit. They also consider radiolysis of ice grains (before their accretion into the comet) in the protosolar nebula (PSN), but even assuming the highest possible irradiation dose and a maximum yield, this mechanism cannot account for the observed O₂ abundances: only 1% of O₂ can be formed in the lifetime of grains in the PSN, assuming a hundred-fold enhanced GCR flux. On the contrary, they find that radiolysis under the same assumptions in a low-density environment (presolar molecular cloud) can form 1 to 10% O₂, relative to H₂O. In this scenario, O₂ would be included into the cometary material by clathration or condensation in the protosolar nebula (PSN). However, the yields, equivalent to 2.5 molecules per 100 eV, may be too high as suggested by recent laboratory experiments that show a maximum yield of 0.5 molecules per 100 eV (Teolis et al. 2017), and values four orders of magnitude lower for penetrating projectiles; temperatures below 100 K, which can be expected in the molecular cloud, may also lower the yield (Teolis et al. 2009).

Regardless of the production mechanism, clathration of O₂ in the PSN can lead to O₂/H₂O ratios as high as those measured by ROSINA in 67P/C-G but the incorporated O₂ clathrate would form a solid phase distinct from H₂O ice. This is in contrast with the conclusions of Bieler et al. (2015), who observe O₂ to be strongly correlated with H₂O, implying it is homogeneously distributed within the water ice.

In this study, we investigated the possible contribution of endogenic radiolysis of cometary ice to O₂ formation. We consider the comet to be mostly made up of chondritic dust grains covered with an ice mantle. We use the term endogenic to qualify radiolysis of the grains' ice mantle caused by the radionuclides embedded in the chondritic material. These radionuclides potentially represent a spatially homogeneous source of radiation. Short-lived radionuclides (²⁶Al, ⁶⁰Fe) could have significantly heated the nucleus early in the history of the comet (Mousis et al. 2017); the same radiation could have affected its chemical composition. Long-lived radionuclides (⁴⁰K, ²³²Th ²³⁵U, ²³⁸U) would have steadily deposited additional energy over the solar system history. Whereas short-lived nuclides emit their energy through β particles (electrons) and γ-rays, long-lived ones deposit a significant fraction of their energy through α particles (IAEA 2017), more efficient at producing O₂ (Teolis et al. 2017). Considering these two kinds of nuclides, we use the radiolysis model presented in Bouquet et al. (2017) to evaluate whether, within the possible range of yields, there is a credible scenario leading to the formation of O₂ quantities in 67P/C-G comparable to the observations. To make our study relevant to the cometary case, we improved the model by taking into account the possibility that very short range radiations (α particles) can be fully absorbed in the dust grains they are generated in, without ever reaching the surrounding ice mantle.

Section 2 briefly reviews the values of yields given by the experiments described in the literature. Section 3 describes our radiolysis model and its basic assumptions. Section 4 describes the results, followed by Section 5's discussion.

Multiple experiments have been performed to understand the chemistry occurring during, and after, irradiation of water ice. The majority of these experiments feature irradiation by nuclei or electrons (see, e.g., Zheng et al. 2006 and references therein). Commonly identified products are O₂, H₂ and H₂O₂, but their abundances, and even presence, depend on the ice temperature and irradiation type and energy. Experiments point to products quickly reaching "steady states," their abundances being controlled by the equilibrium between their formation and destruction rates (Gomis et al. 2004; Zheng et al. 2006). Therefore, the usual definition of "yield" G_k,i, namely the number of molecules of species k created by the deposition of 100 eV by radiation of type i, should be employed with caution as it is only valid before the equilibrium is reached. Abundances featured in this study are low enough to overlook this issue.

Laboratory experiments have revealed that O₂ yields vary over four orders of magnitude. An explanation for this variability has been recently proposed by Teolis et al. (2017). They found that the yield is dependent on the penetration depth of the projectile and the temperature of the ice target.

Temperatures over 100 K (Teolis et al. 2009) favor the escape of H₂ and therefore the production of O₂ over H₂O₂. Comet 67P/C-G may have reached such temperatures during its history (Mousis et al. 2017), but only for up to 10 Myr, as a result of heating by short-lived radionuclides.

In experiments, O₂ formation occurs mostly in the uppermost 40 angstroms of ice, where H₂ escape is facilitated (Teolis et al. 2017), irrespective of the penetration depth of the radiation. Therefore, the O₂ yield depends on how much of the projectile's energy is deposited in these uppermost 40 angstroms. This has two implications: First, cometary ice may be relatively favorable to O₂ production compared to a laboratory sample of crystalline ice, since it is highly porous and features cavities due to GCR bombardment (Mousis et al. 2016), increasing the yield by offering pathways for hydrogen escape (Grieves & Orlando 2005). Second, projectiles with a short range achieve higher yields. Short range is usually due to the projectile having low energy, and/or high Linear Energy Transfer (LET), as is the case of heavy nuclei. The maximum yield of ${G}_{{{\rm{O}}}_{2}}=0.5$ has been measured for argon nuclei, transferring 100 keV over 100 nanometers (Teolis et al. 2017). The decay of long-lived radionuclides considered in our study generates helium nuclei (α particles) at 4 to 7 MeV¹. This makes ${G}_{{{\rm{O}}}_{2},\alpha }=0.5$ an upper end member for the yield that we used to determine the upper limit of O₂ production.

Electrons (β particles) are another type of radiation emitted by the radionuclides we are considering (typically ⁴⁰K). Irradiation of ice by electrons has been observed experimentally to produce O₂ (Sieger et al. 1998; Orlando & Sieger 2003; Johnson et al. 2005; Petrik et al. 2006; Zheng et al. 2006). However, most experiments that quantified ${G}_{{{\rm{O}}}_{2},e-}$ have been conducted at energies from a few eV to a few keV, while electrons from radioactive decay have higher energies, usually a few tens of keV and up to 560 keV in the case of the β decay of ⁴⁰K.⁷ These experiments show an increase of ${G}_{{{\rm{O}}}_{2},e-}$ with electron energy, from 7 × 10⁻³ for 100 eV electrons impacting ice at 110 K (Sieger et al. 1998), to 0.057–0.228 for 5 keV electrons impacting ice at 12K (Zheng et al. 2006). Higher energy (0.1 to 10 keV) electron irradiation experiments of ∼100K ices found that the sputtered O₂ yield stays at about 1 to 4 O₂ molecules per impactor (Galli et al. 2017). This translates into ${G}_{{{\rm{O}}}_{2},e-}=0.1$ at 1 keV and ${G}_{{{\rm{O}}}_{2},e-}=0.04$ at 10 keV. As an upper end member value for our model, we used ${G}_{{{\rm{O}}}_{2},e-}=0.5$, of the same order as the higher experimental yields. In this work, we cover the intermediary values from ${G}_{{{\rm{O}}}_{2},e-}=0$ to 0.5 O₂/100 eV.

Short-lived radionuclides ²⁶Al and ⁶⁰Fe emit the majority of their decay energy through γ-rays¹. Measurements of ${G}_{{{\rm{O}}}_{2},\gamma }$ are sparser and have been performed at low temperatures favoring H₂O₂ production (Siegel et al. 1961). The notion of LET does not strictly apply to γ-rays. The attenuation of a γ-ray beam is not due to each individual photon progressively losing its energy over its travel distance through many interactions, but rather to the photon traveling mostly unhindered through the material until it interacts with a molecule. As such, the interaction of γ-rays with matter is best described by a mean-free-path, which we find to be 12 cm in olivine and 34 cm in water for a 1 MeV beam (see Appendix A.2). At the considered energies, electron ejection by the Compton effect is the dominant mechanism (Nelson & Reilly 1991), effectively subjecting the ice to electron irradiation. These electrons' energies may be distributed over a wide spectrum (Nelson & Reilly 1991), depending on the scattering angle. In the absence of more experimental data pertaining to γ-rays, we considered the same range of yields as for electrons, i.e., ${G}_{{{\rm{O}}}_{2},\gamma }=0$ to 0.5 O₂/100 eV.

In the case of 67P/C-G, radiation is emitted from inside micron-sized dust grains (McDonnell et al. 1991) and propagates through the ice-dust matrix. If the radiation range in both mediums is much longer than the scale of the grains, a simple relation using the ratios of volume of ice and dust and those of their stopping powers is enough to describe the deposition of energy in ice (see Section 3.1). However, if the range has a scale shorter or equal to those of the grains (case of α particles, see Section 3.2), then a significant portion of radiation emitted by the dust grain may be absorbed without reaching the ice mantle, depending on the core's size (Figure 1). A more complex model, described in Appendix A.1, takes this effect into account.

3.1. Production by $\beta$ Particles and $\gamma$-Rays

If the average range of radiation is longer than the typical grain radius, we can consider the energy is deposited as (Bouquet et al. 2017):

$$\begin{eqnarray}&&{D}_{s,i}=\displaystyle \frac{{\rho }_{r}{A}_{s}{\left(\sum ,E\right)}_{i,s}}{\tfrac{1}{1-\phi }+\tfrac{1}{{S}_{i}^{{\prime} }\phi }}\end{eqnarray} \tag{ 1 }$$

where D_s,i is the rate of energy deposition in ice due to radiation type i emitted by radionuclide type s (energy per unit of volume per unit of time), ρ_r is the grain density, A_s is the activity of species s (decays per unit of time per unit of mass of rock), ${(\sum E)}_{i,s}$ is the energy emitted by species s as radiation type i over one chain of decay (energy, see values in Table 1), ϕ is the ratio of the volume occupied by ice to the total volume (ice + rock), and ${S}_{i}^{{\prime} }$ is the ratio of stopping powers of radiation i between water and rock. The determination of ${S}_{i}^{{\prime} }$ values for an ice/fayalite medium is detailed in Appendix A.2. Our calculations allowed us to derive ${S}_{\alpha }^{{\prime} }=0.52$, ${S}_{\beta }^{{\prime} }=0.26$ and ${S}_{\gamma }^{{\prime} }=0.36$.

Table 1.  Concentration of Radionuclides Used for the Model (4.6 Gyr ago), Half Lives, and Energy Imparted by Decay Sequence for Each Type of Radiation

	Concentration in ordinary	Half-life	α decay energy sum	β decay energy sum	γ decay energy sum
	chondrites (ppb)	(years)	${(\sum E)}_{\alpha ,s}$	${(\sum E)}_{\beta ,s}$	${(\sum E)}_{\gamma ,s}$
			(MeV/decay)	(MeV/decay)	(MeV/decay)
40K	1347.5a	1.25 × 109b	0c	1.1760c	0.1566c
232Th	50.23a	1.4 × 1010b	35.95c	2.8408c	2.2447c
235U	8.06a,d	7.04 × 108b	34.03c	10.4470c	0.55c
238U	24.53a,d	4.46 × 109b	42.97c	6.0935c	1.7034c
26Al	875e	7.27 × 105e	0d	0d	2.6744d
60Fe	2.35d,f	2.6 × 106g	0d	0.148d	2.51d

Notes.

^aVance et al. (2016). ^bIAEA (2017). ^cBlair et al. (2007). ^dLodders (2003). ^eCastillo-Rogez et al. (2009). ^fTang & Dauphas (2015). ^gRugel et al. (2009).

Download table as:  ASCIITypeset image

3.2. Production by $\alpha$ Particles

4.1. Production by $\alpha$-particle Irradiation

We first considered the favorable end member for O₂ production by α particles (energy deposited into the ice given by Equation (1)). We set ${G}_{{{\rm{O}}}_{2},\alpha }$ at 0.5, an upper end member value (see Section 2). We assumed that this value remained constant over the 4.5 Gyr of our simulation. This is also a favorable end member case since it doesn't account for the decline of radiogenic heating by short-lived radionuclides, which allows the comet's temperature to decrease below values favorable to O₂ production. The dust-to-ice mass ratio is considered to be four (Rotundi et al. 2015). We found the final O₂ content of the ice, ignoring contribution from other sources, is only 0.24%, well below the 3.8% ± 0.85% reported by Bieler et al. (2015). Therefore, even under the most favorable assumptions, irradiation by α particles cannot account for the O₂/H₂O ratio observed by ROSINA at comet 67P/C-G.

Zoom In Zoom Out Reset image size

Since irradiation of ice can also produce H₂O₂ with a yield comparable to the maximum O₂ yield (see Section 5), and 0.24% H₂O₂ would be two orders of magnitude larger than observations, we assessed the effect of grain size on the amount of energy available to the ice mantle following the model described in Appendix A.1. We considered a range of grain radii from 1 to 100 μm. The effect of grain size on radiolytic production is shown in Figure 2. While this effect is low for 1 μm grains (it decreases the available energy by 3.4% only), production is reduced by a factor of 10 for a 100 μm grain. In the case of H₂O₂, this would still represent an overabundance, an order of magnitude larger than suggested by observations (Bieler et al. 2015).

Zoom In Zoom Out Reset image size

4.2. Production by $\beta$-particle and $\gamma$-ray Irradiation

Both β particles and γ-rays are decay products of long-lived radionuclides (⁴⁰K, ²³²Th,²³⁵U,²³⁸U) and short-lived radionuclides (²⁶Al and ⁶⁰Fe). The latter represent a high amount of energy imparted to the cometary ice in the first million years of its existence. Due to the absence of a firm value on the yield in O₂ of γ-rays, and the uncertainty on the yield of electrons at energies typical of decay (see Section 2), we ran calculations with values ranging from 0 to 0.5 molecules/100 eV (maximum observed for heavy nuclei) for both radiation types. The energy imparted to ice is deduced from Equation (1).

Figure 3 shows the results of our calculations. The final O₂ percentage in the ice exceeds 1% only in the case of the highest yields considered. We also display an example of time evolution of the O₂ fraction in Figure 4. Short-lived radionuclides produce the majority of radiolytic O₂. The relative contribution of long-lived radionuclides is likely to be even lower if we consider that after the initial heating by short-lived radionuclides, the comet cooled below 100 K, lowering the O₂ yield. Therefore, only very high yields, not supported by experimental data, would allow for a radiolytic production able to match the O₂/H₂O ratio observed at 67P/C-G.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The values of yield required to eventually reach a concentration of 1%, coupled with experimental evidence, lead to the conclusion that endogenic radiolysis is unlikely to account for the quantities of O₂ observed at 67P/C-G. It is therefore more likely that O₂ was already present in the material that formed the comet, possibly after being produced in the interstellar cloud Taquet et al. (2016). Endogenic radiolysis could nonetheless have contributed several tenths of percent, which would have to be considered when investigating possible sources for the steady emission of O₂. This contribution cannot be precisely quantified using the current state of knowledge on ice radiolysis. Additional experimental data on ice irradiation by particles and γ-rays in the range of energies emitted by radionuclide decay are required.

Another product of ice radiolysis is H₂O₂; measured yields range between 0.1 and 0.6 molec/100 eV (Teolis et al. 2017). Observations from ROSINA (Bieler et al. 2015) show an H₂O₂ abundance of at most 0.0023%. From Figure 2, it follows that quantities of H₂O₂ one to two orders of magnitude (depending on dust grain size) higher than suggested by observations would have been produced in 67P/C-G by α particles alone. The H₂O₂ saturation dose (steady-state H₂O₂ concentration reached by an ice sample under constant irradiation) experimentally observed for projectiles of comparable stopping powers is compatible with such a production (Teolis et al. 2017, Figure 8). Additional production by β and γ radiation is harder to quantify. H₂O₂ saturation doses reached under high-energy electron irradiation is lower than with atomic nuclei irradiation (Teolis et al. 2017). This could indicate that β and γ irradiation would have a lower H₂O₂ yield than α particles.

The suppression of H₂O₂ formation due to temperatures above 100 K (Teolis et al. 2017) is unlikely. Mousis et al. (2017) showed that the outer layer of 67P/C-G (several kilometers), to conserve the volatiles observed by Rosetta, would need to have stayed below 100 K for most of the comet's history. Therefore, the part of the comet at the origin of the observations has been under thermal conditions favorable to H₂O₂ production. Another possible factor in suppressing H₂O₂ production is that high-energy radiation may have a lower yield than predicted from stopping power alone. Alternatively, H₂O₂ could have been destroyed by a chemical process yet to be determined. Calculations with the Vienna Ab Initio Simulation package (VASP) (Kresse & Furthmüller 1996) show that H₂O₂ is stabilized within the ice matrix with an energy of −1.4 eV, to be compared with −0.3, −0.5, −1.1, −1.2 eV for O₂, S₂, OH, and HO₂, respectively. This makes it relatively stable and unlikely to diffuse over a long time period, i.e., until it is destroyed by chemical processes. Among them is a multiple process, globally referred to as dismutation (Dulieu et al. 2017):

$$\begin{eqnarray}&&2{{\rm{H}}}_{2}{{\rm{O}}}_{2}\to {{\rm{O}}}_{2}+2{{\rm{H}}}_{2}{\rm{O}}.\end{eqnarray} \tag{ 2 }$$

Ending up with O₂ production, this mechanism has been considered as a possible source for the O₂ observed at comet 67P. However such processes would have to be more than 98% efficient to reduce even a H₂O₂ abundance of 0.1% to the observed value. The outcome of the full chemical network inside the ice critically depends on the activation barriers, which are currently under investigation using first principle solid-state methodology.

On the other hand, experimental yields for the considered range of energies are needed to quantify more accurately how much H₂O₂ would be expected, and narrow down the cause of its near absence in 67P/C-G.

Here we have considered a radionuclide abundance in the comet based on ordinary chondritic values, for all our calculations. The possibility of inhomogeneities of short-lived radionuclides produced in the protosolar disk has been raised (Lee et al. 1998), which could lead to a comet enriched or impoverished in radionuclides. However, the alternate scenario of supernova injection, an efficient mechanism to enrich the PSN in radionuclides, involves extensive mixing of radionuclides into the disk (Ouellette et al. 2009) and would not lead to significant inhomogeneities.

The radiolytic production per unit of volume can be influenced by the size of the body in which it occurs, only if its dimensions are shorter than the typical radiation range. This could have been the case for the building blocks of the comet shortly after CAI, when short-lived radionuclides were most active. Most of their radiation energy would become available for radiolytic production once the boulders become meter-sized. The accretion scenario of Davidsson et al. (2016), even though its timescale is considerably longer compared to the original scenario of Weidenschilling (1997), yields meter-sized boulders within 0.56 Myr after CAI, at most. This delay is insufficient to fully negate the production by short-lived radionuclides (see Figure 4).

Overall, endogenic radiolysis is unlikely to have produced more than a fraction of the O₂ observed in 67P/C-G. Better experimental constraints on the radiolytic production of O₂, as well as its subsequent emission, are required for a more definitive assessment of this fraction. Additionally, we have investigated the potential production of H₂O₂, where our model yields a much higher final abundance (about 0.1%) than that observed by Rosetta (upper limit of 0.0023%). A better understanding of H₂O₂ production yields and of activation barriers in its dismutation is required to understand how the thermal and/or chemical history parameters led to the comet's present-day composition.

Plots were drawn with Matplotlib (Hunter 2007). This work is supported by the Cassini Project through JPL subcontract 1405853 and by the Rosina project through JPL subcontract 1296001. O.M. acknowledges support from CNES. O.M. and T.R. acknowledge support from the A*MIDEX project (n° ANR-11-IDEX-0001-02) funded by the "Investissements d'Avenir" French Government program, managed by the French National Research Agency (ANR).

A.1. Attenuation of $\alpha$ Particles Exiting a Grain

A radionuclide type s inside a dust grain emits, in the form of emission j, a power dP_s,j/dV per unit of volume. "Emission j" designates a given emission in the decay chain (e.g., the 5.013 MeV α particle emission from the ²³¹Pa to ²²⁷Ac transition in the ²³⁵U decay chain). dP_s,j/dV is given by:

$$\begin{eqnarray}&&\displaystyle \frac{{{dP}}_{s,j}}{{dV}}=\displaystyle \frac{{n}_{s}}{{\tau }_{s}}{E}_{0,j}{I}_{j}\end{eqnarray} \tag{ 3 }$$

where n_s is the density of radionuclides in the grain (atoms per unit of volume, deduced from the mass fractions in chondrites shown in Table 1, the grain density and the molar masses of the radionuclides), τ_s is the average decay time of radionuclide s (deduced from the half-life values in Table 1), E_0,j is the energy of emission j, and I_j is the absolute intensity of emission j (i.e., the probability of emission j occurring per decay). n_s is a function of the time t elapsed since CAIs:

$$\begin{eqnarray}&&{n}_{s}(t)={n}_{s,0}\exp \left(-\displaystyle \frac{t}{{\tau }_{s}}\right)\end{eqnarray} \tag{ 4 }$$

where n_s,0 is n_s at the time of CAIs.

A small element of volume dV, situated at a distance r_g from the center of the grain, emits a power dP_s,j isotropically (see Figure 1). We define a z-axis along the radius of the grain passing through dV. Depending on the direction of the emission (angle δ between the trajectory of the radiation and the z-axis), the radiation will travel a distance d_g(δ, r_g) in the grain before exiting. It will undergo an attenuation following the relationship (Dzaugis et al. 2015):

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{j}({d}_{g})}{{E}_{0,j}}={\left(1-\displaystyle \frac{{d}_{g}}{{R}_{\mathrm{stop},j}}\right)}^{\tfrac{1}{{b}_{j}}}\end{eqnarray} \tag{ 5 }$$

where E_j(d_g) is the energy of the particle after traveling d_g, R_stop,j is the stopping distance of particle j in the medium constituting the grain, and b_j is a coefficient determined empirically and dependent on the particle type, energy, and the medium considered.

To obtain values of b_j and R_stop,j, we used the SRIM-TRIM suite⁹ (Ziegler et al. 2010). Through a quantum mechanical treatment of ion-atom collisions, SRIM-TRIM calculates the motion of a particle through a user-defined medium. Each run consisted of 100,000 α particles traveling along the x-axis and hitting a target made of either ice (density of 0.94 g cm⁻³) or fayalite (Fe₂SiO₄), used as a proxy for comet dust grains (Engrand et al. 2016) with a density of 3.0 g cm⁻³. The range of the particle is given by the average length it travels along the direction of emission before stopping. We assumed here a straight propagation of the radiation. Though α particles get deviated by elastic collisions, SRIM-TRIM simulations show that at the high energies considered (4–6 Mev), the lateral deviation is noticeably shorter than the distance traveled along the direction of emission (on average 5% of the penetration depth). We find 4 MeV α particles have average R_stop values of 28.8 and 15 μm in ice and fayalite, respectively, while 6 MeV α particles reach ranges of 52.2 and 26.5 μm in these two media.

The coefficient b_j is similarly determined by fitting Equation (5) to the loss of energy of the simulated particles. We find the initial part of the trajectory is best fit with b_j = 1.2. Since most of the energy is lost in the first half of the particle's travel, we elect to use b_j = 1.2 to adequately represent the rate of energy loss in the grain.

The power is emitted isotropically, which allows us to put dP_s,j, the power due to radiation s, j emitted by dV, under the form

$$\begin{eqnarray}&&{dP}{{\prime} }_{s,j}={\int }_{\delta =0}^{\delta =\pi }\displaystyle \frac{{{dP}}_{s,j}}{{dV}}\displaystyle \frac{{dV}}{\pi }d\delta .\end{eqnarray} \tag{ 6 }$$

Now taking into account attenuation, the power ${dP}{{\prime} }_{s,j}$ coming out of the grain due to radiation s, j emitted by element dV is

$$\begin{eqnarray}&&{dP}{{\prime} }_{s,j}={\int }_{\delta =0}^{\delta =\pi }\displaystyle \frac{{{dP}}_{s,j}}{{dV}}\displaystyle \frac{{dV}}{\pi }{\left(1-\displaystyle \frac{{d}_{g}(\delta ,{r}_{g})}{{R}_{\mathrm{stop},j}}\right)}^{\tfrac{1}{{b}_{j}}}d\delta .\end{eqnarray} \tag{ 7 }$$

The total power due to radiation s, j coming out of the grain, after integration over dV, is given by

$$\begin{eqnarray}&&{P}_{s,j}^{{\prime} }=4\displaystyle \frac{{{dP}}_{s,j}}{{dV}}{\int }_{\delta =0}^{\delta =\pi }{\int }_{{r}_{g}=0}^{{r}_{g}={R}_{g}}{r}_{g}^{2}{\left(1-\displaystyle \frac{{d}_{g}(\delta ,{r}_{g})}{{R}_{\mathrm{stop},j}}\right)}^{\tfrac{1}{{b}_{j}}}d\delta {{dr}}_{g}.\end{eqnarray} \tag{ 8 }$$

To perform the calculation, the explicit dependence of d_g on δ and r_g is required.

To do so, we define a Cartesian referential as shown in Figure 1.

The trigonometric functions of δ give us

$$\begin{eqnarray}&&\cos \delta =\displaystyle \frac{{d}_{g}(0)-({R}_{g}-{z}_{E}(\delta ))}{{d}_{g}(\delta )}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\sin \delta =\displaystyle \frac{{x}_{E}(\delta )}{{d}_{g}(\delta )}\end{eqnarray} \tag{ 10 }$$

it follows that:

$$\begin{eqnarray}&&{z}_{E}(\delta )={r}_{g}+{d}_{g}(\delta )\cos \delta \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{x}_{E}(\delta )={d}_{g}(\delta )\sin \delta .\end{eqnarray} \tag{ 12 }$$

By putting Equations (11) and (12) to the square, adding them together and using ${x}_{E}^{2}+{z}_{E}^{2}={R}_{g}$ we get

$$\begin{eqnarray}&&{d}_{g}^{2}(\delta )+2\cos \delta {r}_{g}{d}_{g}(\delta )+{r}_{g}^{2}-{R}_{g}^{2}=0\end{eqnarray} \tag{ 13 }$$

which leads to

$$\begin{eqnarray}&&{d}_{g}(\delta ,{r}_{g})=\sqrt{{{r}_{g}}^{2}({\cos }^{2}\delta -1)+{{R}_{g}}^{2}}-{r}_{g}\cos \delta \end{eqnarray} \tag{ 14 }$$

allowing one to perform numerically the integration in Equation (8) and obtain ${P}_{s,j}^{{\prime} }$.

A.2. Determination of the Ratios of Stopping Powers

The values of S_i commonly used are ratio of mass stopping powers; it is the case in the formula used as a basis for our calculations (Hoffman 1992). The attenuation of a γ-ray beam going through a material can be described by:

$$\begin{eqnarray}&&I(L)={I}_{0}{e}^{-{\mu }_{l}L}={I}_{0}{e}^{-\mu \rho L}\end{eqnarray} \tag{ 15 }$$

where I₀ is the original intensity, I(L) is the intensity reaching depth L in the material, μ is the mass stopping power, μ_l is the attenuation factor and ρ the density of the material (Nelson & Reilly 1991). S_i values correspond to a ratio of μ values, while our ${S}_{i}^{{\prime} }$ values are a ratio of μ_l values. Equation (15) makes apparent that μ_l is the inverse of the mean-free-path λ. We calculated, for a 1 MeV beam, values of μ in fayalite and ice using the tables of Storm & Israel (1970). The μ values are respectively 0.027 and 0.031 cm²/g. With densities ρ_rock = 3.0 g cm⁻³ and ρ_water = 1.0 g cm⁻³, and using λ = 1/(μρ), we find λ_(fayalite) = 12.4 cm and λ_{(water ice)} = 34.6 cm. Using notations related to our case of interest (water versus a "rock" made of fayalite), we can write these relations:

$$\begin{eqnarray}&&{S}_{i}^{{\prime} }=\displaystyle \frac{{\mu }_{l,\mathrm{water}}}{{\mu }_{l,\mathrm{rock}}}=\displaystyle \frac{{\mu }_{\mathrm{water}}}{{\mu }_{\mathrm{rock}}}\displaystyle \frac{{\rho }_{\mathrm{water}}}{{\rho }_{\mathrm{rock}}}={S}_{i}\displaystyle \frac{{\rho }_{\mathrm{water}}}{{\rho }_{\mathrm{rock}}}=\displaystyle \frac{{\lambda }_{\mathrm{rock}}}{{\lambda }_{\mathrm{water}}}.\end{eqnarray} \tag{ 16 }$$

Values of μ are tabulated for individual elements (Storm & Israel 1970). The mass stopping power of compound materials is given by (Nelson & Reilly 1991):

$$\begin{eqnarray}&&{\mu }_{\mathrm{compound}}=\displaystyle \sum _{j}{\mu }_{j}{w}_{j}\end{eqnarray} \tag{ 17 }$$

where μ_j is the μ value of the element j, and w_j is its mass fraction in the compound.

We calculated S_γ from the tabulated values of Storm & Israel (1970) for water ice and olivine, for 1 MeV γ-rays; we find S_γ = 1.15, and ${S}_{\gamma }^{{\prime} }=0.36$.

Similarly, we calculated ${S}_{\alpha }^{{\prime} }$ and ${S}_{\beta }^{{\prime} }$ through the ratios of their average range in water ice and olivine. We obtained these ranges through SRIM-TRIM and CASINO calculations, for α particles at 4 MeV and electrons at 500 keV, respectively. We find ${S}_{\beta }^{{\prime} }=0.26$ and ${S}_{\alpha }^{{\prime} }=0.52$.