Marion Scarafagio

Alexandre Tallaire

Marie-Hélène Chavanne

M Cassir

A Ringuedé

Diana Serrano

Philippe Goldner

Alban Ferrier
email: alban.ferrier@chimieparistech.psl.eu

Improving the luminescent properties of Atomic Layer Deposition Eu:Y2O3 thin films through optimized thermal annealing

Keywords: Atomic layer deposition, sesquioxide, photoluminescence, lanthanide, quantum technologies

Crystalline rare-earth (RE) doped Y2O3 films are an attractive system for a wide range of photonics applications including quantum technologies which aim at harnessing optical or spin transitions with long coherence times to achieve new functionalities like quantum storage or information processing. Here we present Atomic Layer Deposition (ALD) of Eu doped Y2O3 thin films with improved optical properties. A crucial post-treatment step to obtain high quality films is annealing at elevated temperatures (> 900 °C). However the main drawback of this approach is the formation of unwanted parasitic phases due to reaction at the interface with the substrate, especially with silicon. In this paper, we discuss this issue for different kinds of substrates and buffer layers. The use of such modified substrates allows advantageously extending the maximum thermal treatment up to 1150 °C without being limited by interface reactions. We demonstrate that emission of the 5 D0  7 F2 transition for Eu 3+ in Y2O3 film can be as narrow as that of bulk materials when optimized thermal treatments and a thin undoped Y2O3 buffer layer are used. We thus propose a versatile method in order to reduce the impact of the substrate-film reaction on the optical properties.

Introduction

The development of low cost, compact, highly integrated optical devices will greatly benefit a large range of applications including optical sensing, telecommunication, spectroscopy and quantum technologies (QT). In the framework of QTs, rare earth (RE) doped thin films appear as promising systems along other solid-state materials such as NV centres in diamond or quantum dot semiconductors. Indeed promising results have already been obtained on bulk single crystals as well as on nanoparticles. [START_REF] Goldner | Handbook on the Physics and Chemistry of Rare Earths[END_REF][START_REF] Bussières | [END_REF][3][4][5][6] Extremely long spin coherence times (6 h) have for example been measured at cryogenic temperatures with Er:Y2SiO5 bulk crystals [7] as well as all optical spin coherent control in Eu:Y2O3 nanoparticles. Thin film technologies, on the other hand, are less explored although they could potentially open the way for on-chip integration with other devices such as light sources or detectors. [8][9] To achieve that, a key requirement is the growth of thin crystalline films of rare-earth-doped oxides with low levels of impurities and defects. Besides, the development of high-quality films can also benefit a broader community since new applications of Y2O3 have recently emerged such as catalyst for biodiesel production, wastewater treatment, or anti-corrosion. [10--12] Atomic Layer Deposition (ALD) enables the deposition of luminescent films with a precise control of thickness at the atomic level. This method relies on exposing the substrate surface alternately to different vaporized precursors until surface self-limiting regime is reached (i.e. an almost complete monolayer of reactant is fully adsorbed at the surface). Thanks to that sequential approach, ALD allows unique uniform and conformal deposition as well as dopant distribution engineering that prevents concentration quenching even for relatively large doping contents. [13][14][15] One of the inherent drawbacks of this technique though is the low growth rates associated with this sequential layer-by-layer growth mode together with low deposition temperatures required to remain within this self-limiting regime (e.g. 50°C -400°C). [START_REF] Fanciulli | Rare Earth Oxide Thin Films[END_REF][START_REF] Ritala | M Leskelä Atomic layer deposition Handbook of Thin Films[END_REF] Recently, we demonstrated for RE:Y2O3 10 nm-thin films that a high deposition temperature and post treatment annealing above 900°C are key parameters for optimizing the luminescent properties. In fact, by improving the crystalline environment of the ions, narrow inhomogeneous lines were obtained for the 7 F0 ↔ 5 D0 transition of Eu 3+ even for nanoscale films . [START_REF] Scarafagio | [END_REF] Other motivations to work with the cubic structure of Y2O3 are the isotropic refractive index and for Eu 3+ doped Y2O3 the longest optical and spin coherence lifetime reported for any nanomaterials. [5,6] However some limitations exist. On the one hand, deposition temperature cannot be increased easily since it must remain within the ALD temperature windows where the deposition process is surface limited (temperature range of 250-375 °C for Y(tmhd)3 and O3 precursors). On the other hand the annealing post treatment temperature can only be increased until about 950 °C on silicon above which parasitic silicate phases are formed and modify the emission properties of the RE ions. [START_REF] Scarafagio | [END_REF] In this work, we present strategies in order to further improve the optical properties of thin Eu:Y2O3 ALD films and go beyond those limits. First, we investigated the deposition on alternative refractory transparent oxide substrates such as sapphire (Al2O3) and yttria (Y2O3) transparent ceramics in an attempt to further extend the post-treatment temperature. Another approach is based on the use of buffer layers between the silicon and the luminescent layer. The main idea is to space the reactive interface of silicon from the active RE:Y2O3 layer. We tested two different buffer layers: Si3N4 or in situ grown undoped Y2O3. This approach aims at developing devices compatible with silicon-integrated photonics. [START_REF] Scarafagio | [END_REF]

Results

Rare Earth sesquioxides are polymorphic materials with at least 5 different crystallographic structures. [START_REF] Adachi | The Binary Rare Earth Oxides[END_REF] In the case of Y2O3, two main crystal structures are found: a monoclinic (space group C2/m) and a cubic (Ia-3) one. [START_REF] Adachi | The Binary Rare Earth Oxides[END_REF]20] Unlike bulk crystals for which the cubic structure is predominant, some competition occurs for nanoscale materials and a monoclinic contribution may be detected. In a previous study, we indeed observed that both phases are deposited by ALD on silicon. [START_REF] Scarafagio | [END_REF] The selection during the nucleation between the cubic or the monoclinic structure is complex to predict since it is connected to deposition parameters (temperature, pressure) and to the interaction between the surface and the yttria nuclei. Indeed a correlation between the surface energy and the polymorphism have been demonstrated for nanocrystals. [20,21] ALD process is expected to be particularly sensitive to the nature of the surface since the growth is surface self-limited. Modification of the orientation or the nature of the substrate could lead to texture or new phase crystallisation in comparison to silicon. In order to probe the effect of the substrate on the crystalline environment we used XRD and photoluminescence (PL) of europium ion transitions. In particular, the 5 D0  7 F2 emission is used as a structural probe in order to distinguish the cubic from the monoclinic phase. Indeed, a narrow line centred at 611 nm is expected for the cubic one whereas PL is broader and redshifted for the monoclinic structure (centred around 625 nm). [22][23][24] This large variation is observed because the forced electric dipole (ED) transition 5 D0  7 F2 is hypersensitive (J=2).

Those hypersensitive transitions are excellent environment probe since their intensities are much more influenced by the local symmetry around the RE 3+ ion than the intensities of the other ED transitions. [25,26] 2.1. Deposition on sapphire substrates

The main advantages of sapphire in comparison to silicon are the chemical compatibilities for oxide films (e.g insensitive to oxidation), the transparency and the lower refractive index (1.77 for sapphire against 1.89 for Y2O3 at 600 nm) useful for waveguide technology. [27][28] Furthermore the higher melting point and thermal stability of sapphire in comparison to silicon implies a lower reactivity. The main drawback in comparison to silicon is the larger lattice mismatch with Y2O3 around 5 % (2 % on Si) but a lower thermal mismatch is expected since the thermal expansion of Y2O3 and Al2O3 are closed. [29] Among the different Al2O3 orientations, we selected the a-plane (11-20) since it shows cubic symmetry and can be commercially obtained easily.

Figure 1a presents the XRD pattern of 87 nm-thick films doped with 5 % Eu 3+ and grown on a-plane sapphire for different post-annealing temperatures. The position of the 2 main diffraction peaks: (222) and (400) are indicated. For as-grown and low annealing temperatures (< 700 °C) the film is poorly crystallized and preferentially textured along the <100> directions. [START_REF] Scarafagio | [END_REF] When the annealing temperature is increased, crystallinity improves and a change in texture occurs since a significant contribution of the (222) planes appears. Other groups have already observed similar changes in orientation during the post-annealing treatment. [28,31,32] Until 1000°C no additional phases are observed on the XRD pattern. At 1100 °C, the XRD pattern changes significantly and new broad lines become visible, in particular at 29° and 35°. Those new lines are attributed to the monoclinic yttrium aluminate (YAM) Y4Al2O9 (JCPDS card 00-034-0368). [33] This indicates that the deposited yttria layer has reacted by solid-state reaction with the substrate at this temperature and that the ions are located in a different crystalline environment. This solid-state reaction has already been reported for Y2O3 sol gel film deposited on sapphire. [34] Figure 1. Effect of the annealing temperature on the XRD pattern (a) and on the normalized room-temperature PL emission spectra of the 5 D0  7 F2 transitions (b) of 87-nm-thick Y2O3 films doped with 5 % Eu 3+ and deposited on a-plane Al2O3 (11-20) substrates. The XRD data in (a) are scaled to facilitate comparison. The Y4Al2O9 (YAM) peaks present a good agreement with the JCPDS 00-034-0368. The PL spectrum of the sample annealed at 1100°C (in red) has been vertically translated for better clarity. Nanoparticles and a ceramic with the same composition are respectively used as a reference for the XRD pattern and the PL spectra.

The PL spectra presented in Figure 1b confirms this analysis. In a similar fashion to films on Si substrates, a narrowing of the PL peaks is observed when the annealing temperature is increased. The spectra become closer to that measured for a reference transparent ceramic. Furthermore emission arising from the monoclinic phase of yttria (indicated by an arrow in Figure 1b) is reduced with an increase of the annealing temperature. When the temperature reached 1100 °C however, the FWHM increased significantly indicating solid-state reaction with the substrate. A blue shift of the main emission line from 611 to 607 nm suggests the formation of a yttrium-aluminium mixed phase such as Eu:YAM . [35] In this phase, Eu 3+ ions substitute Y 3+ ions located in a different C1 symmetry site. This multiplication of substitutional sites together with the lowering of crystal symmetry can be responsible for the broadening of the line observed in PL. [36] Therefore, although the use of sapphire as a substrate for Eu:Y2O3 thin film growth allows slightly increasing the maximum annealing temperature as compared to a standard Si substrate, a limitation also exists at about 1100 °C due to the formation of parasitic phases.

Deposition on Y2O3 transparent ceramics

Undoped Y2O3 transparent ceramics have been selected since, due to their having the same composition as the layer, the absence of lattice and thermal mismatch, they may allow limiting stress and interface reactions. The main drawback of this substrate is its noncommercial availability and difficult preparation. [37] The substrates were cut into small pieces of typically 5×5 mm² and mirror-polished. Structural analysis by XRD is difficult because of the predominant signal from the substrate that cannot be easily distinguished from that originating from the layer. For this reason, only PL results are presented. In Figure 2 When the annealing temperature is raised above 900 °C, a significant narrowing of the main line at 611 nm is observed. Simultaneously disappearance of the monoclinic contribution occurs.

Interestingly in this case, we were able to raise the annealing temperature up to 1200 °C with obviously no limitation in parasitic phase unlike other substrates considered previously. At the highest temperature, the PL spectrum is very narrow and the different transitions between stark sublevels at 622 and 627 nm are resolved. Therefore the use of such substrates seems to present no real limits in terms of maximum annealing temperatures due to the absence of interfacial reactions. However one has to take into account that diffusion of ions into the bulk ceramic might occur for higher temperatures and longer annealing times. [38] This can possibly be a limiting factor for the foreseen application in QTs for which spatial localization of the emitters in a thin layer at the surface may be required. [39] Although the transparent ceramic proves to be an ideal candidate for obtaining optimized PL properties of ions in a very high crystalline quality environment, the difficulty of preparation and poor availability of the substrates led us to turn to other strategies. In particular, silicon remains an attractive platform with easier scaling-up and post-processing capabilities. The possibility to limit surface reactions through the intercalation of a buffer layer was thus attempted. Effect of the annealing temperature on the room-temperature PL emission spectra of 82-nm-thick Y2O3 films doped with 5 % Eu 3+ and deposited on an undoped Y2O3 transparent ceramic. A 5 % Eu:Y2O3 transparent ceramic is used as a reference.

Deposition on silicon with a 100 nm-thick Si3N4 buffer layer

The first buffer layer studied is a 100 nm-thick low stress Si3N4 deposited on silicon by PECVD (purchased from University Wafer). Deposition of Eu:Y2O3 luminescent films have already been reported on Si3N4 membrane by e-beam evaporation. Furthermore, Si3N4 is a particularly attractive option for capping passivation or buffer layer due to its high barrier properties, low porosity and high chemical resistance. [40,41] Reactivity between silicon nitride and Y2O3 have been largely reviewed since yttria is used as sintering additive of Si3N4 ceramics. [42][43][44] During the high temperature sintering process, an intergranular crystalline secondary phase of yttrium silicate like -Y2Si2O7 was shown to appear at the grain interface of -Si3N4.

[ 42-44] The XRD pattern and PL spectra are presented in Figure 3 for different annealing temperatures. The XRD pattern for the 900°C annealing presents broad and weak peaks indicating a poor crystalline quality. Furthermore, ratio of the (222) peak area on the (400) one presents a value of about 2 whereas the expected ratio is about 3.6 for random orientation. The use of Si3N4 buffer layer gives a polycrystalline film weakly textured and crystallized. From 1000°C, a significant change in the XRD pattern occurred and new lines appeared. Indexation of these new lines suggested the presence of two yttrium silicate phases indicating a solid-state reaction between the silicon substrate, the buffer and the Y2O3 layer.

The PL spectra confirms the XRD analysis. Indeed for a 900°C annealing post treatment the characteristic PL of europium in the cubic and in the monoclinic phase are both present.

When the post treatment-annealing rose up to 1000°C a significant change in the PL spectrum was observed, in a similar way as for a bare silicon wafer, indicating that a reaction occurred between the oxide and the buffer layer or substrate. Therefore, the use of a Si3N4 buffer was not sufficient to prevent the formation of parasitic phases.

Deposition on silicon with an undoped Y2O3 buffer layer

In order to move away the reactive interface between the substrate and the luminescent layer, we interspaced a Y2O3 undoped layer. This 20 nm-thick buffer layer was deposited in situ by ALD before the active luminescent layer of Eu:Y2O3. XRD pattern for different annealing temperatures is presented in Figure 4. As already observed for deposition on silicon a strong texture along the <100> directions is observed. [START_REF] Scarafagio | [END_REF] This texture remains whatever the annealing temperature in contrast to growth on sapphire. Furthermore, we observed that the diffraction peaks are shifted towards higher angles with respect to the reference which suggests the presence of strain as already described in ref [18]. At 1000°C a careful analysis of the XRD pattern revealed small additional peaks due to the formation of yttrium silicate compounds.

Those new silicate phases (such as Y2SiO5 and Y2Si2O7) are due to chemical reaction between the Y2O3 film and the silicon wafer. For 1150°C annealing temperature, the intensity of the XRD peaks of the silicate phase increased significantly and are indicated by a star on Figure 4a.

On the contrary, no additional lines were observed on the PL spectrum (Figure 4b) even for annealing treatment up to 1100°C whereas on a bare silicon wafer (i.e. without the buffer), emission of RE ions would start to be modified at temperatures as low as 1000 °C. In addition, the Stark sub-levels were resolved for the highest annealing temperatures. This result indicates that the silicate parasitic phase is well localized at the interface between the silicon wafer and the undoped Y2O3 buffer layer but does not extend to the top Eu-doped layer in which the optical properties of the ions are preserved. In this case the maximum temperature could be increased to values as high as 1150 °C leading to a significant improvement of the crystalline environment.

Discussion

Figure 5 and Figure 6 summarize the effect of the annealing temperature on the PL intensity and area of the 5 D0  7 F2 transition for deposition on different substrates. To better visualize the effect we normalized the data with those collected for a 900°C annealing. The results are presented up to temperatures at which parasitic phases appear. We can appreciate that the annealing temperature has two effects. For low temperature annealing (less than 900°C) a significant increase of the integrated emission (PL area of Fig 5b) is observed by a factor of 5 to 8. We attribute this larger emission to an increase of the number of Eu 3+ ions that participate to the PL spectrum. Indeed in this temperature window, the post treatment annealing is expected to remove quenching centres that would otherwise limit total luminescence emission. This hypothesis has been confirmed for deposition on silicon for which a substantial change in the decay time is observed (not presented here). Indeed, for as grown or low-temperature annealed films, the decay time presents a short lifetime with a non-exponential behaviour indicating a strong quenching of the europium emission or energy transfer, whereas at higher annealing temperatures the decay time approaches that measured for bulk crystals. [START_REF] Scarafagio | [END_REF] It is remarkable though on Figure 5b that the PL area remains rather constant when the annealing temperature is further increased. In contrast, the maximum PL intensity continues to rise (Fig 5a). This suggests that crystalline quality improves and that, although the total number of ions contributing to luminescence remains constant, their emission falls within a smaller inhomogeneously broadened line. This is particularly striking for the highest annealing temperatures that show very high PL intensity. The results are finally summarized in Figure 6 which compares the FWHM of the main PL emission at 611 nm for films grown on different substrates. We observe that it strongly decreases with increasing the annealing temperature until the maximum temperature of

Conclusion

In this work we fabricated Eu 3+ doped Y2O3 thin films on different substrates by ALD and compared the structural and photoluminescence properties with standard deposition on Si (100).

For all the substrates, crystalline films were obtained exemplifying the great versatility of the ALD deposition technique. In order to maximize the PL intensity, a high temperature annealing post treatment enhancement is required but the appearance of parasitic interfacial phases is usually a strongly limiting factor. For growth on Al2O3 substrates, the maximum annealing temperature is slight improved (1050°C) as compared to that of bare silicon (950°C).

Furthermore, the use of a buffer layer like Si3N4 did not help in increasing this temperature.

When a transparent ceramic is used as a substrate, no interfacial reactions are formed allowing extending post-treatment up to 1200 °C. We note though that diffusion of the dopants deep into the matrix might occur and could be an issue for QTs applications for which localized dopants are needed. Given the difficulty of preparation and poor availability of the material, transparent Y2O3 ceramics also do not constitute an ideal platform. The use of a thin undoped Y2O3 buffer layer was found to be the most appealing strategy allowing annealing up to 1150°C without modifying Eu 3+ photoluminescence whereas significant changes in the XRD pattern were observed. This suggests that this thin layer acted as a sacrificial interface that reacted with the silicon substrate and prevented the development of the silicate parasitic phase further into the active Eu-doped film. The use of this approach of sacrificial undoped layer is a versatile method that can be applied to other photonics systems.

Experimental Section

Film depositions were carried out with a a Picosun Sunale R200 using conventional βdiketonate precursors: Y(tmhd)3 and Eu(tmhd)3 from (99.9% Strem Chemicals). These stable precursors were held at 140-160 °C and delivered using 300 sccm N2 carrier gas. The vaporized precursors were sent sequentially with an oxidant gas (O3) into the thermalized deposition chamber at 300-350 °C. In order to maximize the film crystallinity, we used ozone as a strong oxidizing agent. A typical ALD sequence consisted of a 3 s Y(tmhd)3 pulse followed by a 3 s purging pulse of N2; a 3 s ozone pulse also followed by a 3 s purge with N2. Eu doping of the oxide films was easily tuned by sequentially introducing pulses of the dopant element in the standard yttrium oxide cycle. In our case, a doping level of 5 % was employed by introducing 5 pulses of Eu(tmhd)3 every 95 pulses of Y(tmhd)3. More details about the impact of the different deposition parameters and their optimization are discussed in details in Scarafagio et al. [START_REF] Scarafagio | [END_REF] The following different substrates were considered for deposition: silicon (100), a-plane Al2O3 sapphire (11-20), Y2O3 transparent ceramic and low thermal stress 100 nm-thick Si3N4 deposited on Si(100). All substrates were used as received after a cleaning with ethanol, acetone with ultrasonic bath. They were all purchased from private companies except for the Y2O3 ceramics that were provided by Dr Ikesue (World Lab Co, Nagoya, Japan). In order to remove experimental bias, deposition on a reference silicon (100) wafer was done for systematic comparison for each experiment. In order to improve the crystallinity of the films and the optical properties of the emitters, we investigated a post-growth annealing treatment at different temperatures under air for 2 hours. The crystallinity of the films was evaluated using X-ray Diffraction (XRD) with a Panalytical XPert Pro diffractometer using an incident beam Ge monochromator, and a stainless-steel sample holder. For all XRD patterns, baseline from Fullprof free software was subtracted.

Film thickness was measured by white light interferometry in the range 250-1000 nm with an Ocean Optics NanoCalc system using Y2O3 optical constants. Finally, photoluminescence (PL) measurements were performed in a Renishaw InVia micro-PL apparatus with a 50 × objective and a 532 nm laser as the excitation source.

  the 5 D0  7 F2 PL for different annealing temperatures are shown. For 900°C annealing and below an important contribution from the monoclinic phase (indicated by an arrow) is visible showing again that crystallinity remains relatively poor when low annealing temperatures are used.

Figure 2 .

 2 Figure2. Effect of the annealing temperature on the room-temperature PL emission spectra of 82-nm-thick Y2O3 films doped with 5 % Eu 3+ and deposited on an undoped Y2O3 transparent ceramic. A 5 % Eu:Y2O3 transparent ceramic is used as a reference.

Figure 3 .

 3 Figure 3. Effect of the annealing temperature on the XRD pattern a) and the room-temperature PL emission spectra (b) of 100-nm-thick Y2O3 films doped with 5% Eu 3+ and deposited on a Si (100) substrate with a 100 nm buffer layer of Si3N4. Nanoparticles and a ceramic with the same composition are respectively used as a reference for the XRD pattern and the PL spectra. The data are scaled to facilitate comparison. Indexation of the diffraction pattern have been done using the JCPDS 04-015-5862 for Y4Si3(SiO4)O10 and 01-086-1105 for Y2Si3O3N4.

Figure 4 .

 4 Figure 4. Effect of the annealing temperature on the XRD pattern (a) and room-temperature PL emission spectra (b) of Y2O3 films doped with 5 % Eu 3+ and deposited on a Si (100) substrate with a 20 nm-thick buffer layer of undoped Y2O3. Nanoparticles and a ceramic with the same composition are respectively used as a reference for the XRD pattern and the PL spectra. Stars indicate the new XRD peaks due to the formation of silicate phase.

Figure 5 .

 5 Figure 5. Comparison of the 5 D0  7 F2 PL for films deposited on silicon, on Al2O3 (11-20) and silicon with a 20 nm thick Y2O3 buffer layer. (a) Comparison of the variation of the normalized integrated PL signal in the range 580-640 nm as a function of the annealing temperature. (b)Variation of the normalized 611 nm peak intensity for different annealing temperature. For each substrate the value of the PL area and PL peak intensity were normalized to that obtained for a 900°C annealing.

Figure 6 .

 6 Figure 6. Comparison of the FWHM of the 5 D0  7 F2 red emission of Eu:Y2O3 as a function of the annealing temperature. A drastic decrease in the FWHM is observed when the annealing temperature is increased. The inset represents a zoom into the high temperature area (850°C-1250°C). The reference value is measured on a transparent ceramic of the same composition.

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 712721 (NanOQTech) and from Région Ile de France through the SIRTEQ framework. Dr Ikesue is thanked for providing transparent ceramic samples.

Received: ((will be filled in by the editorial staff))

Revised: ((will be filled in by the editorial staff)) Published online: ((will be filled in by the editorial staff))


Christopher N Shingledecker
email: cshingledecker@benedictine.edu

Sébastien Incerti

Alexei Ivlev

Dimitris Emfietzoglou

Ioanna Kyriakou

Anton Vasyunin

Paola Caselli

Shi

Cosmic-Ray Tracks in Astrophysical Ices: Modeling with the Geant4-DNA Monte Carlo Toolkit

Keywords:

come   L'archive ouverte pluridisciplinaire

Within the last decade, a number of observational studies have revealed that cold, prestellar cores are far more chemically complex than has been previously assumed. For example, first Öberg et al. (2010), andlater, e.g., Bacmann et al. (2012), Cernicharo et al. (2012), andJiménez-Serra et al. (2016) detected a number of species, referred to as Complex Organic Molecules (COMs), in cold cores, including acetaldehyde (CH 3 CHO), dimethyl ether (CH 3 OCH 3 ), methyl formate (CH 3 OCHO), and ketene (CH 2 CO). More recently, Scibelli & Shirley (2020) found, in a survey of 31 starless and prestellar cores in the Taurus Molecular Cloud (TMC), that 70% contained observable gas-phase abundances of acetaldehyde, and moreover, that methanol (CH 3 OH) was observable towards all sources in their sample. The startling chemical complexity of cold cores was further emphasized by the detection of the aromatic molecule, benzonitrile (C 6 H 5 CN), in TMC-1 by McGuire et al. (2018).

These observational findings are remarkable, in part, because they challenge conventional notions about how such COMs form. It has typically been assumed that COM production occurs mainly within a brief window of time during core collapse in which warming temperatures facilitate the diffusion of radicals on the surfaces of dust-grain ice mantles, as well as the subsequent desorption of COMs, thus produced, into the surrounding gas (Garrod & Herbst 2006;Herbst & van Dishoeck 2009). However, the observations of COMs at earlier, colder stages of star-formation show that the ability of such species to form at low temperatures has been significantly underestimated, with astrochemical models being only partially successful in shedding light on the underlying formation mechanisms in these regions. For example, the gas-grain code of Vasyunin & Herbst (2013) and [START_REF] Vasyunin | [END_REF] was able to qualitatively reproduce the observed abundance of O-bearing COMs by accounting for the increased reactive desorption efficiency on CO-rich ices as well as neutral-neutral reactions efficient at low temperatures. However, this model did not reproduce the observation of N-bearing COMs and overproduced CH 3 OH compared to observations. Cosmic rays provide a likely explanation for the aforementioned conundrum. These energetic particles consist mainly of protons with energies of MeV -GeV (Indriolo & McCall 2013), and are a ubiquitous feature of nearly all astrophysical environments, with the possible exception of protoplanetary-disk midplanes (Cleeves et al. 2013), though recent work by Padovani et al. (2018) suggests that cosmic rays might be important there as well. A large body of experimental work has now shown that the interaction between cosmic rays and ices similar to those coating interstellar dust grains can result in both (a) the production and desorption of COMs such as those observed toward cold cores, as well as (b) drive a variety of interface-dynamical mechanisms that can introduce them into the gas (see, e.g,, reviews by Hudson & Moore (2001), Rothard et al. (2017), andArumainayagam et al. (2019)). Below, we discuss each of these topics in more detail.

Cosmic-ray driven chemistry

During the bombardment of some target material (such as a dust-grain ice mantle) by an energetic primary ion (such as a cosmic ray), the primary ion will collide with the atoms and molecules that comprise the material. Following Bohr (1913), it is customary to divide these collisions into two main categories, namely, inelastic (electronic) and elastic (nuclear) components. For a particle of some energy, E, moving through some material of mass density, ρ, one can thus describe the energy lost per unit path length (dE/dx), as

dE dx = ρ(S e + S n ), (1)
where S e and S n are, respectively, the electronic and nuclear loss functions, with units of cm 2 eV/g. The nuclear component S n is substantially smaller that S e at energies relevant to cosmic rays, and is implicated in changes to the physical structure of the target through, e.g., sputtering and the formation of lattice defect sites, more so than changes to the composition of the target through chemical reactions. Thus, we ignore the elastic component of the primary ion energy loss in this work, and instead focus on the contribution of the inelastic component (Sigmund 1969;[START_REF] Johnson | Energetic Charged-Particle Interactions with Atmospheres and Surfaces[END_REF].

In general, inelastic collisions between the primary ion and atoms in the material result in the ionization and excitation of target species [START_REF] Spinks | An introduction to radiation chemistry[END_REF]Shingledecker & Herbst 2018). Ionizing collisions result in the formation of secondary electrons, which have a broad energy spectrum, but the average energies of which do not exceed around 50-70 eV, depending on the target material [START_REF] Spinks | An introduction to radiation chemistry[END_REF]. Along the trajectories of these secondary electrons, further ionizations and excitations occur. Collectively, the trajectories of the primary ion and all secondary electrons in the target are referred to as the track. Most secondary electrons are stopped near the site of their formation, and thus, the track can approximately be pictured as a cylinder characterized by some radius, r cyl , which we will refer to as the track "core."

Within this cylindrical region surrounding the path of the primary ion, the short-lived excited (suprathermal) species drive a rich variety of reactions at even very low temperatures (< 10 K) that can result in the formation of COMs and even prebiotic molecules (Holtom et al. 2005;Lafosse et al. 2006;Hudson et al. 2008). In Abplanalp et al. (2016) it was shown for the first time that reactions involving these suprathermal species are critical for reproducing the chemistry of cold cores, with later investigations showing that their inclusion in astrochemical models results in significant enhancements of the abundance of COMs such as methyl formate under TMC-1 conditions (Shingledecker et al. 2018).

The energy deposited in the track core also results in a sudden, sharp increase in the temperature of this region (Leger et al. 1985;Bringa & Johnson 2004;Ivlev et al. 2015). This rise in temperature further stimulates chemical changes in the target, and in particular, drives reactions with energy barriers that otherwise could not occur at the equilibrium temperature of the ice mantle. Thus, taken together, the combination of suprathermal reactions, and thermal chemistry in the hot track core represent two promising mechanisms that can help explain the observations of COMs in cold prestellar cores.

Cosmic-ray driven desorption mechanisms

At the interface between the track core and the surrounding vacuum, i.e. the point at which the primary ion enters the target, a hot spot forms in the material with an area of ∼ πr 2 cyl . In this process, known as impulsive spot heating, the increased temperature of the ice surface within this area significantly increases the rate of thermal desorption. As shown by Leger et al. (1985), assuming r cyl < r grain , where r grain is the radius of the grain, this process is independent of the actual grain size.

Conversely, the subsequent process of whole grain heating is not independent of grain size, and occurs as a result of the distribution of the heat deposited in the core throughout the rest of the ice mantle and underlying dust grain. The timescale of this heating has an a 2 dependence on the grain radius, a, such that for interstellar grains with average radii of a ≈ 10 -5 cm, it occurs on the order of nanoseconds (Leger et al. 1985). Fast, exothermic radical-radical recombinations triggered by such 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Energy (MeV) heating could results in the catastrophic loss of the ice mantle through a so-called grain explosion, first noted by Greenberg & Yencha (1973).

A separate mechanism that could likewise trigger grain explosions is impulsive spot heating. This process was studied in detail by Ivlev et al. (2015) who showed that, depending on (a) the value of r cyl , and thus, the volume of the core, as well as (b) the amount of energy deposited therein, a dramatic bifurcation in the fate of ice mantles can occur, characterized by the value of a dimensionless parameter referred to as λ by Ivlev and coworkers. For λ greater than some critical value λ CR = 9.94, the bombardment by an energetic ion will similarly result in the sudden loss of the ice mantle via a grain explosion.

Determining r cyl

From the preceding discussion, it should hopefully be clear that knowledge of r cyl is required for accurate considerations of both the chemistry and interfacial dynamics which occur as a result of the bombardment of a dust-grain ice mantle by a cosmic ray. To date, previous astrochemical works dealing with these topics have typically relied on the estimation from Leger et al. (1985) of r cyl ≈ 5 nm, independent of primary ion velocity (Shen et al. 2004;Ivlev et al. 2015;Kalvāns 2016). Conversely, Bringa & Johnson (2000) proposed that r cyl should increase with dE/dx due to fast (< 1 ps) transport of excitation energy over a few lattice spacings. Based on fits to experimental data, they proposed an expression for r cyl that has a linear dependence on dE/dx and an overall density dependence of ρ 1/3 (Bringa & Johnson 2000).

In principle, however, the radius of the track core qualitatively described by the stopping range of an electron with the average energy of the ejected secondary electrons. Assuming the continuous slowing down approximation, where all particles with the same energy are assumed to travel the same average distance [START_REF] Johnson | Energetic Charged-Particle Interactions with Atmospheres and Surfaces[END_REF] this range for secondary electrons with an average energy W av is given by

r cyl ≈ Wav 0 dE dx -1 electron dE = 1 ρ Wav 0 1 S e,electron
dE.

(2)

Here, (dE/dx) electron is the energy deposited per unit path length for electrons, and similarly, S e,electron is the electronic loss function, also for electrons. In principal, one could use Eq. ( 2) to estimate the track core radius, however such an approach would require accurate analytical expressions for the electronic stopping losses, the derivation of which are beyond the scope of the current work. Nevertheless, Eq. ( 2) is still useful, since it allows us to qualitatively understand the results obtained using the Monte Carlo methods utilized here, for example, by comparing the dependence of our results on material density with the 1/ρ dependence one would expect from Eq. ( 2). Thus, we are presented with three conflicting predictions as to the value of r cyl . In this work, we seek to resolve this confusion and establish more explicitly the value of this critical parameter based on a leading-edge Monte Carlo code, designed to yield astrochemically relevant values for amorphous solid water over a range of cosmic ray proton energies. The rest of this work is arranged as follows: in §2 we provide details of our model and computational approach, in §3 we present the results of calculations and discuss their astrophysical significance, and finally, in §4 we summarize our conclusions.

METHODS

Geant4-DNA

For this work, we have employed the GEANT4 v10.6 Monte Carlo simulation toolkit (Agostinelli et al. 2003;Apostolakis et al. 2009;Allison et al. 2016), which was initially designed to simulate systems relevant for high-energy physics. However, the flexibility of the code allows for its application to problems in a wide variety of fields, including medical physics and astrophysics. The toolkit was later extended by the Geant4-DNA project to simulate microdosimetry through the addition of additional physical processes, such as excitation and charge exchange, that were not included in the original Geant4 code, and allow for the accurate modeling of collisional events down to energies of a few eV (Incerti et al. 2018;Bernal et al. 2015;Incerti et al. 2010a,b). The original motivation for this extension, as indicated by the addition of "DNA" to the name of the toolkit, was to investigate the effects of ionizing radiation in biological systems, including especially DNA and RNA damage.

For this study, we utilized the G4EmDNAPhysics option2 physics list. A full description of the processes and valid energy ranges for particles considered in Geant4-DNA can be found in Incerti et al. (2018), and includes elastic electron scattering, shell ionization cross sections (5 shells), excitation cross sections (5 levels), full secondary electron cascade generation from individual shells (using shell-specific differential ionization cross sections), vibrational excitation, and molecular attachment. In the context of Geant4-DNA, and indeed, of all similar MC codes the track is defined as the collection of the above-mentioned interaction "points," which occur at a given set of x,y, and z coordinates in our simulated volume. Among other restrictions of quantum origin, the spatial extent of the "point" cannot be smaller than the dimensions of the target molecule which, in the case of water, is governed by the ∼ 0.3 nm diameter of the molecule.

Simulation of Ice Bombardment

In the interstellar medium, water ice forms on dust grains through the adsorption and subsequent reaction of, e.g., O, H, OH, and O 2 (Ioppolo et al. 2010;Cuppen et al. 2010). This ice exists mostly as amorphous solid water (ASW). The properties of this glassy metastable material depend on the physical conditions under which it formed, combined with the effects resulting from any changes of these conditions and of any subsequent processing. The two main types of ASW of astrophysical relevance are low-density amorphous (LDA) and high-density amorphous (HDA) ice, which have densities of 0.94 and 1.1 g cm -3 , respectively (Narten et al. 1976;Jenniskens & Blake 1994). In the ISM, the bombardment of LDA by energetic particles similar to cosmic rays has been found to result in its compactification, leading possible to HDA (Palumbo 2005;Mitterdorfer et al. 2014).

By default, the Geant4-DNA simulations include data only for liquid water, however, given the structural similarity of it with the glassy ASW, we here approximate ASW by scaling the density of the material in our model. Since the mean-free-path of electrons is dependent on the density of the material, we here perform calculations at densities relevant to both LDA and HDA. In the context of MC transport of energetic charged particles, several studies have examined the energy-loss function properties of solid water (amorphous and hexagonal ice) and found to be very similar to that of liquid water (Emfietzoglou et al. 2007a(Emfietzoglou et al. ,b, 2008)). Therefore, we expect the energy loss of charged particles in ASW and liquid water to be quite similar and as such, the inclusion of ASW-specific interaction cross sections will not appreciably change the results or conclusions described in the following sections. Perhaps the main uncertainty in the present work comes from the lack of rigorous corrections to the first Born approximation for inelastic electron scattering below about 100 eV. Geant4-DNA has already implemented such corrections in some of its physics models, including the one used in this work, but they are mostly phenomenological. It is possible that these uncertainties may influence very low energy electron transport at the few nm scale, however, an investigation in this matter is beyond the scope of our study.

In our simulations, we represent the ice as a cube with edges 1 µm in length. This ice is then bombarded with protons -the major constituent of cosmic rays -with energies between 100 keV and 100 MeV, which covers both the peak and subsequent fall-off of the electronic stopping power, as depicted in Fig. 1. Incident primary ions are assumed to collide with the ice normal to the surface in the center of the topmost side. Since the stopping length of protons in the energy range considered here is larger than the 1 µm thickness of our ice, they are able to pass completely through, at which point they are considered to have left the system and are not followed any further. For each model run, the simulation begins with the first collision of the primary ion, and ends when all secondary electrons reach energies of around 7.4 eV, below which Geant4-DNA does not currently simulate them.

Finally, in order to aid in the calculation of track core widths, we have disabled elastic scattering processes for the incident protons 1 as well as neutral hydrogen 2 , since the incident H + can capture an electron from the target material. Since, as one can see from Fig. 1, the nuclear component of the stopping power is more than ∼ 2 orders of magnitude less than the electronic component, and moreover, since it does not effect the width of the track core, this assumption should not hinder the accuracy of the resulting calculated values of r cyl .

RESULTS & DISCUSSION

For each initial proton energy in the considered energy range of 0.1 -100 MeV, 1000 simulations were performed, with different randomly-chosen seeds used for each model run. In order to aid in the determination of track core radii, each three-dimensional track was projected onto a (61 nm) 2 2-D surface in the y-z plane of the simulated volume (with the x coordinate here giving the vertical component), represented as a grid of 601 × 601 bins, chosen to encompass the entirety of the track core. For the total track, each bin in our grid was assigned a value equal to the total number of collisions with y-z coordinates within the area covered by the bin, averaged over the 1000 model runs.

One advantage of Geant4-DNA is that the type of each collision simulated in our model is recorded, as well as the resulting energy loss from the primary ion or any other particle generated (e.g. secondary electrons). These data were used to estimate the following four radii:

• Total Ionizations (r ion ): Here, we mean the subset of inelastic collisions which result in the ionization of a bulk species and the concomitant formation of a secondary electron. Electrons formed via collisions of bulk species with the incident primary ion are referred to as firstgeneration secondary electrons, and these, in turn, can be formed with sufficient energy to ionize yet more bulk species, thereby forming second-generation secondary electrons, which, depending on their energies, can generate still higher-generation secondaries. As noted in §2.1, all such electrons are followed in the code until they reach energies of ∼ 7.4 eV.

• Total Excitations (r exc ): Here, we mean the subset of inelastic collisional processes which result in the excitation of a bulk species to some higher bound state. In the code, these collisions can involve either first-generation secondary electrons, or any higher generation of secondary electron. A brief description of the inelastic collisional processes considered in our code is given in §2.1 and the references mentioned there.

• Total Collisions (r tot ): Here, we mean all collisions, both elastic as well as inelastic, e.g. ionization and excitation, by all generations of electrons. This representation of the track is interesting, since it most closely corresponds to the "real" track, as described in §2.1.

• Average Energy Deposition (r energy ): In each inelastic collision, some amount of energy is lost by either the primary ion or secondary electron, and these energy losses are explicitly recorded as outputs of a Geant4-DNA simulation. By thus calculating the average total energy deposited in collisions occurring with y-z coordinates covered by each 1 Å2 bin in our grid, we obtain energy deposition maps from which we can calculate track core radii as described below.

1 proton G4DNAElastic 2 hydrogen G4DNAElastic

These radii were determined by calculating the cumulative distribution of the averaged values for circles of radii from 0.1 -30.0 nm from the origin, i.e. where the primary ion enters the ice, in steps of 1 Å. Based on these cumulative distribution maps the track core was defined as the circular region inside of which 1σ (68.27%) of the events occurred. This value was chosen because it was found to best approximate the fairly stable track core region across the range of energies considered, unlike the much larger 2σ (95.45%) or 3σ (99.73%) radii, which better traced the less numerous secondary electrons of the track "penumbra." Note that in the following discussion we assume that r cyl = r tot . To show the three-dimensional structure of a track, as well as to illustrate a magnified view of our total 1 µm 3 volume, we show in Fig. 2 a portion of the track of a 0.1 MeV proton in LDA ice, along with both the projection onto the y-z plane and the cylindrical, 1σ track-core region. We note that, when averaged over 1000 simulation runs, the irregularities formed by the individual secondary electron paths visible in the projection onto the y-z plane in Fig. 2 are averaged out, and an approximate radial symmetry emerges.

LDA Results

The results of our simulations for LDA ice are listed in Table 1. There, in addition to listing the total radius counting all collisions along the track (r LDA tot ), we further list radii in which only ionization or excitation collisions of the total track were counted. These radii are plotted as a function of initial particle energy in the left panel of Fig. 3. In addition, In Fig. 6 of Appendix A, we show the 1-and 2D representations of the total track, as well as the several sub-tracks considered here. To better illustrate our method of obtaining the values given in Table1, we show in Fig. 4 a representative example of the cumulative distribution plots for a 0.1 MeV H + , as well as contours indicating the track core radii.

From Fig. 3, one can see that the calculated radii for ionization, r ion , are the smallest among those shown. This result is expected since secondary electrons only have the energy to undergo ∼ a few ionization collisions before falling below the ionization threshold of the material. Moreover, ionization cross-sections are generally on the order of ∼1 order of magnitude larger than excitation cross-sections [START_REF] Johnson | Energetic Charged-Particle Interactions with Atmospheres and Surfaces[END_REF], and therefore, the former are expected to occur with a similarly higher probability than the latter (Shingledecker et al. 2017). Thus, for secondary energies above the ionization threshold of the material, ionizing collisions will dominate, and thus, will occur within a shorter distance to the point of formation than for excitation collisions.

Since we track the average energy deposited per primary ion, rather than the total energy deposition, the fact that more energy is lost in ionizing collisions (being equal to the ionization energy of the material) than excitation collisions (< 10 eV) means that r energy , closely follows r ion , though at a slightly larger value due to the contribution of energy deposited during excitation collisions.

One can also see that the total radius of the track core, r tot , is nearly identical to the radius of excitation collisions, r exc , from 0.1 -100 MeV. Over this energy range, these radii display values of ∼ 9.0 ± 1.0 nm. The close association between the total and excitation radii is again not surprising, since, once falling below the ionization threshold, excitation collisions of various kinds, e.g. electronic and vibrational, are the dominant inelastic processes until the electron falls below the excitation energy threshold of the material.

These radii characterize the behaviour of secondary electrons in the material and are thus sensitive the their energies upon formation. As noted by Rudd et al. (1992), for an incident proton with energy E ion , momentum, p, and mass m p , the maximum secondary electron energy, W max which can be produced by that ion is given by

W max (p) = 2m e v 2 0 = 4 m e m p E ion , (3)
where m e and v 0 are, respectively, the electron mass and incident ion velocity. Thus, the maximum secondary electron energy capable of being produced by a 0.1 MeV proton is 218 eV, while for a 100 MeV proton the value is 218 keV. However, the probability of producing such a secondary electron with energy W max is small, with the majority of secondary electrons having a broad energy spectrum but with average energies of not more than 50-70 eV, depending on the target material. The secondary electron energies indicated by our radii do not show the linear dependence on E ion implied by Eq. ( 3). In fact, our results are in agreement with semi-empirical estimations of secondary electron energies by Rudd (1988), who found that the average secondary electron energy, W av , increased until it saturates at around incident proton energies of a few tenths to a few MeV, again depending on the target (cf. Fig. 10 of Rudd (1988)). The results presented in Fig. 10 of Rudd (1988), along with the expression for r cyl given in Eq. ( 2), together provide a very reasonable description of the energy dependence of the radii shown in Fig. 3, where the track core radii will likewise follow changes in W av . This ability to qualitatively describe the findings presented in Fig. 3 in light of Eq. ( 2) and the work of Rudd (1988) is thus a reassuring indicator of the reliability of our approach. We note that the relationship between our derived r cyl values and incident proton energies, though in agreement with the expected average secondary electron energies calculated by Rudd (1988), is not what one would expect based on the formula derived by Bringa & Johnson (2000), which is linearly dependent of the stopping power (dE/dx). As shown in Fig. 1, the electronic stopping power varies by around 3 orders of magnitude over the 0.1 -100 MeV range considered here; however, as shown in Fig. 3, our track core radii vary by only a factor of a few over that range.

HDA Results

Calculated radii for high-density amorphous (HDA) ice are also listed in Table 1, and plotted as a function of energy in the right panel of Fig. 3. The 1-and 2D representations of the total track, as well as the several sub-tracks considered here, are given in Fig. 7 of Appendix B. As with the LDA ice, we see a similar trend where the smallest predicted radii are for ionization, with values that level off at ∼ 5.5 ± 0.5 nm. Likewise, following roughly the same trend, but at slightly larger values, the energy deposition radii have values of ∼ 4.2 nm. Also similar to the LDA ice, the total track core radii closely follow the radii of excitation collisions, again with an average radius of ∼ 8.0 ± 0.5 nm, which is smaller than the LDA value.

Equipped now with results for both LDA and HDA ices, we next turn our attention to a consideration of the dependence of the track radii on the material density. Recall from Eq. ( 2) in §1 that, in principle, the track radius is determined by the stopping range of an electron with the average energy of the ejected secondary electrons, resulting in a density dependence of r cyl ∝ 1/ρ. Using the densities we employed here for LDA and HDA ice of 0.94 and 1.1 g cm -3 , we predict a ratio of r HDA cyl /r LDA cyl = 0.85. Taking the ratio of the r tot values at 6 MeV, we similarly obtain of r HDA tot /r LDA tot = 0.85, thereby nicely recovering the expected density dependence. These results are somewhat at odds with the findings of Bringa & Johnson (2000), whose expression for track core radii has an overall density dependence of r cyl ∝ n 1/3 . Similarly, the peak radii we obtain in our simulations of 9.9 nm and 8.4 nm for LDA and HDA ice, respectively, are roughly double the constant 5 nm value assumed in Leger et al. (1985) and widely used in astrophysical models. However, as with the LDA results, the track core radii we obtain for HDA ice qualitatively follow the behavior we would expect from Eq. (2) and Fig. 10 of Rudd (1988).

CONCLUSIONS

Here, we have reported on the results of calculations we have carried out on the track core radii of energetic protons (characteristic of cosmic rays) in both low-density amorphous (LDA) and highdensity amorphous (HDA) ices. These calculations were performed using the leading-edge microscopic Monte Carlo toolkit, Geant4-DNA, which, despite its initial focus on radiobiological effects, has also proven to be a useful tool in understanding astrophysically-relevant phenomena due to the commonality of the underlying physics. Our main conclusions are the following:

1. track core radii show the weak energy dependence in the range of 0.1 -6.0 MeV expected from previous calculations of secondary electron energies, but that within the range of 6.0 -100.0 MeV the radii values stabilize, 2. the peak track core radii, r cyl , for LDA and HDA ices are, respectively, 9.9 nm and 8.4 nmapproximately double the radii of 5 nm assumed in Leger et al. (1985) and to increase somewhat for incident proton energies below a few MeV, and finally 3. in agreement with the radii predicted from the electron stopping ranges, our results show a density dependence consistent with 1/ρ.

As summarized in §1, an accurate knowledge of track core radii is essential for understanding the chemistry and interfacial dynamics of low-temperature irradiated materials. Moreover, a knowledge of r cyl is essential for predicting the possible importance of grains explosions (Greenberg & Yencha 1973;Ivlev et al. 2015). Thus, these results should be of great importance in further improving how astrochemical models simulate cosmic ray-irradiated interstellar dust-grain ice mantles (Shingledecker & Herbst 2018;Shingledecker et al. 2018Shingledecker et al. , 2020)).

Given the foregoing, there are several key implications of our results, related to the larger values of r cyl obtained here (for incident protons with energies below a few MeV) compared with, e.g., the constant r cyl = 5 nm assumed in Leger et al. (1985) and widely adopted in astrochemical models. First, regarding radiation-chemical changes induced by cosmic ray bombardment, our larger radii imply that the volume in which fast reactions involving suprathermal species occur is larger than previously thought, and that a greater fraction of the ice mantle is involved in such chemistry per collision event. Moreover, as noted in the introduction, desorption at the top of cylinder -i.e. at the ice/vacuum interface -is greatly enhanced due to, for instance, the sharp rise in the temperature immediately following cosmic ray bombardment. Thus, the r cyl values predicted by this work imply somewhat more efficient desorption via this (and related) mechanisms than what might have been estimated in prior studies.

Finally, we note that this initial study proves the utility of the Geant4-DNA Monte Carlo toolkit in understanding processes of interest in astrochemistry. Given the flexibility of the code, future studies using Geant4-DNA could investigate, e.g., how track radii change with different ice compositions or structures.

Figure 1 .

 1 Figure 1. Electronic (S e ), nuclear (S n ), and total (S e + S n ) mass stopping power of protons in liquid water, calculated using PSTAR (https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html)

Figure 2 .

 2 Figure 2. A portion of the track for a 0.1 MeV proton, shown in grey, with the cylindrical track-core region depicted in red and the projection onto the y-z plane shown in light blue. Note: here, the incident ion enters the ice through the center of the top surface.

Figure 3 .Figure 4 .Figure 5 .

 345 Figure 3. Radii calculated with Geant4-DNA for both LDA (left) and HDA (right) ices.

Figure 6 .

 6 Figure 6. 1-and 2D representations of the track for 0.1 MeV protons in LDA ice. The complete figure set (28 images) showing the corresponding plots for all energies considered is available in the online journal.

Figure 7 .

 7 Figure 7. 1-and 2D representations of the track for 0.1 MeV protons in HDA ice. The complete figure set (28 images) showing the corresponding plots for all energies considered is available in the online journal.

Table 1 .

 1 Track radii of energetic protons in amorphous water ice

	Primary Ion Energy	r tot		r exc		r ion		r energy
	(MeV)	(nm)	(nm)	(nm)	(nm)
		LDA HDA LDA HDA LDA HDA LDA HDA
	0.1	3.7	3.1	3.7	3.1	0.8	0.6	1.3	1.0
	0.2	4.3	3.5	4.3	3.5	1.5	1.2	2.0	1.7
	0.3	4.8	4.0	4.8	4.0	2.0	1.6	2.6	2.2
	0.4	5.3	4.5	5.3	4.5	2.4	2.0	3.1	2.6
	0.5	5.8	4.8	5.8	4.8	2.7	2.3	3.5	2.9
	0.6	5.7	4.8	5.7	4.8	2.8	2.4	3.6	3.0
	0.7	6.0	5.1	6.0	5.1	3.1	2.6	3.9	3.3
	0.8	6.3	5.3	6.3	5.3	3.4	2.8	4.2	3.5
	0.9	6.5	5.5	6.5	5.5	3.6	3.0	4.4	3.7
	1.0	6.7	5.7	6.7	5.7	3.8	3.2	4.6	3.9
	2.0	8.2	6.9	8.2	6.9	5.4	4.5	6.1	5.1
	3.0	9.1	7.8	9.0	7.8	6.3	5.6	7.0	6.1
	4.0	9.4	8.0	9.3	7.9	6.8	5.8	7.4	6.3
	5.0	9.7	8.8	9.6	8.7	7.3	6.7	7.8	7.1
	6.0	9.9	8.4	9.9	8.3	7.2	6.2	7.8	6.6
	7.0	9.6	8.1	9.6	8.1	7.0	5.9	7.6	6.4
	8.0	9.4	8.2	9.4	8.2	6.8	6.0	7.4	6.5
	9.0	9.8	8.4	9.7	8.3	7.2	6.4	7.8	6.8
	10.0	9.3	8.0	9.3	8.0	6.8	5.9	7.4	6.4
	20.0	9.3	8.1	9.3	8.0	6.7	5.9	7.3	6.4
	30.0	9.2	7.5	9.2	7.4	6.8	5.3	7.4	5.8
	40.0	8.4	7.5	8.4	7.4	5.6	5.2	6.3	5.8
	50.0	8.3	7.3	8.2	7.2	5.4	4.8	6.1	5.4
	60.0	9.5	7.6	9.4	7.6	7.1	5.3	7.6	5.9
	70.0	8.8	7.6	8.7	7.6	6.4	4.9	6.9	5.6
	80.0	9.0	7.9	8.9	7.9	6.2	5.8	6.9	6.4
	90.0	8.3	7.8	8.2	7.8	5.4	5.9	6.3	6.4
	100.0	8.9	8.4	8.8	8.4	5.9	6.1	6.7	6.7
		4.

ACKNOWLEDGMENTS C.N.S. thanks the Alexander von Humboldt Stiftung/Foundation for their generous support. The work by A.V. is supported by the Russian Ministry of Science and Higher Education via the State Assignment Project FEUZ-2020-0038. A.V. is the head of Partner Group of the Max Planck Institute for Extraterrestrial Physics, Garching, at the Ural Federal University, Ekaterinburg, Russia. I. Kyriakou and D. Emfietzoglou acknowledge financial support from the European Space Agency (ESA Contract No. 4000126645/19/NL/BW)

APPENDIX

A. PROTON TRACKS IN LDA ICE In Fig. 6, we show the 1-and 2D representations of the total track, as well as the several sub-tracks considered here, namely, those showing specifically the excitation or ionization collisions, as well as the average energy deposition. In Fig. 7, we show the 1-and 2D representations of the total track, as well as the several sub-tracks considered here, namely, those showing specifically the excitation or ionization collisions, as well as the average energy deposition.


Arnaud Mangeret
email: arnaud.mangeret@irsn.fr

Jean-Louis Reyss

Marina Seder-Colomina

Lucie Stetten

Guillaume Morin

Antoine Thouvenot

Marc Souhaut

Pieter Van Beek

Marina Seder-Colomina

Early diagenesis of radium 226 and radium 228 in lacustrine sediments influenced by former mining sites

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Uranium extraction activities generate high amounts of potentially hazardous byproducts including waste rocks, mill tailings, and can also impact lands, river and lake sediments [START_REF] Abdelouas | Uranium mill tailings: geochemistry, mineralogy, and environmental impact[END_REF][START_REF] Morin | Mononuclear U (IV) complexes and ningyoite as major uranium species in lake sediments[END_REF], which can result in their enrichment in natural radionuclides such as 226 Ra [START_REF] Carvalho | Radioanalytical assessment of environmental contamination around non-remediated uranium mining legacy site and radium mobility[END_REF][START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF]. Due to the radiotoxicity of this radionuclide, contaminated materials must therefore be appropriately managed in the long term in order to minimize its potential migration to and accumulation in human and ecological targets.

To date, the knowledge on radium mobility and dynamics in low-temperature natural systems is firstly derived from studies conducted in the marine environment. Radium isotopes are indeed frequently used as proxies for the determination of groundwater discharge fluxes into the coastal ocean [START_REF] Moore | Using the radium quartet for evaluating groundwater input and water exchange in salt marshes[END_REF][START_REF] Krest | Using natural distributions of short-lived radium isotopes to quantify groundwater discharge and recharge[END_REF] and transport of pore water solutes (Sun and Torgensen, 2001). In these environments, the distribution of radium is known to be indirectly controlled by redox processes. For example, radium-bearing phases can be sensitive to redox variations, especially iron and manganese (oxyhydr)oxides and oxides, respectively [START_REF] Gao | Mobility of radium and trace metals in sediments of the Winterbeek: Application of sequential extraction and DGT techniques[END_REF][START_REF] Nirdosh | Radium in uranium mill tailings -Some observations on retention and removal[END_REF][START_REF] Nirdosh | Adsorption-desorption studies on the 226 Rahydrated metal oxide systems[END_REF] and sulfate minerals such as barite [START_REF] Doerner | Co-Precipitation of Radium and Barium Sulfates[END_REF][START_REF] Henkel | Diagenetic barium cycling in Black Sea sediments -A case study for anoxic marine environments[END_REF].

The current knowledge of radium geochemistry in natural systems is also derived from studies conducted in continental systems [START_REF] Greeman | Geochemistry of radium in soils of the Eastern United States[END_REF][START_REF] Rihs | Used shortlived nuclides of the U-and Th-series to probe the kinetics of colloid migration in forested soils[END_REF] or in radium-enriched waste such as uranium mill tailings and adjacent surface waters and ponds [START_REF] Sebesta | Behavior of radium and barium in a system including uranium mine waste waters and adjacent surface waters[END_REF][START_REF] Wiles | The radiochemistry of radium and thorium in uranium-mine tailings[END_REF][START_REF] Landa | US Geological Survey-Research on the environmental fate of uranium mining and milling wastes[END_REF][START_REF] Peacey | Field study of geochemistry and solute fluxes in flooded uranium mill tailings[END_REF][START_REF] Martin | The mobility of radium-226 and trace metals in pre-oxidized subaqueous uranium mill tailings[END_REF]. Although a large number of studies have addressed the radium behavior in marine and continental systems, less is known about the processes governing the fate of radium isotopes from reducing environments, except for the study of [START_REF] Martin | The mobility of radium-226 and trace metals in pre-oxidized subaqueous uranium mill tailings[END_REF] conducted in a tailing pond. Notably, the influence of early diagenesis processes on the stability of radium-bearing species in reduced lake sediments is poorly documented. This knowledge is however crucial to choose the appropriate strategy to manage radium-enriched sediments.

The lack of data is notably significant regarding processes that could control radium mobility in anaerobic systems. This gap is partly related to the challenges encountered for i) sampling significant volumes of anaerobic sediment pore waters without modifying radium distribution and ii) measuring low radium isotopes activities in these pore waters. As an example, the study of [START_REF] Gao | Mobility of radium and trace metals in sediments of the Winterbeek: Application of sequential extraction and DGT techniques[END_REF] leads to the quantification of low 226 Ra fluxes from sediments to the overlying waters based on the use of Diffusive Gradient in Thin Films (DGT) techniques and concludes that labile 226 Ra can be related to the reductive dissolution of iron and manganese (oxyhydr)oxides without clearly demonstrating the close association between Ra and Fe association in the solid phase.

The present study aims at providing relevant clues on the geochemical processes controlling the radium distribution in the sedimentary column of a lake. The behavior of 226 Ra is compared with analogous solutes such as 228 Ra and barium. The obtained results reveal important geochemical factors that may be considered for better predicting the post-depositional mobility of radium isotopes in impacted lacustrine systems, as well as, by extension, in underwater radioactive waste storage.

MATERIALS AND METHODS

Site description and sampling procedures for sediment and pore water

The lake Saint-Clément is an artificial reservoir created in 1931 for electricity generation. This lake is located at nearly 20 kilometers downstream from the former uranium mining site of Bois-Noirs-Limouzat and its tailing storage in the Massif Central (Figure 1). A detailed description of the site can be found in [START_REF] Morin | Mononuclear U (IV) complexes and ningyoite as major uranium species in lake sediments[END_REF], [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF] and Stetten et al. (2018a). Lake sediment cores were sampled within PVC liners using a 90 mm diameter Uwitec ® hand corer. Sediments cores were collected at location C2 situated about 100 m upstream from the dam [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF]), in October 2015and May 2016, i.e. cores 2015C2G and 2016C2I, respectively. Immediately after sampling, the cores were placed in a glove bag previously purged three times with N2. Volumes of 15-25 mL of pore waters were extracted from sediments every 5 cm using 0.15 µm pore size Rhizons® samplers. Each of them was connected to a glass vial that had been previously purged with N2 in an anaerobic glove box in the laboratory. Field and laboratory sampling and storage procedures are detailed in Stetten et al. (2018a). Solids were sampled every 5 cm along the core and dried at 30°C before analysis.

Analyses of major and minor elements in sediments and pore waters

The analyses of major and trace elements in the solid sediment samples were performed at the IRSN/LUTECE laboratory with a ICP-OES Thermo Fisher ICap 7600 Duo and an Agilent 8800 ICP-MS/MS (Agilent Technologies, Japan), after acid digestion according to the same procedures as previously described by Stetten et al. (2018a).

The measurements of pH, alkalinity and dissolved Fe(II) in pore waters were performed in an anoxic glove-box (< 5 ppm O2) at the IMPMC laboratory within 48 hours after sampling.

Alkalinity measurements were performed following the procedure described previously by Sarazin et al. (1999). Dissolved Fe(II) was determined by a colorimetric method adapted from the classical orthophenanthroline method (Fadrus and Malỳ, 1975). A detailed description of the two latter analyses can be found in Stetten et al. (2018a).Uranium and trace metals analyses in pore waters were performed at the IRSN/LUTECE laboratory with ICP-MS Thermo Fisher X7 Serie 1. Major and trace elements analyses were performed using an ICP-OES ICap 7600 Duo Thermo Fisher. Measurements were performed using external standard solutions and uncertainties were estimated from 5 replicate measurements on the standard solutions.

Quantification of radionuclides in sediments and pore waters

For sediment samples, aliquots of 1-4 g were used for measuring 238 U and 232 Th-series nuclide activities using the High Purity Germanium well-type detector of the IRSN/LUTECE laboratory (GWL-220-15). 210 Pb, 234 Th (considered in equilibrium with 238 U), 214 Pb and 214 Bi ( 226 Ra), 228 Ac ( 228 Ra), 212 Pb and 208 Tl ( 228 Th) were analyzed after a counting time of one day.

The measurement procedures were similar to those reported by [START_REF] Reyss | Large, low background well-type detectors for measurements of environmental radioactivity[END_REF][START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF]. The gamma detector was calibrated using several international IAEA standards such as RGU-1, RGTh-1, IAEA-312 and 314.

The determination of 226 Ra and 228 Ra activities in pore water samples was performed at the "Laboratoire Souterrain of Modane" in the French Alps [START_REF] Cazala | Improvement in the Determination of 238U, 228-234Th, 226-228Ra, 210Pb, and 7Be by γ Spectrometry on Evaporated Fresh Water Samples[END_REF][START_REF] Reyss | Large, low background well-type detectors for measurements of environmental radioactivity[END_REF].

The pore water samples (15-30 mL) were thereafter evaporated until a final volume of 5 mL was reached. A counting time of 3-5 days was required to achieve a 1statistical error lower than 10% [START_REF] Reyss | Large, low background well-type detectors for measurements of environmental radioactivity[END_REF]. Activities of 232 Th activities in the sediment solid phase were assumed to be at an average level of 140 Bq.kg -1 , as previously reported in the sediments of the two sampling stations located close to the dam of the lake Saint-Clément [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF].

Selective chemical extraction procedures

Selective chemical extraction procedures were applied to the sediment to identify the different phases that carry Ra. The proportions of 226 Ra and 228 Ra in the extractable (FE), reducible (FRED) and oxidizable fractions (FOX) of the sediments were assessed. Two grams of samples representative of sparsely (10-15 cm) and mining-impacted (110-115 cm) sediments were used for each extraction condition. The boundary between the two types of sediments samples was fixed at 92.5 cm depth, due to the simultaneous increase of solid 238 U, 226 Ra and 210 Pb as was previously reported [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF].

First, the sediment was placed with 40 mL of 1M Mg(NO3)2 stirred during 2 h at ambient temperature and then washed with deionized water. This leaching procedure is expected to extract ionically and loosely bound radium, i.e. exchangeable (FE fraction). Second, radium sorbed to or coprecipitated with Fe and Mn, nano-crystalline and crystalline (oxyhydr)oxides (FRED fraction) was extracted using the Citrate-Dithionite-Bicarbonate procedure [START_REF] Mehra | Iron oxide removal from soils and clays by a dithionitecitrate system buffered with sodium bicarbonate[END_REF]. Third, radium associated with organic matter (FOX fraction) was separated by oxidation with H2O2 under nitric media [START_REF] Aguado | Sequential extraction of 226Ra in sediments from an estuary affected historically by anthropogenic inputs of natural radionuclides[END_REF]. Details regarding the protocols for extracting these above-mentioned exchangeable (FE), reducible (FRED), and oxidizable (FOX) fractions can be found in Table 1. After each type of selective extraction, the mixtures were centrifuged at 10,000 rpm and the residue was washed once with the reagent solution. All supernatant solutions were then acidified with 2% HNO3 prior being analysed by gamma spectrometry as described previously.

The Ra activities determined following each sequential leaching were then compared to the Ra activities measured in the bulk fraction.

RESULTS AND DISCUSSION

Vertical evolution of U and Th-series nuclides activities

Solid phase Ra concentration profiles. The activities of the 238 U-serie nuclides in the sediment solid phase are reported in Figure 2 and Table 2 for core 2016C2I and are roughly similar to those measured in a previous study on the same lake [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF]. The activities of 238 U, 226 Ra and 210 Pb are fairly constant until approximately 90 cm depth, reflecting the deposition of the same sediment type until this depth. From 90 to 120 cm depth, the 238 U, 226 Ra and 210 Pb activities increase sharply up to 3862, 1154 and 1245 Bq.kg -1 , respectively. These increases in 238 U-series nuclide activities at depth have been previously interpreted as the legacy of former uranium mining activities that ceased at the beginning of the 1980's [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF]Stetten et al., 2018a). In addition, activities of the 232 Th-series have been formerly characterized in the sediments of the lake Saint-Clément, with an average value of 140 Bq.kg -1 for 232 Th activities in the solid phase [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF]. Figure 1 also displays the profile of 228 Ra activity that is characterized by a continuous decrease with depth, from 260 to 140 Bq.kg -1 . Due to little variation of 232 Th activities, it gave an excess of 228 Ra, which was previously identified in the sediments of the lake Saint-Clément. This excess 228 Ra decrease with depth and this trend was previously interpreted as a radioactive decay and used for dating these sediment deposits [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF]. According to this method, the sedimentation rate for core 2016C2I analysed in the present study is estimated at 5.2 cm.y -1 ; which is exactly the same value obtained previously by [START_REF] Reyss | Estimation of sedimentation rates based on the excess of radium 228 in granitic reservoir sediments[END_REF] on another core that was taken in the close vicinity of our sampling station in the lake. activities of 226 Ra and 228 Ra in pore waters increase with depth, respectively ranging from 0.3 to 9.8.10 -8 and from 0.1 to 0.7.10 -10 µM (Figure 3, Table 3). These activities are higher than in the overlying waters of the lake, where 226 Ra and 228 Ra activities range from 1.4 to 2.3.10 -9 and 1.4 to 2.1.10 -12 µM respectively at 1 and 10 meters depth below the lake surface. Deeper in the sediment, i.e. below 40 cm, the behavior of 226 Ra in sediment pore waters differs from that of 228 Ra (Figure 3). The decrease of pore water 228 Ra activities with increasing depth below 40cm, can be related to its progressive radioactive decay, with its shorter half-life (5,75 years) compared to 226 Ra (1600 years). In contrast, dissolved 226 Ra activity roughly increases to 5-6.10 -8 µM and remains constant from 40 to 90 cm depth and then increases again below 90 cm, up to 9.10 -8 µM (Figure 3).

Radium distribution between mineral and organic phases in the sediment

In order to identify the mechanisms controlling the release of 226 Ra and 228 Ra from sediment to pore water, selective chemical extraction experiments were performed to assess the solid-state speciation of radium in the sediment core. For this purpose, a recent sediment sample and a mining-impacted one were chosen at 30 and 115 cm depth in core 2016C2I; these samples exhibited contrasting 226 Ra activities of 850 and 2100 Bq.kg -1 respectively.

Selective chemical extractions (Table 1) performed on these representative samples indicate that for the surface sediments, the exchangeable fraction contains 6.1 and 10% of total 226 Ra and 228 Ra respectively while the reducible and oxidizable fractions contain 0.7-1.3 % and 0.3-3.7% of total Ra, respectively (Figure 4). For the deep sediments, the exchangeable and oxidizable fractions display few measurable differences compared to the surface sediment samples. The reducible fraction was 6.7 and 10.3% of 226 Ra and 228 Ra, which is slightly higher than for the surface sediments. In soils and sediments, the retention of radium in the solid phase is generally reported to be controlled by adsorption onto ferric (oxyhydr)oxides and some clay minerals, as well as by coprecipitation of mineral sulfates such as barite and strontianite [START_REF] Chen | Radium Sorption to Iron (Hydr)oxides, Pyrite, and Montmorillonite: Implications for Mobility[END_REF]. According to our chemical extractions results, the fraction of Ra associated with Fe and Mn (oxyhydr)oxides or sorbed to clay minerals or organic matter does not exceed 7% of total Ra in the sediment. Due to similar chemical properties of the Ra 2+ and Ba 2+ ions, the 226 Ra/Ba ratio can be proposed as a potential proxy for the mineralogical association between these two elements. Unlike radium isotopes, the distribution of Ba content in sediment exhibits a nearly constant trend for the core studied (Table 2, Figure 5). Accordingly, the vertical evolution of the 226 Ra/Ba ratios in the solid phase of the sedimentary column clearly displays two different patterns, suggesting that i) 226 Ra and Ba may be carried by the same phase down to 87.5 cm depth (nearly constant 226 Ra/Ba ratios) and ii) between 87.5 and 120 cm, another 226 Ra-bearing phase that is likely related to former U mining inputs occurs (increasing 226 Ra/Ba ratios, Figure 5).

Influence of early diagenesis on Ra mobility

Previous measurements of dissolved oxygen at the SWI in the lake Saint-Clément provide evidence of the swift development of anoxic conditions below 1-3 millimeters depth (Stetten et al., 2018a). Redox proxies such as dissolved nitrate and sulfide were previously measured in sediment pore waters and were found to be below the detection limits (Stetten et al., 2018a).

Moreover, pyrite crystals were rare and of submicron size in the entire sedimentary column.

These observations suggest that nitrate and sulfate reduction processes are minor in the sedimentary column.

The vertical profiles of dissolved Fe and Mn in pore waters (Figure 3) are characteristic of the reductive dissolution of iron and manganese (oxyhydr)oxides (Froelich et al., 1979). The occurrence of ferrihydrite in the sediments of the lake Saint-Clément was suggested by the measurement of 8-14% of ascorbate extractable iron in a sediment core collected at the same sampling station, in the previous study by Stetten et al. (2018a). In the same study, XANES and EXAFS analyses at the Fe K-edge indicated that Fe(III) was reduced to > 95% as Fe(II) in the solid phase below 40 cm depth, with iron being dominantly hosted in the structure of clay minerals. Minute iron (oxyhydr)oxide grains were however observed by SEM-EDXS in the 59 -62 cm depth sample (Figure 6), in the 2012-core studied by [START_REF] Morin | Mononuclear U (IV) complexes and ningyoite as major uranium species in lake sediments[END_REF]. Iron was found to be entirely in the form of aqueous Fe(II) in the pore waters (Stetten et al., 2018a).

Dissolved barium in sediment pore waters follows the same trend as dissolved Fe and Mn ; with a maximum concentration of 10 µM at 90 cm depth. The close correlation between dissolved Fe, Mn and Ba profiles support a possible occurrence of adsorbed Ba onto Fe and Mn (oxyhydr)oxides (Carbonell et al., 1999). Similarly, 226 Ra and 228 Ra activities seem to be correlated with Fe(II) in sediment pore waters, which suggests that the release of 226 Ra and 228 Ra from the SWI down to 40 cm depth may be related to the reductive dissolution of Fe and Mn (oxyhydr)oxides (Figure 3). The release of radium in the dissolved phase due to microbial Fe(III) reduction was previously identified using laboratory incubation experiments of uranium mill tailings under anaerobic conditions [START_REF] Landa | Release of226Ra from uranium mill tailings by microbial Fe(III) reduction[END_REF]. However, the contribution of Mn-oxides dissolution to the release of radium cannot be excluded due to the high affinity of these minerals for Ra sorption [START_REF] Goonnea | New perspectives on radium behaviour with a subterranean estuary[END_REF].

Because dissolved 226 Ra and 228 Ra are correlated with dissolved Ba in the sediment pore waters in the first 40 cm under the SWI, we infer that this correlation is best explained by the common release of these elements via the reductive dissolution Fe and Mn (oxyhydr)oxides hosts phases.

It is unlikely that the remobilization of radium in the sediments would be mainly controlled by sulfate reduction and subsequent dissolution of sulfate minerals such as strontianite, gypsum or barite. Indeed, SEM-EDXS observations of the sediments [START_REF] Morin | Mononuclear U (IV) complexes and ningyoite as major uranium species in lake sediments[END_REF]Stetten et al. 2017) indicated the absence of gypsum and the presence of rare but massive barite crystals (Figure 6), and dissolved sulfide was found to be below detection limit of colorimetric methods in the study by Stetten et al. (2017).

Be 40 and 90 cm depth, the apparent correlation between dissolved Ba and Mn remains consistent with the reductive dissolution of Mn oxides, which are known to be excellent scavengers for Ba [START_REF] Carbonell | Effect of phosphogypsum and barite amendments on heavy metals and trace elements chemistry in Mississipi river alluvial sediment[END_REF] This correlation, that is observed along the whole sediment column, suggests that the fraction potentially mobile Ba is dominantly hosted by Mn oxides. Besides, 226 Ra appears to be better correlated to dissolved Fe(II), which suggests that the reductive dissolution of Fe(III) (oxyhydr)oxides may be the dominant process involved in 226 Ra release. This observation may indicate that the potentially mobile fraction of 226 Ra is mainly hosted by Fe(III) (oxyhydr)oxides, without excluding a fraction of 226 Ra associated to Mn oxides.

Below 90 cm, the increase of 226 Ra activity in pore waters is not related to Fe but possibly to the deposition of Ra-enriched particles, as demonstrated by the elevation of solid 226 Ra activity (Figure 2). ). On the contrary, the Fe concentration in pore waters remain globally constant below this depth. The two peaks of Fe concentrations at 42.5 and 87.5 cm in the sediment, correlated with Mn, Zn and Cu (Figure 5), does therefore not induce significant changes in the geochemical behavior of Fe.

Altogether, our results suggest that the increase of radium activities in pore waters down to 40 cm depth below the SWI (Figure 3, Table 3) can be explained by the dissolution of Fe and Mn (oxyhydr)oxides. The release of Ra may however be limited because these phases contain relatively low Ra concentrations, as evidenced by selective chemical extractions (< 7% of total Ra) (Figure 4). Similar behaviour was previously shown in sediments affected by the discharge of 226 Ra associated with phosphate fertilizer plants where radium was also mostly associated with refractory phases [START_REF] Aguado | Sequential extraction of 226Ra in sediments from an estuary affected historically by anthropogenic inputs of natural radionuclides[END_REF]). In the following section, we estimate the radium flux that diffuses through the SWI in order to evaluate the extent of Ra release into the water column.

Diffusive radium fluxes determination through the sediment-water interface

The radium distribution in the pore waters was used to calculate the fluxes of radium that diffuses across the SWI. To do so, we used the Fick's First Law: 𝐽 = -𝐷 𝜑 (Eq. 1)

( [START_REF] Martin | The mobility of radium-226 and trace metals in pre-oxidized subaqueous uranium mill tailings[END_REF][START_REF] Gao | Mobility of radium and trace metals in sediments of the Winterbeek: Application of sequential extraction and DGT techniques[END_REF]Handcock et al., 2010). In these latter studies fluxes calculations were made with 5-7 times higher resolution in the Ra pore water profiles than those in this study. However, our approach contributes to increase the knowledge on radium environmental geochemistry since i) to our knowledge, the data regarding Ra diffusive fluxes in fresh sediments impacted by former uranium mining activities are scarce and ii) this calculation is an useful step for the assessment of risks regarding the degree of radium mobilization. The 226 Ra diffusive fluxes calculated here (3.7.10 -25 mol.cm -2 .s -1 ) are significantly lower than the values previously determined in the sediment of the Winterbeek river impacted by discharges from a phosphate plant in Belgium [START_REF] Gao | Mobility of radium and trace metals in sediments of the Winterbeek: Application of sequential extraction and DGT techniques[END_REF] and in an uranium mill tailings pond at Rio Algom's Quirke Mine in Canada [START_REF] Martin | The mobility of radium-226 and trace metals in pre-oxidized subaqueous uranium mill tailings[END_REF]. These two sites constitute two rare field cases where 226 Ra diffusive fluxes were determined in the context of former industrial activities involving uranium extraction. We note that the diffusive fluxes of 226 Ra and 228 Ra reported here for lake Saint-Clément sediments are in the same order of magnitude as those determined from sediments of the Bega river estuary in Australia [START_REF] Hancock | Using Ra isotopes to examine transport processes controlling benthic fluxes into a shallow estuarine lagoon[END_REF].

The differences discussed above in terms of diffusive 226 Ra fluxes cannot be attributed to the amount of 226 Ra in the sediments. Indeed, the Ra activity in sediments of the lake Saint-Clément is in the range of the Ra activity of the Winterbeek river sediments whereas the radium fluxes reported here are three to four orders of magnitude lower (Table 5). Secondly, for Rio Algom's Quirke Mine, the diffusive 226 Ra fluxes for this site were higher by two orders of magnitude than in the lake Saint-Clément, notably for the station Q1, in which dissolved oxygen concentrations in the bottom overlying waters were close to those encountered in the lake Saint-Clément. However, the 226 Ra activities in mill tailings are not reported [START_REF] Martin | The mobility of radium-226 and trace metals in pre-oxidized subaqueous uranium mill tailings[END_REF] but could be possibly one to two orders of magnitude larger than in the sediments of the lake Saint-Clément [START_REF] Abdelouas | Uranium mill tailings: geochemistry, mineralogy, and environmental impact[END_REF]. Finally, the total solid 226 Ra in Bega estuarine sediments (50-105

Bq.kg -1 ) are 5 to 10 times lower than in those of the lake Saint-Clément. It therefore confirms that the sedimentary Ra activity level alone does not explain the differences of 226 Ra fluxes for each site.

Overall, this short comparative analysis of various sites suggests that the combined effects of the sediment or tailing redox state and of the Ra solid speciation, notably the amount of Ra associated to redox-sensitive phases, can exert an important role in the 226 Ra diffusive flux.

Indeed, the study of [START_REF] Martin | The mobility of radium-226 and trace metals in pre-oxidized subaqueous uranium mill tailings[END_REF] This study reported original data regarding the effects of early diagenesis on the mobility of radium isotopes in the sediments of the lake Saint-Clément, located downstream from a former uranium mining site. Determination of the radium activities was successfully performed for the first time in small pore water volumes, carefully sampled under inert atmosphere, using welltype Germanium detectors, which is challenging.

Our results suggest that the Ra released in sediment pore water of the lake Saint-Clément is associated with the reductive dissolution of Fe and Mn (oxyhydr)oxides. This study demonstrates that even if Fe and Mn (oxyhydr)oxides and possibly natural organic matter are minor Ra-bearing phasesma,, the Ra released in sediment pore waters seems to be significantly driven by reductive dissolution of Fe and Mn (oxyhydr)oxides due to early diagenesis processes. The resulting diffusive Ra fluxes through the SWI remain low in the lake Saint-Clément. These fluxes are significantly lower than those previously reported in Ra-enriched sediments, suggesting that the Ra fluxes in lake Saint-Clément are not driven by the dissolution of Ra-bearing minerals such as e. g. barite (BaSO4).

Further studies are needed for i) constraining the fate of natural organic matter and its association with radium isotopes when reducing conditions develop in the sediment and ii)

improving the characterization of the other mechanisms (e. g. biodiffusion, bioturbation, gas bubbling and resuspension events) that may be involved in the early diagenesis of natural radionuclides..

Fig. 1

 1 Fig. 1 Map of the lake of Saint-Clément (Allier, France) and the core sampling station d from the French

Fig. 2 :

 2 Fig. 2: 238 U and 232 Th-serie nuclide profiles in solid sediment samples (Bq.kg -1 ). 226 Ra/ 228 Ra Activity

Fig. 3 :

 3 Fig. 3: 226 Ra, 228 Ra, Fe, Mn and Ba distributions in pore waters. Solid and empty red and blue circles

Fig. 4 :

 4 Fig. 4: Amount of 226 Ra and 228 Ra extracted from two different sediment samples: exchangeable (FE), oxidizable (FOX) and reducible fraction (FRED). Total 226 Ra and 228 Ra in the pristine samples are

Fig. 5 :

 5 Fig. 5: Fe, Mn, Ba, Cu and Zn concentrations associated in dry sediment samples of the lake Saint-

Fig. 6 :

 6 Fig. 6: Backscattered SEM images and EDX spectra of an iron (oxyhydr)oxide grain in the 59 -62 cm

  indicated that sub-oxic conditions exhibited minor 226 Ra fluxes while the development of anoxia contributed to a larger flux of 226 Ra. Secondly, th increased mobility of 226 Ra was found to be related to the dissolution of redox-sensitive 226 Ra solid phases such as barite.

  Tab. 1: Detailed protocols of selective chemical extraction of radium in sediments

	Reagent	Solid/Liquid Ratio	Experimental conditions	Radium fraction extracted
	1M Mg(NO3)2	1: 20	2h at ambient temperature	Loosely and ionically bound Ra
	0.27M C6H5Na3O7, 2H2O			Ra sorbed and coprecipitated with
	0.1M NaHCO3	2 : 75	40 min at 80°C	amorphous/poorly ordered and well
	2 g Na2S2O4			crystallised iron oxides
	0.02M HNO3 (1) + 30% H2O2; further 30% H2O2 (2)	1: 20	1 => 2 h at 85°C 2=> 3 h at 85°C	Radium associated with organic matter

  Tab. 3: Activities of Ra and concentrations of Fe, Mn and Ba in pore waters for the 2016C2I sediment core Activities of Ra and concentrations of Fe, Mn and Ba in pore waters for the 2015C2G sediment core

	Sample	226 Ra ( 214 Pb + 214 Bi)	228 Ra ( 228 Ac)	Fe(II)	FeT	MnT	BaT
		[× 10 -8 ] µM	[× 10 -10 ] µM		mM		µM
	-10	0.227 ± 0.004	0.021 ± 0.001	-	-	-	-
	7	-	-	0.020 ± 0.002	0.39 ± 0.04	0.014 ± 0.001	1.11 ± 0.01
	17.5	1.064 ± 0.108	0.129 ± 0.026	0.44 ± 0.04	0.56 ± 0.06	0.023 ± 0.002	1.78 ± 0.18
	27.5	1.866 ± 0.144	0.310 ± 0.056	0.93 ± 0.09	1.31 ± 0.14	0.054 ± 0.005	4.22 ± 0.44
	37.5	2.966 ± 0.191	0.572 ± 0.103	1.55 ± 0.15	1.77 ± 0.18	0.070 ± 0.007	5.68 ± 0.58
	47.5	5.023 ± 0.179	0.598 ± 0.064	2.14 ± 0.21	2.20 ± 0.23	0.083 ± 0.008	7.79 ± 0.80
	57.5	4.724 ± 0.239	0.606 ± 0.090	2.34 ± 0.23	2.44 ± 0.25	0.091 ± 0.009	8.52 ± 0.87
	67.5	4.999 ± 0.227	0.520 ± 0.077	2.38 ± 0.23	2.51 ± 0.25	0.095 ± 0.009	8.45 ± 0.87
	77.5	4.843 ± 0.287	0.482 ± 0.103	3.06 ± 0.31	2.60 ± 0.27	0.096 ± 0.010	8.23 ± 0.87
	87.5	5.094 ± 0.239	0.473 ± 0.095	3.38 ± 0.34	2.86 ± 0.30	0.106 ± 0.011	9.69 ± 1.02
	97.5	8.120 ± 0.490	0.439 ± 0.125	2.77 ± 0.28	2.86 ± 0.30	0.102 ± 0.01	9.69 ± 1.02
	107.5	9.459 ± 0.323	0.344 ± 0.082	2.33 ± 0.23	2.70 ± 0.29	0.095 ± 0.009	8.96 ± 0.95
	117.5	9.818 ± 0.419	0.348 ± 0.086	2.85 ± 0.28	2.36 ± 0.25	0.082 ± 0.008	7.86 ± 0.80
	Tab. 4 Sample 226 Ra ( 214 Pb + 214 Bi)	228 Ra ( 228 Ac)	Fe(II)	FeT	MnT	BaT
		[× 10 -8 ] µM	[× 10 -10 ] µM		mM		µM
	7.5	0.287 ± 0.048	0.064 ± 0.017	0.03 ± 0.00	0.13 ± 0.01	0.010 ± 0.001	0.7 ± 0.1
	12.5	0.646 ± 0.048	0.146 ± 0.026	0.24 ± 0.02	0.38 ± 0.04	0.017 ± 0.002	1.0 ± 0.1
	17.5	1.818 ± 0.108	0.314 ± 0.034	0.37 ± 0.04	0.51 ± 0.05	0.021 ± 0.002	1.4 ± 0.1
	22.5	2.284 ± 0.167	0.417 ± 0.064	1.02 ± 0.10	-	0.052 ± 0.005	3.7 ± 0.4
	27.5	2.631 ± 0.144	0,490 ± 0.073	1.51 ± 0.15	1.75 ± 0.18	0.070 ± 0.007	5.2 ± 0.5
	32.5	3.767 ± 0.108	0.563 ± 0.043	1.81 ± 0.18	2.06 ± 0.21	0.080 ± 0.008	6.3 ± 0.6
	37.5	4.903 ± 0.203	0.615 ± 0.082	2.12 ± 0.21	2.24 ± 0.22	0.082 ± 0.008	6.7 ± 0.7
	42.5	5.645 ± 0.263	0.731 ± 0.099	2.25 ± 0.22	2.52 ± 0.25	0.090 ± 0.009	7.4 ± 0.7
	47.5	5.884 ± 0.167	0.666 ± 0.056	2.29 ± 0.23	2.51 ± 0.25	0.095 ± 0.009	8.8 ± 0.9
	52.5	5.286 ± 0.263	0.675 ± 0.107	2.53 ± 0.25	2.60 ± 0.26	0.098 ± 0.009	9.1 ± 0.9
	57.5	5.896 ± 0.227	0.589 ± 0.073	2.43 ± 0.24	2.76 ± 0.28	0.095 ± 0.009	8.4 ± 0.8
	62.5	6.063 ± 0.179	0.516 ± 0.064	2.35 ± 0.23	2.70 ± 0.27	0.098 ± 0.010	8.3 ± 0.8
	67.5	6.458 ± 0.179	0.568 ± 0.064	2.56 ± 0.26	2.85 ± 0.28	0.100 ± 0.010	8.7 ± 0.9
	72.5	6.888 ± 0.957	0.546 ± 0.344	2.72 ± 0.27	2.92 ± 0.29	0.104 ± 0.010	9.3 ± 0.9

  Tab. 5: Concentrations of metals in the solid sediment phase (expressed in mg.kg -1 )

	Sample	Fe (%)	Mn	Cu	Ba	Zn	226 Ra/Ba
	12.5	3.6 ± 0.4	672 ± 68	79 ± 8	588 ± 59	214 ± 7	0.9
	17.5	4.5 ± 0.5	884 ± 89	114 ± 12	709 ± 71	263 ± 8	0.7
	27.5	4.5 ± 0.5	1055 ± 106	168 ± 17	599 ± 60	267 ± 7	-
	42.5	6.0 ± 0.6	1547 ± 155	300 ± 30	582 ± 58	325 ± 8	0.9
	52.5	3.4 ± 0.4	609 ± 61	56 ± 6	641 ± 64	206 ± 7	0.7
	57.5	3.4 ± 0.4	639 ± 64	70 ± 1	624 ± 63	207 ± 7	0.9
	67.5	3.4 ± 0.4	573 ± 58	48 ± 5	636 ± 64	201 ± 7	0.8
	77.5	3.6 ± 0.4	690 ± 69	104 ± 11	656 ± 66	213 ± 7	-
	87.5	6.7 ± 0.7	1844 ± 184	402 ± 4	589 ± 59	366 ± 8	0.8
	92.5	3.7 ± 0.4	710 ± 71	100 ± 10	682 ± 68	234 ± 7	1.1
	97.5	4.0 ± 0.4	800 ± 80	114 ± 12	700 ± 70	243 ± 7	-
	102.5	3.8 ± 0.4	736 ± 74	99 ± 9	692 ± 69	247± 8	1.3
	112.5	3.70 ± 0.37	667 ± 67	93 ± 10	737 ± 74	223 ± 7	1.6
	117.5	3.40 ± 0.34	613 ± 62	83 ± 8	692 ± 69	206 ± 7	-

  Tab. 7: Estimation of diffusive fluxes from the sediments of the lake Saint-Clément (mol.cm -2 .s -1 ) compared with available data in the literature in other Ra-enriched sediments.

		2016		Q1	Q5	W0	W2	W5
		This study	Handcock	Martin et al. (2003)		Gao et al. (2010)
			et al. (2000)
		U-mining impacted	Shallow	Uranium mill tailings	River sediment impacted by a
		lake sediments	lagoon	pond			phosphate plant
	Solid
	sediment
	activity	500	37-42	Not determined	< 180	84	2100
	(Bq.kg -1 )
			Diffusive fluxes (mol.cm -2 .s -1 )
	226 Ra	3.7.10 -25	3.5.10 -25	3.5.10 -23	1.2.10 -20	9.3.10 -22 1.9.10 -21 4.3.10 -22
	228 Ra	1.1.10 -27	1.6.10 -28	-	-	-	-	-

Acknowledgements

We thank Anthony Julien, Olivier Diez and Marine Garcia for their help in field sampling and chemical analyses, respectively. This work was supported by the federator project Needs-

where Jz is the radium diffusive flux ; D i S the diffusion coefficient in the sediment calculated from the data reported previously [START_REF] Yuan-Hui | Diffusion of ions in sea water and in deep-sea sediments[END_REF] ;  the porosity of 0.8 and dc/dz the concentration gradient accross the SWI.

Equation 1 neglects physical processes such as advection, wave-induced transport and gas ebullition in the system. As a first approximation, biodiffusion and bioirrigation processes are also ignored due to the absence of data regarding the abundance and the number of benthic animals in our studied site. A mean porosity value of 0.8 is used for the calculation of the diffusive fluxes. The diffusion coefficient Dw0 in free water was corrected for the effect of temperature according to the Stockes-Einstein relation, assuming an in situ temperature of 15°C

(Table 6). The sediment diffusivity was calculated according to [START_REF] Berner | Early diagenesis[END_REF]. Diffusive fluxes of 226 Ra and 228 Ra were finally determinated according to activities measured at 10 cm above and 7.5 cm below the SWI. The 226 Ra and 228 Ra diffusive fluxes at the SWI in the lake Saint-Clément are reported in Table 7 and compared with fluxes reported in the literature for other types of Ra-enriched sediments