Cécile Gautheron
email: cecile.gautheron@u-psud.fr

Mbongo Duval

Djimbi

Jérôme Roques

Hilal Balout

Richard A Ketcham

Eric Simoni

Raphael Pik

Anne-Magali Seydoux-Guillaume

Laurent Tassan-Got

A multi-method, multi-scale theoretical study of He and Ne diffusion in zircon

The quantification of He and Ne diffusion behavior in crystals rich in U and Th such as zircon is key for the interpretation of (U-Th)/ 4 He and (U-Th)/ 21 Ne thermochronometric ages.

Multiple parameters such as chemical substitution, channel obstruction and damage can modify the diffusivity compared to a pristine structure. To investigate the impact of these parameters, we have conducted a theoretical diffusion study combining a series of methods and approaches to address the problem across the necessary range of scales (atomic to crystal size). First, using quantum calculation, we determine the different He and Ne insertion sites, insertion energies and diffusion pathways at the atomic scale for an ideal pristine zircon structure (i.e. damage free). These results serve as input for a 3D random walk simulation of atomic trajectories that provides diffusion coefficients for damage-free zircon crystals. Second, as natural zircon crystals are not perfect, we model the impact of different types of damage and diffusion pathway obstruction at the atomic level on He and Ne diffusion in 3D. The calculated He and Ne diffusion coefficients for pure ZrSiO4 exhibit strongly anisotropic behavior and very high diffusivity along the c-axis, and with 3D, closure temperatures of -197°C and -202 °C respectively. The results for He are comparable to previous DFT studies but strongly different from experimental diffusion results; results for Ne are similar in this respect. Modelling the impact of different types of damage (vacancies, recoil, fission, voids or fluid inclusions) and obstruction on He and Ne diffusion reveals important implications for the (U-Th)/He and (U-Th)/Ne thermochronometers. First, obstruction alone does not significantly modify He and Ne diffusion except to reduce anisotropy. Second, trapping is the primary mechanism altering He and Ne diffusion even at low dose, and we predict the maximal trapping energies for He and Ne to be 164 and 320 kJ/mol, similar to values inferred from experimental data. We also propose that the closure temperature increases non-linearly with damage, with effective trapping energy increasing with dose until a threshold, possibly corresponding to a percolation transition, after which retentivity decreases. Based on field data sets we also anticipate a value for this threshold of around ~2-5×10 17 /g, lower than previously proposed. We show Ne to be highly blocked by damage and predict similar diffusion behavior to He, but with higher retentivity. We demonstrate the importance of investigating rare gas diffusion at the atomic level for comparison with experimental data, in order to build a predictive diffusion law at different scales.

1.

Introduction

In the last two decades, radiogenic 4 He and nucleogenic 21 Ne have been used for (U-Th)/ 4 He and (U-Th)/ 21 Ne thermochronology in different minerals such as apatite, zircon and iron oxides for a broad range of geological applications (e.g. [START_REF] Farley | U-Th)/Ne and multidomain (U-Th)/He systematics of a hydrothermal hematite from Eastern Grand Canyon[END_REF][START_REF] Farley | U-Th)/He Dating of Phosphates: Apatite, Monazite, and Xenotime[END_REF][START_REF] Gautheron | U-Th)/Ne chronometry[END_REF][START_REF] Reiners | He diffusion and (U-Th)/He thermochronology of zircon: initial results from Fish Canyon Tuff and Gold Butte[END_REF][START_REF] Shuster | Weathering geochronology by (U-Th)/He dating of goethite[END_REF][START_REF] Zeitler | U-Th-He dating of apatite: A potential thermochronometer[END_REF]. Robust He and Ne age interpretation strongly depends on quantitative knowledge of diffusion behavior in crystals, which is linked to the closure temperature Tc [START_REF] Dodson | Closure temperature in cooling geochronological and petrological systems[END_REF], representing to first order the temperature below which He and Ne are significantly retained during cooling. He and Ne diffusion in crystals theoretically follows a simple Arrhenius law with isotropic or anisotropic diffusion. However, it has been shown that, for apatite and zircon, diffusion behavior is highly variable. Experimental data sets have Tc values ranging from 40 to 120°C for He in apatite [START_REF] Shuster | The influence of natural radiation damage on helium diffusion kinetics in apatite[END_REF] and from -60 to 190°C for zircon (Guenthner et al., 2013, and references therein). For Ne diffusion in zircon, data are sparse due to analytical difficulties, and Tc estimates come only from indirect assessment and one published abstract [START_REF] Gautheron | U-Th)/Ne chronometry[END_REF][START_REF] Van Soest | Neon diffusion in titanite and zircon and the utility of these minerals for (U-Th)/Ne thermochronometry[END_REF].

Radiation damage is considered to be a primary contributor to this variability, although the impacts of damage type (alpha, recoil and fission damage) and topology (Frenkel-pair, large damage and amorphous zone) on He diffusion in zircon are not well defined. In addition to forming traps for diffusing elements [START_REF] Farley | Helium diffusion from apatite: general behavior as illustrated by Durango fluorapatite[END_REF], damage can also change diffusive behavior for even low-damage zircon by reducing the strong anisotropy along the c-axis [START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF][START_REF] Farley | He diffusion systematics in minerals: Evidence from synthetic monazite and zircon structure phosphates[END_REF]. Similar complications may arise from other imperfections such as inclusions, voids, and defects [START_REF] Zeitler | Solubility and trapping of helium in apatite[END_REF][START_REF] Danišík | Seeing is believing: Visualization of He distribution in zircon and implications for thermal history reconstruction on single crystals[END_REF].

The relationship between He diffusivity and dose in zircon suggests another topological effect. Diffusivity has been shown to decrease with increasing dose until a threshold is reached, after which it strongly increases with further damage [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF], although recent work suggests that this threshold dose may vary [START_REF] Anderson | Empirical constraints on the effects of radiation damage on helium diffusion in zircon[END_REF]. This flip in behavior may be linked to the onset of percolation, in which overlapping damage zones start to form an interconnected network [START_REF] Salje | Is "metamictization" of zircon a phase transition?[END_REF], essentially dividing the crystal into two transport regimes, one through normal crystal (albeit with imperfections) and one through the damage network [START_REF] Ketcham | Geometric analysis of radiation damage connectivity in zircon, and its implications for helium diffusion[END_REF]. Superimposed on the complexity of diffusion is the process of annealing, which changes the landscape in which diffusion occurs, as it is happening. Recent work suggests that damage annealing in zircon also may depend on dose [START_REF] Ginster | Annealing kinetics of radiation damage in zircon[END_REF], further highlighting the challenges to understanding the regime in which diffusion operates.

As a result of this complexity, there are limits to how much progress can be made using exclusively empirical observations and experiments on incompletely-understood natural materials (e.g., [START_REF] Reiners | Zircon (U-Th)/He thermochronometry: He diffusion and comparison with 40 Ar/ 39 Ar dating[END_REF][START_REF] Reiners | Zircon (U-Th)/He thermochronometry[END_REF][START_REF] Cherniak | Diffusion of helium in zircon and apatite[END_REF][START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF][START_REF] Anderson | Empirical constraints on the effects of radiation damage on helium diffusion in zircon[END_REF][START_REF] Ginster | Annealing kinetics of radiation damage in zircon[END_REF][START_REF] Zeitler | Solubility and trapping of helium in apatite[END_REF][START_REF] Danišík | Seeing is believing: Visualization of He distribution in zircon and implications for thermal history reconstruction on single crystals[END_REF]. A rigorous physical understanding of the diffusion process and how it is affected by imperfections is an indispensible complement for the varied and extensive data gathering currently underway.

Although the inevitably chaotic and variable nature of damage zones makes them hard to characterize, transport of He and Ne in the areas between damage zones, and between damaged and undamaged crystal, is accessible through theoretical and computational approaches.

Several studies have been conducted using Density Functional Theory (DFT) calculations to quantify the He diffusion coefficient in damage-free zircon [START_REF] Reich | Low-temperature anisotropic diffusion of helium in zircon: Implications for zircon (U-Th)/He thermochronometry[END_REF]Saadoune and De Leeuw, 2009;Saadoune et al., 2009). Independent studies consistently find that He is not expected to be retained in the zircon structure at Earth-surface temperatures; for example, [START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF] calculate Tc ~ -150°C along the crystallographic c-axis.

However, recent DFT studies on He and Ne diffusion in apatite and hematite give results more similar to experimental data on real specimens, supporting the veracity and applicability of this theoretical approach (Balout et al., 2017a;[START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF][START_REF] Gerin | Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales[END_REF].

In this contribution, we investigate He and Ne diffusion in zircon from atomic to grain scales using a coordinated series of methodologies developed in recent studies [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF], Balout et al., 2017a, b;[START_REF] Gerin | Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales[END_REF] to examine important factors that can modify diffusion rates, including trapping, anisotropy, and path obstruction. As no published Ne diffusion data exist, we use DFT at the atomic scale to build and optimize the ideal zircon unit cell accounting for quantum features, and we inventory the insertion and interstitial sites. Once these sites are identified, we combine the DFT approach with the Nudged Elastic Band (NEB) method to determine the migration energies between them. We complete our atomic-scale characterization using Transition State Theory (TST) to compute jump rates between sites. To extend to a larger scale spanning defects, DFT and TST outputs are used in a Kinetic Monte Carlo (KMC) simulation to determine the diffusion trajectories and the effective diffusion coefficients in 3D. We then compare our results with published experimental data and investigate the processes that can modify He and Ne diffusion in zircon. We simulate and evaluate different hypotheses, such as trapping of He and Ne in damage, and diffusion pathway obstruction, in 3D, and their impact on closure temperature. Finally, we discuss the model implications for the (U-Th)/He and (U-Th)/Ne thermochronometers.

Methods, the multi-scale approach

We first provide a general summary of our methods and how they interact with each other, as illustrated in Figure 1. Density Functional Theory (DFT) is a computational modeling method to calculate the energy of interacting atoms by calculating their external electronic configuration. This is done by solving the Schrödinger equation for the electrons using functionals of the electron density. The DFT just solves static problems, atoms at given positions, and does not include dynamics, contrary to molecular dynamics. However, in contrast to molecular dynamics it includes the quantum effects for the electrons, thereby providing the most accurate energy estimation of an ensemble of interacting atoms. Due to the high calculation time, requiring cluster computing, the number of considered atoms remains small (<150 atoms) compared to a real crystal, and therefore a super-cell containing a small number of motionless atoms is used and an infinite repetition of this super-cell is assumed to represent an infinite lattice. The stable position of the atoms for a crystal lattice is found by nudging atom positions until the total energy is minimized. In the same way, an insertion site for He or Ne is found by inserting it into a pre-computed lattice and finding the positions with the minimum energy with respect to displacement after relaxation of the neighboring atoms.

He and Ne atoms can statistically jump between insertion sites when they vibrate as a consequence of the finite temperature. This is a dynamical process, beyond the scope of DFT, which addresses static configurations only. Therefore, another approach has to be followed, namely the Transition State Theory (TST), which computes the jumps statistically. It is based on a full mapping of the energy field encountered when the inserted atom moves along its path to its neighboring site, as delivered by the DFT analysis. This approach is valid because the Born-Oppenheimer approximation allows consideration of the motion of atoms as quasistationary from the electron viewpoint.

The ideal pristine crystal represented by DFT and TST does not capture all the properties of a natural one, which will always contain defects, distortions, dislocations and substitutions, and typically mineral and/or fluid inclusions. However, even in a highly degraded crystal, diffusion occurs in pristine lattice for at least some length scale. The approach we develop here is to consider the ideal case and see how it should be amended when alterations of the pristine lattice are added. Only when significant amorphisation takes place does this approach become invalid. The alteration of the ideal lattice by vacancies, interstitial atoms, and substitutions breaks the periodicity of the lattice and introduces heterogeneity over length scales much larger than the lattice cell dimension. Their effects on diffusion depend on their energy and their topology. Therefore, the easiest way to study them is the simulation of diffusion by random walk over many lattice cells spanning many defects. For this we use the Kinetic Monte Carlo method (KMC), which operates at the mesoscopic level and delivers diffusion coefficients for a heterogeneous medium, based on the energetics of the DFT and the transition rates of the TST (Figure 1). Finally, to compute He and Ne ages in crystals the diffusion equation has to be solved at the macroscopic level using the diffusivities delivered by KMC, or directly by DFT and TST in case of an ideal crystal.

In summary, DFT calculates the mapping of energy as a function of atom location. TST then computes the transition rates by using this mapping and by making use of thermodynamic equilibrium in phase space. KMC then quantifies the effect of randomly distributed alterations of the crystal-like traps or obstructions, independently of the mechanism producing them (interstitial atoms, substitutions, etc.). This pipeline delivers bulk diffusion coefficients for solving the 3D diffusion equation for a finite grain, as illustrated in Figure 1.

Density Functional Theory

Following the approach described in detail by [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF], we use the Periodic Density Functional Theory (DFT; Hohenberg and Kohn, 1964;[START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF], as implemented in the Vienna Ab-initio Simulation Package (VASP) [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | Ab initio molecular dynamics for liquid metals[END_REF], for all structural optimization and insertion site investigations of both He-and Ne-doped zircon lattices. All parameters were optimized on the zircon unit cell, and the relaxed unit cell geometrical parameters are in good agreement with experimental values. To minimize the volume relaxation effect due to the incorporation of Ne and He atoms, and to avoid interactions between He or Ne with neighboring cells, the zircon unit cell was duplicated along each spatial direction, resulting in a 2×2×2 super cell consisting of 192 atoms (Figure 2A and 2B). We determine the insertion energy of He and Ne in each interstitial site by relaxing the structure with the help of the gradient technique at constant volume. The insertion energy is defined as the energy needed to bring the He or Ne from outside the crystal into the interstitial site. It is calculated as:

𝐸 𝑖𝑛𝑠 = 𝐸 𝐻𝑒+𝑧𝑖𝑟𝑐𝑜𝑛 -(𝐸 𝑧𝑖𝑟𝑐𝑜𝑛 + 𝐸 𝐻𝑒 ) (Eq. 1)
where 𝐸 𝐻𝑒+𝑧𝑖𝑟𝑐𝑜𝑛 , 𝐸 𝑧𝑖𝑟𝑐𝑜𝑛 and 𝐸 𝐻𝑒 are respectively the absolute energy of the system with He or Ne in the zircon host lattice, the absolute energy of the host lattice, and the energy of the isolated He or Ne atom.

We use the Nudged Elastic Band (NEB) [START_REF] Jónsson | Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions[END_REF][START_REF] Mills | Reversible work transition state theory: application to dissociative adsorption of hydrogen[END_REF] method to determine the migration energy between insertion sites. NEB is an efficient approach to estimating the minimum energy path (MEP) and the migration energy between two interstitial sites (𝐸 𝑚𝑖𝑔 𝑖𝑗 ), which is defined as the difference in energy from when the atom is at the saddle point to when it is in the insertion site. It is also called the barrier energy. The NEB is a chainof-states method [START_REF] Elber | A method for determining reaction paths in large molecules: Application to myoglobin[END_REF][START_REF] Pratt | A statistical method for identifying transition states in high dimensional problems[END_REF] found by constructing a series of system images between states. Adjacent images are linked by a spring force to ensure the continuity of the path, thus mimicking an elastic band. In this study, all NEB calculations are performed with a spring constant of -5 eV/Å 2 .

Transition State Theory

Interstitial diffusion occurs through a random walk of an atom from one site in the crystal lattice to another neighboring free one. The atom needs enough energy to accomplish this motion; as there is an energy barrier to overcome, diffusion is a thermally activated process.

As already mentioned, this dynamical process cannot be directly calculated by DFT, which addresses static configurations. However, as the motion occurs in thermodynamic equilibrium the canonical ensemble is assumed in phase space. This leads to the transition state theory [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF]Voter, 1986;[START_REF] Wert | Interstitial Atomic Diffusion Coefficients[END_REF], in which the atomic motion is controlled by the jump probabilities, Γ, which depend exponentially on temperature.

𝛤 𝑖𝑗 = 𝜈 0 𝑒 - 𝐸 𝑚𝑖𝑔 𝑖𝑗 𝑘 𝐵 𝑇 (Eq. 2) With 𝜈 0 = ∏ 𝜈 𝑖 𝑖𝑛𝑖𝑡 3𝑁 𝑖=1 ∏ 𝜈 𝑖 * 3𝑁-1 𝑖=1 Vineyard approx. ⇒ ≈ ∏ 𝜈 𝑖 𝑖𝑛𝑖𝑡 3 𝑖=1 ∏ 𝜈 𝑖 * 2 𝑖=1 (Eq. 3)
And 𝐸 𝑚𝑖𝑔 𝑖𝑗 = 𝐸 * -𝐸 𝑖𝑛𝑖𝑡 (Eq. 4)

Where 𝜈 0 is the attempt frequency; 𝜈 𝑖 𝑖𝑛𝑖𝑡 and 𝜈 𝑖 * are the normal mode frequencies at, respectively, the initial and transition states; 𝐸 𝑚𝑖𝑔 𝑖𝑗 represents the migration energies; and N is the number of atoms in the system. For a molecule with N atoms there are 3N degrees of freedom. For a non-linear molecule three degrees of freedom can be assigned to translation of a body as a whole (Tx, Ty, Tz) and three to rotations (Rx, Ry, Rz). The remaining motions of the atoms are displacements of the atoms from their mean positions. These fundamental vibrations are referred to as "normal modes". Thus, a non-linear molecule has 3N-6 normal modes. In the insertion location, the energy is minimum, and therefore the energy has an upward curvature for any displacement of the inserted atom. This 3D-curvature can be expressed with three eigenvalues, which translate directly to the three normal frequencies. This translation is exactly the same as for any harmonic oscillator like a pendulum. At the saddle point, the same features apply, but only transverse motion relative to the transition path is involved, and therefore there are two frequencies. According to the [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF] approximation, for one atom (He or Ne), the energetics of other atoms are not affected by the position of the migrating atom being considered. The TST allows thus provides a well-understood theoretical approach for calculating the frequency factor (Figure 1).

Kinetic Monte Carlo Method

In a pristine lattice, the activation energy (Ea) corresponds to the migration energy obtained from DFT, and the frequency factor (D0) can be extracted directly from the TST formulas (Fig. 1). However, when obstructions occur the diffusivity depends on their topology, and a simulation over macroscopic crystals has to be run to determine the modified diffusion coefficients. The Kinetic Monte Carlo (KMC) method [START_REF] Bortz | A new algorithm for Monte Carlo simulation of Ising spin systems[END_REF]Gillespie, 1976) is particularly apt for our purposes because it can be used to extend simulation to time scales far beyond the vibrational period and lengths much larger than the cell size. The KMC simulates, at a given temperature, a random walk over the interstitial sites based on the microscopic jump probabilities (Γ), which have been calculated from TST using the DFT results (Fig. 1). In addition, it takes into account potential obstructions and traps, as well as anisotropy.

For a selected temperature, one He (or Ne) atom is placed randomly in one available insertion site. For the possible diffusion directions, with the calculated jumping frequencies given by the TST, we randomly choose a number to simulate the jumping time for each direction. The distribution of jumping times follows an exponential decay. The jump direction with the lowest time is chosen, the time and coordinate are updated, and the time is added to the total travel time. This is repeated m times for each atom, and n trajectories (atoms) are simulated, using the method of [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF]. For each trajectory, the total time is computed as the sum of the residence times in each site between jumps. The diffusion coefficient is obtained by averaging over the ensemble of n trajectories. If <x 2 >, <y 2 >, <z 2 > are the averages over the trajectories of the squares of distances between initial and final points in

x, y and z, <t> is the average time of the trajectories, and if we denote as Dx, Dy, Dz the diffusion coefficient along the axes, then according to the Einstein relation.

𝐷 𝑥 = <𝑥 2 > 2<𝑡> 𝐷 𝑦 = <𝑦 2 > 2<𝑡>
(Eq. 5)

𝐷 𝑧 = <𝑧 2 > 2<𝑡>
From the above relations the diffusion coefficients can be assessed even in the case of anisotropic diffusion, as is necessary for zircon. The temperature dependence is obtained by running the random walk at different temperatures.

To account for highly anisotropic behavior along the c axis and its possible obstruction, the random walk operates on a sub-lattice of 1281284096 interstitial sites. When an atom jumps against an edge of the sub-lattice it is automatically re-inserted into the opposite side, so that the geometry can be considered as infinite with a periodicity given by the size of the sublattice. Obstructions are simulated by inhibiting some fraction of sites randomly distributed over the sub-lattice.

3D diffusion modeling

Finally, these two descriptions provide a bulk diffusion coefficient, which is used to calculate the diffusion over a finite grain where the surface is taken into account by solving in 3D the diffusion equation appropriate to surface boundary conditions. A modified version of the rare gas 3D diffusion code published in [START_REF] Gautheron | A Monte Carlo approach of diffusion applied to noble gas/helium thermochronology[END_REF] and used in subsequent studies [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF][START_REF] Gautheron | Effect of alpha-damage annealing on apatite (U-Th)/He thermochronology[END_REF]2012) has been implemented to simulate He and Ne diffusion in zircon. This code solves the diffusion equation by a Monte Carlo method at the macroscopic level, using the diffusivity extracted from the KMC simulation. This approach accounts for crystal surface boundaries, which is not possible in the KMC code.

Results

Investigation of He and Ne interstitial sites

Zircon crystallizes as a body-centered tetragonal unit cell, with a structure consisting of a chain of alternating edge-sharing (SiO4) tetrahedral and ( ZrO8 The interatomic He -O distances have three different values of 2.02, 2.22 and 2.36 Å. For the Ne insertion case, the Ne atom is located 2.16, 2.30 and 2.42 Å from neighboring oxygen atoms.

These larger distances are due to the larger size of the Ne atom, making the relaxation effect important for Ne.

The insertion energies calculated based on equation 1 for He and Ne in the different insertion sites are listed in Table 1. The derived helium insertion energy is 1.46 eV, which is in good agreement with other DFT calculations performed by Saadoune et al. (2009) and for comparison obtained insertion value in hematite is given (Balout et al., 2017a, b). In the case of Ne, the insertion energy of 3.10 eV is almost twice that calculated for He. This is related to the larger atomic radius of Ne, which must repel the neighboring atoms more to enter the interstitial site.

Migration Energies (Emig)

In the c-channel, since the inter-site distance S1 -S1' is about 1.5 Å, only one image is taken between the two sites for the NEB calculation. For NEB calculation in this channel we start with a linear path composed of a series of 4 images connecting the two S1 (S1') sites. After the relaxation process we found that the migration energy between the sites separated by the rectangular gap is equal to 0.66 and 1.53

eV for He and Ne atoms, respectively (Table 2; Figure 3B). For the rhombic gap, we find a migration energy of 2.60 eV for He, but for Ne the NEB calculation indicated that no migration is possible (Figure 3B'). Figure 4 illustrates the sequence of the interstitial sites along c and how they are connected to the pathways along a and b.

The paradoxical result is that along c the energy barrier for Ne is lower than for He, even though its radius is larger as already recognized in the insertion energy calculation. This can be understood by noting that the height of the barrier is a difference in energy. For Ne both the energy in the interstitial site and at the top of the barrier are higher in comparison to He, reflecting the atom radius effect, but this does not say anything about the difference. As a conclusion, contrary to the insertion energy, the activation energy is not well controlled by the ionic or atomic radius. This outcome was already evident in [START_REF] Farley | He diffusion systematics in minerals: Evidence from synthetic monazite and zircon structure phosphates[END_REF].

According to our results, one can say that, at low temperature, in an ideal crystal He and Ne will diffuse mainly along the c-channel. However, as temperature increases, there will be competition between the probability of jumping the 0.24 eV barrier and the 0.66 eV barrier in the case of He diffusion. For Ne-diffusion, however, even at high temperature diffusion along the c-axis remains strongly favored.

Ideal lattice results

Using Transition State Theory as presented in section 2, the attempt frequencies ν, which allow determination of relative jump probabilities, were calculated for each jump case.

All calculated jump probabilities are provided in Table 3. From the parameters listed in this table, the diffusion coefficient along the c-axis can be derived analytically for a pristine crystal:

𝐷 𝑧 = 𝜈. 𝑎 2 . exp (- 𝐸 𝑚𝑖𝑔 𝑘 𝐵 𝑇 ) (Eq. 6)
with a=1.5 Å. In the same way the diffusivity along a and b for an undamaged crystal and for practical temperatures (T<10 4 K) is obtained as:

𝐷 𝑥𝑦 = 𝜈.𝑎 2 4 . exp (- 𝐸 𝑟𝑒𝑐𝑡 𝑘 𝐵 𝑇 ) (Eq. 7)
with a=3.35 Å. The factor 4 accounts for the fact that only one site over four allows a jump in a given direction (Fig. 4), as crossing the rhombic diamond is quasi-forbidden.

We get respectively for He and Ne:

𝐷 𝑧,𝐻𝑒 = 1.6 × 10 -7 𝑒𝑥𝑝 (- 0.24 𝑘 𝐵 𝑇
) 𝑚 2 𝑠 ⁄ (Eq. 8)

𝐷 𝑥𝑦,𝐻𝑒 = 1.5 × 10 -7 𝑒𝑥𝑝 (-0.66

𝑘 𝐵 𝑇
) 𝑚 2 𝑠 ⁄ (Eq. 9)

𝐷 𝑧,𝑁𝑒 = 1.2 × 10 -8 𝑒𝑥𝑝 (-

0.22 𝑘 𝐵 𝑇
) 𝑚 2 𝑠 ⁄ (Eq. 10)

𝐷 𝑥𝑦,𝑁𝑒 = 1.1 × 10 -7 𝑒𝑥𝑝 (-

1.53 𝑘 𝐵 𝑇
) 𝑚 2 𝑠 ⁄ (Eq. 11)

From these results one can conclude that Ne diffusivity is slightly higher at low temperature than He diffusivity in a perfect zircon host lattice. In addition, the diffusivity is highly anisotropic and diffusion occurs essentially along the c-axis. The anisotropy is even stronger in the Ne case due to the high energy barrier when crossing the rectangular gap.

Discussion

He and Ne diffusion in non-damaged zircon structure

Comparison to previous calculations

Helium diffusivity in zircon, calculated by DFT and molecular dynamics, has already been reported [START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF][START_REF] Reich | Low-temperature anisotropic diffusion of helium in zircon: Implications for zircon (U-Th)/He thermochronometry[END_REF]Saadoune et al., 2009), and for comparison the insertion energy and barriers are listed respectively in Tables 1 and2 2), for a reason which is not elucidated. We will see later that this is barrier is the one that mainly controls the diffusivity in zircon in practical situations.

Concerning He diffusivity, differences exist with respect to the other calculations beyond the migration energy. [START_REF] Reich | Low-temperature anisotropic diffusion of helium in zircon: Implications for zircon (U-Th)/He thermochronometry[END_REF] and Saadoune et al. (2009) used a 1D formula which is known to be inadequate because it assumes that the energetic shape around the interstitial site and the saddle-point is the same. However, this assumption only acts on the preexponential factor and its impact is limited. [START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF] did not use this assumption but assumed that along a/b the higher barrier (~2.6 eV) dominates, disregarding the relatively easy intermediate jump along c that can bypass this barrier. This explains why they obtained a very low diffusivity along a/b.

Comparison to experimental data

For Ne only one abstract reports diffusion data, with an Ea value of ~339 kJ/mol and closure temperature of 396-410°C for a 10°C/Ma cooling rate [START_REF] Van Soest | Neon diffusion in titanite and zircon and the utility of these minerals for (U-Th)/Ne thermochronometry[END_REF]). This value is much higher than that computed in this work. Similarly, He diffusion data obtained on natural zircon crystals present significantly higher apparent activation energies than obtained by theoretical studies. In addition, measured (U-Th)/He and (U-Th)/Ne zircon ages confirm that He and Ne are strongly retained in this mineral. For illustration, Figure 5 presents some He diffusion data obtained on natural zircon crystals (e.g., [START_REF] Reiners | Zircon (U-Th)/He thermochronometry: He diffusion and comparison with 40 Ar/ 39 Ar dating[END_REF][START_REF] Cherniak | Diffusion of helium in zircon and apatite[END_REF][START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF] compared to our data for a pristine zircon. In detail, Figure 5 shows the distribution of activation energy, frequency factor and closure temperature for He obtained for natural zircon, and the deduced activation energy for Ne needed to explain a (U-Th)/Ne age [START_REF] Gautheron | U-Th)/Ne chronometry[END_REF]. The authors based their estimate on the known thermal history of the Gold Butte, Nevada area [START_REF] Reiners | Helium and argon thermochronology of the Gold Butte block, south Virgin Mountains, Nevada[END_REF] should thus be greater than the one for He (i.e., >170 kJ/mol), as Ne appears much more strongly retained than He (Figure 5), in contradiction with the result of our calculation in which a larger diffusivity is found along c. Based on the linear relationship between Ea and closure temperature, for Ne we can anticipate an activation energy value above ~170 kJ/mol, but we cannot clearly define the upper limit, and the lower limit is below the value determined by van [START_REF] Van Soest | Neon diffusion in titanite and zircon and the utility of these minerals for (U-Th)/Ne thermochronometry[END_REF].

Processes modifying He and Ne diffusion

A strong discrepancy exists between results obtained using theoretical and experimental approaches, leading to a misinterpretation of the validity of the theoretical method. However, recent DFT studies on He and Ne diffusion in apatite [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF] and hematite (Balout et al., 2017a;[START_REF] Vacherat | Rift-to-collision transition recorded by tectono-thermal evolution of the northern Pyrenees[END_REF], which produce expectations consistent with experiments, bring a strong credit to this theoretical approach. If the discrepancy is not due to a methodological issue, additional features and processes are likely to be modifying He and Ne diffusion in zircon, as already discussed by [START_REF] Reich | Low-temperature anisotropic diffusion of helium in zircon: Implications for zircon (U-Th)/He thermochronometry[END_REF], Saadoune et al. (2009) and [START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF].

Until now, our calculation has only considered an ideal pristine crystal, but in real crystals defects are present and change local diffusion properties. The defects can be substitutions in which an atom of the ideal lattice (Zr usually) is exchanged for another, such as a rare earth element. Another type is Frenkel pairs, consisting of an atom knocked out from its original lattice site, leaving a vacancy, and inserted somewhere else. Frenkel pairs are thermally activated and exist at any temperature, meaning that they can be created and recombined by thermal motion. They are highly enhanced by irradiation due to the decay of U and Th. The first effect of the alpha particle, and more strongly, the recoiling daughter is to create Frenkel pairs by knocking atoms out of the lattice. For high radiation doses the density of vacancies and interstitial atoms becomes high enough to make the lattice unstable, collapsing it into amorphized zones. The diffusivity in such zones is not known and only hypotheses can be proposed. Dislocations are yet another likely pervasive defect type that can alter diffusional pathways. The detailed description of these defects is far beyond the scope of this work.

Similarly, meso-scale features that may alter diffusivity and retention, such as mineral and fluid inclusions, are also beyond our scope. Our methodology is based on a simple modeling of their effect on the diffusion of inserted He and Ne atoms, by considering that they diffuse in an ideal crystal (Figure 6A) until they find some defect altering the diffusion (Figure 6B and 6C). We address only two phenomena, obstruction and trapping.

Obstruction may occur when an interstitial atom, knocked out due to irradiation for example, is already in an interstitial site, preventing He or Ne from entering the site. This mechanism is efficient for slowing down diffusion along fast pathways along the c-axis in zircon, because it requires the atom to jump over higher barriers to span long distances. The motion is illustrated in Figure 6B, whereas Figure 6A shows the ideal case where the jumps occur primarily along the c-axis due to the low barrier. Obstruction can also result from chemical substitutions that alter the energetics along the diffusion pathway. This has been shown to occur for Cl-rich apatite compared to F-apatite [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF], although the authors demonstrate that a significant effect on diffusion only occurs at a relatively high substitution level (some dozens of percent of the halogen sites) in that mineral. For zircon, Hf substitution can be quite important (some wt% percent, e.g., [START_REF] Hanchar | Rare Earth element behavior in zircon-melt systems[END_REF] and will modify He diffusivity. Dislocations could also lead to obstruction by a local shrinking of the original cell. In Figure 6B for example, the He atom first diffuses along the c-axis and crosses the rhombic gap but is then blocked due to the obstruction. Finally, the He atom will go back as it will demand too much energy to cross the rectangular gap, where the migration energy is 2.6 eV (251 kJ/mol). As already emphasized in [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF], obstruction is very efficient in reducing diffusion anisotropy, which is an important feature in the case of zircon.

The other effect we consider is trapping. This phenomenon has been already considered in several works [START_REF] Farley | Helium diffusion from apatite: general behavior as illustrated by Durango fluorapatite[END_REF][START_REF] Shuster | The influence of natural radiation damage on helium diffusion kinetics in apatite[END_REF][START_REF] Gautheron | Effect of alpha-damage annealing on apatite (U-Th)/He thermochronology[END_REF][START_REF] Flowers | Apatite (U-Th)/He thermochronology using a radiation damage accumulation and annealing model[END_REF][START_REF] Shuster | The influence of artificial radiation damage and thermal annealing on helium diffusion kinetics in apatite[END_REF][START_REF] Gerin | Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales[END_REF] and can result from different kinds of defects.

When vacancies are created, in particular under irradiation, they offer more space to the diffusing atom at lower energy. Therefore additional energy, noted (EaT), is necessary to extract it from its location and bring it to a normal interstitial site. In general, trapping may occur in any situation where the interstitial atom has more space to be accommodated. The amorphized zones could offer this more open space and serve as He and Ne storage. If the zones are not interconnected, diffusion still occurs through an ideal lattice endowed with trapping volumes. Even substitutions in some cases can act as traps. Saadoune and de Leeuw (2009) demonstrate that when U 4+ and Pu 4+ are present at the percent level in the zircon lattice in the vicinity of an He insertion site, U and Pu increase He solubility by decreasing He insertion energy, leading to a trapping effect.

Radiation damage is the favored mechanism producing trapping because it creates vacancies, and in the case of zircon amorphized zones (Figure 6C). This mechanism was initially proposed for apatite, where radiation damage strongly influences He diffusion by trapping He inside damage [START_REF] Farley | Helium diffusion from apatite: general behavior as illustrated by Durango fluorapatite[END_REF][START_REF] Shuster | The influence of natural radiation damage on helium diffusion kinetics in apatite[END_REF]. As a guideline, the maximum available trapping energy corresponds to a case where the diffusing atom is far from its neighbors, such as in a cavity (large damage, voids, fluid inclusions). Any short distance to its neighbors will repel it and increase its energy, thereby reducing (EaT). The extreme case where the atom is far from any other atom is equivalent to an atom outside of the crystal, and in this case the trapping energy is merely the insertion energy. Therefore, the maximal trapping energy is the insertion energy. The only exception to this rule would occur if atoms in the crystal would interact with the He or Ne by attracting it. [START_REF] Gerin | Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales[END_REF] described the energetics of this process for apatite but similar conclusions can be drawn for zircon. As the insertion energy is 1.46 and 3.1 eV for He and Ne respectively, one can infer that the maximal trapping energy for large damage such as recoil damage (EaT) will be respectively 1.46 eV (with global Ea 1.7 eV = 1.46+0.24 eV; 164 kJ/mol) and 3.31 eV (with global Ea 3.32 eV = 3.10+0.22 eV; 320 kJ/mol). In crystals with real damage, different geometries of defects are present, entailing a distribution of EaT. In addition, when radiation damage accumulates vacancies can connect creating larger spaces, increasing the trapping energy. Figure 6C illustrates a diffusion trajectory with different topologies of damage, with small damage (such as point defects), single recoil damage or recoil damage clusters. The EaT will be different for each case according to the available free space. An important feature of trapping is that, unlike obstruction, it has no effect on the anisotropy of diffusion, so that it cannot be responsible for the reduction of the anisotropy in case of zircon, although it has the capability to reduce diffusivity strongly, and may be spatially associated with other defects.

Impact of blocking pathways, i.e. obstruction, on He and Ne diffusion in zircon

It is expected that atoms displaced into the c-axis channel as a result of radiation damage strongly affects diffusion along c, which is quasi-free at normal temperature in an ideal lattice.

In this section, we test the obstruction model using the KMC method and discuss the impact of obstructed pathways on He and Ne diffusion coefficients, by using the model described in section 2.3; additional details are given in Appendix A. especially at low temperature where the slope of the diffusion behavior changes significantly (Figure 7B). An important result is that, surprisingly, the impact of obstruction on 3D diffusion is higher at low temperature with an important deviation from the ideal lattice. As the fraction of obstruction increases the dependence moves closer to the a/b diffusivity. For T<300°C all the curves are almost straight lines parallel to the a/b Arrhenius law, implying similar activation energy, but with D0 scaling down as the level of obstruction increases.

As soon as the latter exceeds 1%, He diffusion along the c-axis is significantly blocked leading to a change in diffusivity. One particular outcome is that the anisotropy, which was very strong in the ideal lattice, is reduced by a factor depending of the level of obstruction, but almost independent of the temperature. For example, at 10% obstruction the anisotropy drops down to 20. It is worth noting that the above key trends are observed experimentally by [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF], who show that He activation energy orthogonal to c is only slightly affected by damage, whereas with increasing damage diffusion along c approaches a constant Ea close to that of the transverse direction. This experimental effect is reproduced by our calculation. It reduces the anisotropy, compared to ideal crystal, and cannot be explained by a trapping effect alone. However, our modeled activation energy in any direction remains much smaller (64 kJ/mol) in comparison to those measured in natural zircon [START_REF] Reiners | Zircon (U-Th)/He thermochronometry: He diffusion and comparison with 40 Ar/ 39 Ar dating[END_REF][START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF], which are around 169 kJ/mol.

In conclusion, although obstruction strongly slows down He diffusion along c, it is not the primary parameter that drives retentivity of He in natural zircons. However, obstruction is the main parameter acting on the anisotropy of diffusion, reducing it strongly with respect to the ideal lattice case.

For Ne the change of diffusivity with obstruction is not progressive but total. At obstruction levels as low as 1%, Ne diffusion is already frozen because the activation energy along the a-and b-axes is so high (147 kJ/mol, Table 3). Nevertheless, as with He, diffusion pathway obstruction alone cannot reproduce the inferred natural Ne diffusion behavior of zircon; as described in section 4.2, an activation energy >180 kJ/mol is necessary to reproduce Ne ages [START_REF] Gautheron | U-Th)/Ne chronometry[END_REF], well above the maximum activation energy of ~147 kJ/mol obtainable from obstruction. Thus, for He and Ne, pathway obstruction mechanisms are not sufficient to block He and Ne diffusion in zircon to the degree observed in natural specimens.

Impact of traps on He and Ne diffusion in zircon

In this section, the trapping model is tested and the amount of trapping energy necessary to reproduce the inferred natural He and Ne diffusion coefficients is discussed. For this purpose, the model published in [START_REF] Gerin | Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales[END_REF] is used where the diffusion coefficient modified by damage, denoted 𝐷 ̃, follows a formula (eq. 12) that includes a trapping term:

(Eq. 12)

Where D(T) is the diffusion coefficient of undamaged zircon (data from this study), f is the damage fraction ranging from 0 to 1, representing the normalized amount of displaced atoms in the structure, and EaT is the additional activation energy necessary for an He or Ne atom to jump out of the trap to the next insertion site. Appendix B provides a complete explanation and demonstration of the equation, and the conversion between damage fraction and alpha dose. In this study, the trapping model is correct only if damage stays as an isolated trap and does not form a percolating network connecting to the grain surface, as then it would not be a trap anymore.

The evolution of the simulated He closure temperature using Eq. 12 to account for equivalent alpha dose is reported in Figure 8A, for a 60-µm radius spherical grain. The simulation was done for Ea values of 40, 80, 120 and 140 kJ/mol, where 140 kJ/mol is the maximum predicted value from the insertion energy (EaT) (see section 4.2). One can observe from Figure 8A that a minimum Ea of 120 kJ/mol is necessary to reproduce retentive behavior in low-damage zircon of [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF]. Interestingly, the activation energies and closure temperatures reported in Figure 8B are similar to those reported for natural samples in For He, closure temperature will increase with the equivalent alpha-dose for relatively low trapping energy Ea (<40 kJ/mol). It seems paradoxical that the closure temperature depends only slightly on the damage fraction, as illustrated in Figure 8A, but this is well understood as this fraction shows up in the pre-exponential factor f in Eq. 12, whereas the activation energy is in the exponential, giving rise to a much stronger effect. In other words when the trapping energy is high most of the helium is trapped, whatever the number of traps provided it is higher than the number of helium atoms, and only a high temperature (some hundreds of °C during geological time) can extract them from the traps. This effect has been seen with the DFT calculation in zircon (Saadoune and De Leeuw, 2009) and apatite [START_REF] Gerin | Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales[END_REF], yet the damage topology between the two minerals is different. However, in the case of zircon a higher trapping energy is needed to reproduce experimental data with respect to expectations and calculations for apatite (~30-40 kJ/mol for apatite [START_REF] Gerin | Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales[END_REF][START_REF] Shuster | The influence of artificial radiation damage and thermal annealing on helium diffusion kinetics in apatite[END_REF] compared to 120-140 kJ/mol for zircon). The difference can be explained by the higher insertion energy for He and even higher for Ne (Table 1) in zircon compared to apatite (~0.6 to 1 eV; [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF], as modification of the crystal structure during damage creation will produce a site more energetically favorable compared to the lattice. In this case, the insertion energy strongly drops, leading to an increase of the trapping energy until a threshold EaT. The higher necessary trapping energy can be explained by the damage topology in zircon, as the crystal structure behaves quite differently than apatite. It is well known that natural apatites are almost always crystalline in contrast to zircon, which is often reported to be partially or fully metamict (e.g. [START_REF] Ewing | Radiation damage in zircon[END_REF]; additional infomation on damage can be found in Appendix B.

Damage percolation threshold

For zircon, at high dose, several authors have reported that damage overlap or percolation may create amorphous zones (e.g. [START_REF] Ewing | Radiation damage in zircon[END_REF]Pidgeon, 2014;[START_REF] Salje | Is "metamictization" of zircon a phase transition?[END_REF]. The more recent amorphous-crystal threshold has been estimated to be at ~2.210 18 /g (Pidgeon, 2014), which has been interepreted as the recoil percolation damage threshold [START_REF] Trachenko | Large swelling and percolation in irradiated zircon[END_REF]. The transition to rapidly rising diffusivity occurs at ~210 18 /g in the diffusivity data set used by [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF], although data by [START_REF] Anderson | Empirical constraints on the effects of radiation damage on helium diffusion in zircon[END_REF] suggest a lower value at 2-510 17 /g. In addition, [START_REF] Ketcham | Geometric analysis of radiation damage connectivity in zircon, and its implications for helium diffusion[END_REF] calculated that alpha recoil and fission damage percolate at different thresholds than previously assumed (~210 16 and ~210 17 /g for alpha-recoil damage in 6-8-step decay chains and isolated decays (nochain) respectively, and ~210 18 /g for fission tracks), due to the elongate nature of recoil and fission damage not taken into account by [START_REF] Trachenko | Large swelling and percolation in irradiated zircon[END_REF].

Recent ZHe data measured in Hercynian granites from the Pyrenees [START_REF] Bosch | Timing of Eocene-Miocene thrust activity in the Western Axial Zone and Chaînons Béarnais (west-central Pyrenees) revealed by multi-method thermochronology[END_REF][START_REF] Vacherat | Rift-to-collision transition recorded by tectono-thermal evolution of the northern Pyrenees[END_REF] exhibit age versus maximum alpha dose correlations that suggest that this retentivity maximum may be reached at ~2-510 17 /g (Figure 9). It is not straightforward to estimate the effective alpha dose needed to explain an age-eU pattern, so we have provided two solutions. First, as an end-member over-estimate, we simply retain all damage for each grain throughout its entire history since crystallization at 305 Ma (U/Pb age); the dose can be no higher than this (Figure 9A). Second, we posit that the divergence of measured ages was due to a reheating and partial resetting event at some time in the past that caused more He loss in the less retentive grains, and so we should consider only the possible damage production prior to that event. A conservative estimate for its timing is the age of the youngest grain, leading to the relationship shown in Figure 9B. Even this dose is likely to be a significant over-estimate of the pertinent dose for explaining the spread in ages, as it presumes that after crystallization damage retention began immediately.

The implied transition dose from increasing to decreasing He retentivity is significantly lower than to the one considered in the RDAAM model of [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF] and is closer to the no-chain alpha recoil percolation level calculated by [START_REF] Ketcham | Geometric analysis of radiation damage connectivity in zircon, and its implications for helium diffusion[END_REF]. The no-chain case is a possibly more appropriate approximation for connectivity that leads to fast pathways, as the connectivity due to U and Th decay chains is "tip to tip", and thus dominated by potentially obstructed regions at the boundaries of tracks. The dose range over which the nonchain percolating network expands to encompass the majority of recoil damage roughly coincides with the negative dose-age correlation in the Pyrenean samples, in the range 2 to 510 17 /g, significantly below the dose range where the Guenthner et al. ( 2013) model postulates a fall in retentivity (above ~2×10 18 /g). Such a lower alpha dose threshold for percolation can also explain why samples in old and stable crustal blocks [START_REF] Reiners | Zircon (U-Th)/He thermochronometry[END_REF] often displays only negative correlations between ZHe and eU (a proxy used for the alpha dose).

In such cases the maximum potential retentivity that connects negative and positive correlations could not be documented because it would correspond to very low and atypical U concentrations for zircon.

We infer from the Pyrenean samples that the recoil damage threshold will be at a value ranging around ~2-510 17 /g, as already proposed by [START_REF] Anderson | Empirical constraints on the effects of radiation damage on helium diffusion in zircon[END_REF], far below the ~210 18 /g proposed by [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF].

Implication for zircon (U-Th)/He and (U-Th)/Ne thermochronometry

He and Ne diffusion coefficients in 3D obtained in this study and the investigation of mechanisms that can modify He and Ne diffusivity demonstrate that obstruction and trapping alter diffusivity significantly (Figure 10). Based on atomic investigation of He diffusion energetics, we can predict how He will diffuse in the zircon structure. Figure 10A provides a schematic illustration of a diffusion pathway as postulated by our work. The He atom will diffuse primarily along the c-axis following a random diffusion pathway. He diffusion is facilitated along the c-axis because of the low migration energy; obstruction could slightly change the trajectory but the He atom can diffuse easily along the a-axis to get around any obstructed pathway. Following the example pathway of Figure 10A, the He atom is trapped in a large defect (i.e. recoil damage, voids, fluid inclusion), with a very low insertion energy (potentially approaching 0). The He atom will need a maximum activation energy of 1.7 eV (<164 kJ/mol) to go back to the crystal lattice (Fig. 10A). After that He diffuses again easily in the lattice from site to site along the a-axis and c-axis (passage through the rhombic gap is insignificant due to high migration energy, Fig. 4), until it reaches a second large energy well.

In this case the damage is connected to other damage extending to the crystal surface (Fig. 10A). In terms of energy, we can infer from our model that in the connected area, because the lattice is distorted, the migration energy may be lowered. He insertion energy is thus very low and the He atom will diffuse in the connected damage zone with lower energy barriers. In Figure 10A, we predict that the He migration in the damage ranges from 1.4 to 0.9 eV (135 to 87 kJ/mol). Indeed, several geological studies show that for high damage dose, ZHe ages are similar to AFT ages (see Fig. 9 for the Pyrenean case example), and in some cases are even younger or equivalent to AHe ages [START_REF] Johnson | Inverted" zircon and apatite (U-Th)/He dates from the Front Range, Colorado: High-damage zircon as a lowtemperature (<50°C) thermochronometer[END_REF]. This implies that at high dose (>10 18 /g), He retentivity is low in damage zone, with Tc of 110 to 50°C and thus Ea of 135 to 90 kJ/mol. We thus anticipate that, as opposed to the harmonic averaging model by [START_REF] Guenthner | red circles) and calculated with the DFT results (this study, red square) are reported in addition to the trapping model with trapping energies[END_REF], diffusivity in zircon after the percolation transition may be controlled by the energetics and topology of the damage network.

Figure 10B combines our modeling of trapping and obstruction with natural data in terms of the closure temperature evolution as a function of the alpha dose. First, obstruction of diffusion pathways reduces diffusivity, not as efficiently as trapping but possibly enough to play a role in addition to it. Indeed, at sufficient obstruction levels (> 20% for He and 1% for Ne), diffusion for He and Ne will be only possible when it includes a component of motion along the a-and b-axis (Table 3). The resulting diffusion behavior will still not be retentive enough to reproduce natural data, however. He and Ne global activation energy Ea values are respectively 63 kJ/mol and 147 kJ/mol, while values of 120-140 kJ/mol and >180 kJ/mol are needed (Figure 10B for He). The 20% obstruction level corresponds to 13% of atoms displaced (and so to an equivalent dose of ~10 18 /g) as explained in Appendix B, and corresponds to an equivalent closure temperature of approximately -75°C assuming a D0 of 1.510 -7 m 2 /s. This value is not very tightly constrained, but provides a reasonable, illustrative first estimate of He retention.

If our modeling is correct, it is possible to make some predictions. First, for He and Ne, a quite high trapping energy is necessary to reproduce closure temperatures inferred from natural samples (Figure 8 and Figure 10B), and we predict maximum possible trapping energies with EaT of 140 kJ/mol (Emig=1.46 eV) and 300 kJ/mol (Emig=3.10 eV) for He and Ne respectively. These trapping energies (EaT) are within the limits provided by our DFT calculations, as discussed in section 4.2. It means that the total energy needed to go back in the crystal structure will be 164 kJ/mol (Einsertion + Emig=1.7 eV) and 320 kJ/mol (Einsertion + Emig=3.32 eV) for He and Ne respectively. Second, based on published very young volcanic zircon of ZHe ages coherent with other methods (e.g. [START_REF] Farley | U-Th)/He Dating of Phosphates: Apatite, Monazite, and Xenotime[END_REF], we can infer that for low alpha dose (~1-2x10 15 /g), He is retentive to surface temperature implying a Tc of ~60-100°C (Fig. 10B). Between low damage and the threshold, we could anticipate that Tc may increase following a non-linear (i.e. cubic or exponential) law as is the case for apatite [START_REF] Flowers | Apatite (U-Th)/He thermochronology using a radiation damage accumulation and annealing model[END_REF], as Ea increases until the EaT threshold (Fig. 10B). As damage fraction increases, clustering increases also and damaged domains with different sizes will form. In this case, the trapping energy can increase as proposed by [START_REF] Recanati | Helium trapping in apatite damage: Insights from (U-Th-Sm)/He dating of different granitoid lithologies[END_REF] for apatite. It is interesting to note that the highest closure temperatures obtained by [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF] can be explained by increasing trapping energy when damage dose increases up to a maximum (Fig. 8A). For Ne, the maximum activation energy predicted in this study is similarly very close to the one reported by van Soest et al. ( 2013), and we can anticipate a complex relationship of Ne diffusion with damage dose similar to He.

After some damage threshold, He diffusion increases (Tc decreases), although the Pyrenees samples [START_REF] Bosch | Timing of Eocene-Miocene thrust activity in the Western Axial Zone and Chaînons Béarnais (west-central Pyrenees) revealed by multi-method thermochronology[END_REF][START_REF] Vacherat | Rift-to-collision transition recorded by tectono-thermal evolution of the northern Pyrenees[END_REF] and work of [START_REF] Johnson | Inverted" zircon and apatite (U-Th)/He dates from the Front Range, Colorado: High-damage zircon as a lowtemperature (<50°C) thermochronometer[END_REF] indicate that the threshold dose may be lower than indicated by [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF]. The source of this discrepancy may be a complex function of damage accumulation and annealing, and more work on this issue is essential. For example, Ginster et al. (2018) document a complicated relationship between damage dose and annealing rate, implying a very complex system.

The other implication for a maximum of retentivity located in the vicinity of the nochain alpha damage percolation level (Fig. 9) would be that the dependence of closure temperature on damage fraction increase can be significantly steeper than modeled with Eq. 12 (Fig. 10B) if a single EaT value is assumed. This additional level of complexity will also have implications for developing very precise future quantitative models to predict individual TC values for single zircon grains. Such a potential scenario supported by natural data will require new experimental diffusion data to document carefully this key window of alpha-dose in the range of the no-chain percolation level.

Conclusions

In this contribution, we have used a series of computational and theoretical tools to provide a mechanistic view of the manner in which He and Ne diffuse in zircon. In particular, we have characterized the undamaged crystalline state and its energy structure at the atomic scale with respect to the diffusing atoms, and how trapping and obstruction as may result from damage or inherent defects affect retention at the grain scale. The set of insights gained can Our results corroborate previous theoretical and analytical work implying that pristine zircon is non-retentive for He due to low migration energies along the c-axis channels in the crystal structure. We also obtain the unexpected result that Ne moves even more freely than He through an unblocked c-axis channel, due to the distortion it causes in the insertion site. Both systems therefore require additional factors to explain their natural behavior. Using the Kinetic Monte Carlo method we show that obstruction of c-axis pathways forces the migration energetics to become more similar to the slower a/b direction as blockages increase, causing diffusion to become effectively isotropic on the grain scale and characterized by rectangularpath-jump energetics. This effect is insufficient to explain slow diffusion of He in zircon, but has a larger impact on Ne, as rectangular a/b transitions are much higher-energy, and rhombic ones are impassable. Of the mechanisms we inspect, we find that trapping is required to explain the high retention of He and Ne observed in natural zircon.

We consider the topology of radiation damage in the context of a data set from Pyrenees granites that indicates a transition from falling to rising diffusivity at a much smaller dose than posed in the [START_REF] Guenthner | Helium diffusion in natural zircon: radiation damage, anisotropy, and the interpretation of zircon (U-Th)/He thermochronology[END_REF] diffusivity model, closer to ~2-5×10 17 than 2×10 18 /g.

The transition dose may correspond to the no-chain alpha recoil percolation threshold of [START_REF] Ketcham | Geometric analysis of radiation damage connectivity in zircon, and its implications for helium diffusion[END_REF], suggesting that tip-to-tip recoil track connectivity is not sufficient to form fast diffusion pathways, but more interpenetrative damage intersections can. Increasing track intersections may also create more and deeper traps, and the evolution of not just trap abundance but also trap energetics may be necessary to arrive at a quantitative understanding of rare gas diffusion in natural zircon.

Tables captions

Table 1: Helium and neon insertion energies in zircon, hematite and apatite crystal structure.

Table 2: Migration energies in eV for both helium and neon in zircon compared with the theoretical studies of He diffusion [START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF][START_REF] Reich | Low-temperature anisotropic diffusion of helium in zircon: Implications for zircon (U-Th)/He thermochronometry[END_REF]Saadoune et al., 2009). Rect. and Rhom. are for the rectangular and rhombic gap. from this study are used for illustration. These results are used in a KMC code to calculate in 3D diffusion paths that accounts for the effects of obstruction and trapping. The example of diffusion in apatite is shown [START_REF] Djimbi | Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view[END_REF]. At the macroscopic level, for anisotropic diffusion and to account for crystal shape, a finite geometry Monte Carlo code [START_REF] Gautheron | A Monte Carlo approach of diffusion applied to noble gas/helium thermochronology[END_REF]) is used to calculate the (U-Th)/He age for a determined thermal history and effective closure temperature. The example of apatite and zircon crystals are shown. Finally, all these data can serve as the basis for a damage model. Pyrenean massifs [START_REF] Bosch | Timing of Eocene-Miocene thrust activity in the Western Axial Zone and Chaînons Béarnais (west-central Pyrenees) revealed by multi-method thermochronology[END_REF][START_REF] Vacherat | Rift-to-collision transition recorded by tectono-thermal evolution of the northern Pyrenees[END_REF]. The values for the 6-8-chain, nochain and fission-track percolation thresholds of [START_REF] Ketcham | Geometric analysis of radiation damage connectivity in zircon, and its implications for helium diffusion[END_REF] Emig is the migration energy, that can be for He 0.24, 0.66 eV for regular lattice (24 and 66 kJ/mol respectively). In a damage zone, He insertion can be lower and the resulting trapping energy will increase and can be up to 1.46 eV, resulting in a migration energy of 1.7 eV [START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF][START_REF] Reich | Low-temperature anisotropic diffusion of helium in zircon: Implications for zircon (U-Th)/He thermochronometry[END_REF]Saadoune et al., 2009). Rect.

and Rhom. are for the rectangular and rhombic gaps.

Rect.

Along c-axis diffusion between S1 and S1' sites [START_REF] Farley | U-Th)/He Dating of Phosphates: Apatite, Monazite, and Xenotime[END_REF] Ini.

Figure10

Electronic Annex Click here to download Electronic Annex: Appendix_sept2019.docx

  ) triangular dodecahedra parallel to c-axis(Figure 2A and 2B). The lateral connection of these chains provides two types of possible sites (S1 and S1') into which He or Ne atoms can be inserted; these insertion site positions, and the diffusion pathways between them, are identical He and Ne. The connection of these interstitial sites provides open channels along each direction in which He or Ne could diffuse(Figure 2A' and 2B'). In this part of our analysis, a periodic-DFT is run first to build an optimized model of zircon. Next, the model is used to identify diffusion channels and all possible interstitial insertion sites for both He and Ne atoms. As mentioned above, a 2×2×2 supercell was optimized as the host lattice (Figure2), which was large enough to simulate and study the He and Ne insertion in the zircon crystal lattice free of boundary effects. The c-axis channel is formed by alternating edge-sharing SiO4 tetrahedron and ZrO8 dodecahedron units (Figure 2A), and Figure 2A' illustrates the different He (or Ne) insertion sites along this pathway. In each c-channel four sites are inventoried for each lattice cell and they are energetically identical due to the symmetries of the crystal. They are spaced by 1.5 Å. Two types of sites are recognized: S1, which allows a jump to an adjacent c channel along a (jump along b forbidden), and S1', which allows a jump along b to an adjacent c channel. Along c, the sites form a sequence S1-S1'-S1-S1' … (Figure 2A'). Therefore, along a, diffusion pathways connect S1 sites. The same applies along b with the S1' sites. The channels along the a-and b-axis are energetically identical and they will be noted a/b-(axis or channel) for simplicity. To move along the a/b-axis, He (or Ne) atoms have to pass through a rhombic (Rhom.) or rectangular (Rect.) gap as shown on Figure 2B'. The former is formed by the intersection of Zr -Zr and O -O diagonals. The latter is formed by the intersection of Si -Si and O -O diagonals, and its area is much larger than the rhombic gap. Each channel consists of an alternating series of the two gap types (Figure 2B'); for example, along a the sequence is S1-Rhom.-S1-Rect.-S1-Rhom.-S1-Rect. This configuration makes the a/b-channels considerably less diffusive than the c-axis due to the narrow rhombic gaps. After channel determination, the two possible interstitial sites, called S1 and S1', were optimized. They are structurally symmetrical and energetically equivalent. Only one interstitial site environment exists for He and Ne, an octahedral-like site formed by the six oxygen atoms.

  Figure 3A and 3A' present the migration energies path results along the c-channel for He and Ne migration, respectively. The migration energies, reported in Table 2, are 0.24 eV for He (Figure 3A) and 0.22 eV (Figure 3A') for the Ne atom. According to a (b)-channel structures, the S1 (S1') sites are separated on one side by a rhombic gap and on the other by a rectangular one. Hence, two NEB calculations were performed for He and Ne migration, corresponding to the energy necessary to pass through each gap.

  . It is interesting to point out that only[START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF] recognizes the 1.5 Å inter-site distance along the c-axis, which results from the crystal symmetries.[START_REF] Reich | Low-temperature anisotropic diffusion of helium in zircon: Implications for zircon (U-Th)/He thermochronometry[END_REF] andSaadoune et al. (2009) missed this feature although a dip is visible in the energy landscape for the former, indicating that the minimal energy path was not found along the c-axis. Concerning the migration energies, which have the strongest impact on the diffusivity, we are in good agreement withSaadoune et al. (2009) along the c-axis, and at variance with the other calculations. Although we found the same energy pattern along this channel as[START_REF] Bengtson | He diffusion and closure temperatures in apatite and zircon: A density functional theory investigation[END_REF], our migration energy is lower. As the variations of the energy along this channel are small the computations of the barrier are very sensitive to the finding of the minimum energy path. When comparing to the previous calculations along the a/b directions we find the same barrier across the rhombic gaps, but a higher one through the rectangular gap (0.66 instead of ~0.45 eV, Figure4, Table

Figure 7

 7 Figure 7 presents the dependence of the diffusion coefficient Da,b and Dc (respectively

Figure 5B ,

 5B Figure 5B, lending support to the idea that our calculations have reproduced the natural range

  help leverage ongoing analytical efforts to unravel the complex behavior of the (U-Th)/He and (U-Th)/Ne systems in natural settings.Using Density Functional Theory and Transition State theory, we have calculated the insertion energies for He and Ne in ideal zircon, as well as the migration pathways along the a/b and c axes and the migration energies for each jump. For diffusion in the a/b direction, we document the energetics of the rectangular and rhombic pathways to neighboring interstitial sites. Our insertion energies for He are similar to previous work, but our migration energies and their interrelationships are different; notably, our work finds that the relative migration energies along the c-axis channel and through the rectangular and rhombic a/b gaps are approximately 0.24 vs. 0.66 and 2.60 eV for He, and 0.22 vs. 1.53 eV and infinite for Ne. One important finding is that our insertion energies, in combination with our migration energies, predict the maximum activation energies possible for the zircon He and Ne systems, which are consistent with literature values.

Table 4 :

 4 Diffusion coefficients (D0 and Ea) along a/b and c axes for different percentages of obstruction, as shown in Figure 8.

Figure captions Figure 1 :

 captions1 Figure captions

Figure 2 :

 2 Figure 2: Representation of the zircon crystal supercell (2×2×2) showing the channels

Figure 3 :

 3 Figure 3: Evolution of the minimal energy path (MEP) between insertion site S1 to S1'

Figure 4 :

 4 Figure 4: Schematic illustration of the different diffusion pathways in zircon identified

Figure 5 :

 5 Figure 5: Correlation of the activation energy Ea with frequency factor D0 (A) and

Figure 6 :

 6 Figure 6: Schematic representation of mechanisms that may affect diffusivity. (A) He

Figure 7 :

 7 Figure 7: Evolution of the diffusion coefficient as a function of inverse temperature

Figure 8 :

 8 Figure 8: Models of helium and neon closure temperature dependence on alpha dose

Figure 9 :

 9 Figure 9: ZHe age evolution as a function of the equivalent alpha-dose, for two

  are reported for comparison. Apatite fission track (AFT) and apatite (U-Th)/He (AHe) ages range from Bosch et al. (2016) and Vacherat et al. (2016) are also shown for comparison.

Figure 10 :

 10 Figure 10: (A) Schematic representation of the energetic levels of a diffusion pathway

  Figure1

Figure2S

  Figure2

  Figure4

Table 3 :

 3 Diffusion parameters for He and Ne in ideal lattice.

Table 1 . Insertion energies for He and Ne into interstitial sites in this and other studies. Insertion energy Zircon, this work Zircon, Saadoune et al. (2009) Hematite, Balout et al. (2017a, b) Apatite, Djimbi et al. (2015)

 1

					1.15 (a,b-plan);
	He (eV)	1.46	1.45	1.51
					0.66 (c-axis)
	Ne (eV)	3.10	-	3.38	-

Table

Table 2 . Migration energies (in eV) for He and Ne in zircon.

 2

			Helium			Neon
			a and b-channel		a and b-channel
		c-channel			c-channel
			Rect. Rhom.		Rect. Rhom.
	This work	0.24	0.66	2.60	0.22	1.53	none
	Reich et al. (2007)	0.14	0.46	1.05	-	-	-
	Saadoune et al. (2009)	0.22	0.44	2.68	-	-	-
	Bengtson et al. (2012)	0.44	0.43	2.64	-	-	-
	Migration energies in eV for both helium and neon in zircon compared with the theoretical
	studies of He diffusion

Table 3 . Diffusion parameters for He and Ne in ideal lattice. Saddle point E mig (eV) E a (kJ/mol) ν (THz) a (Å)

 3

		Rect.	0.66	63.7	5.50	3.35
	a and b-channel
	He	Rhom.	2.6	250.9	3.44	3.35
	c-channel		0.24	23.2	7.06	1.5
		Rect.	1.53	147.6	3.90	3.35
	a and b-channel
	Ne	Rhom.	none	none	none	3.35
	c-channel		0.22	21.2	1.36	1.5

Table 4 : Diffusion coefficients

 4

Acknowledgment

This project and salary for Duval Mbongo-Djimbi were funded by the Agence National de la Recherchegrant no. ANR-12-BS06-0005-01. DIM OXYMORE funded by the region île de France is thanked for salary for Hilal Balout. We thank Christophe Diarra for his precious help in managing the IPNO cluster GRIF (http:// www.grif.fr). William Guenthner is warmly thanked for sharing the diffusion data of his 2013 contribution. Daniele Cherniak, Peter Reiners and Martin Danišík are warmly thanked for their constructive reviews.

Appendix A and B -Supplementary data

Supplementary data associated with this article can be found, in the online version.


Séverine Le Roy

Fulbert Baudoin

Christian Laurent

G Teyssedre

G Teyssedre Introduction

Ionic Contri- Butions To

Introduction of Ionic Contributions to a Charge Transport Model for Dielectrics

come   L'archive ouverte pluridisciplinaire

Introduction of Ionic Contributions to a Charge

Transport Model for Dielectrics S. Le Roy 1 , F. Baudoin 1 , C. Laurent 1 , G. Teyssèdre 1

I. INTRODUCTION

Polymers are increasingly used as electrical insulation in systems of the electrical engineering, because of their properties (thermal, electrical, mechanical), their cost, and their possibility to be recycled for some of them. It is for example the case for high voltage direct current (HVDC) links, where cross-linked polyethylene (XLPE) has become a challenging material compared to oil impregnated paper. However, these materials can store charges when they are submitted to an electrical stress, increasing the electric field locally, and potentially leading to the material ageing and at last to its dielectric breakdown [START_REF] Montanari | Electrical degradation threshold of polyethylene investigated by space charge and conduction current measurements[END_REF]. The generated charges can originate from injection processes at the electrodes or from dissociation and ionization of species in the volume of the material. Dissociation processes can give rise to electronic species that can travel across the material and be extracted at the electrodes. It can also be ions, that can also travel if they are not too big, and accumulate at the opposite electrode. This charge accumulation can locally increase the electric field. However, it is difficult experimentally to differentiate electronic and ionic species, as space charge measurements only give access to a global positive or negative charge. Then, it is even more difficult to decipher the role of ions in the global space charge behavior.

In this contribution, we propose to simulate the impact of the presence of ions using a bipolar charge transport model, adding elementary generation and transport processes related to ions and combining it to electronic phenomena. In a first approach, only ions are simulated, and the impact of ion concentration and mobility are observed. Simulations are then performed for a more realistic case, where electronic and ionic species coexist, giving rise to already observed experimental pattern from space charge measurements. The objective of this paper is however not to reproduce experimental data, but to show that ions are important and should be taken into account in the further development of charge transport models.

II. MODEL DESCRIPTION

The one-dimension bipolar charge transport model (BCTM) used has already been described in the literature [START_REF] Roy | Description of charge transport in polyethylene using a fluid model with a constant mobility: fitting model and experiments[END_REF], and features electronic charge injection at each electrode, conduction between shallow traps via a constant mobility, trapping into deep traps, detrapping and recombination between electronic charges of opposite polarity. Ionic species, positive (A + ) and negative (B -), have also been taken into account. They are thought to be present inside the dielectric prior to any voltage application, arising from an initially neutral molecule AB. Hence, the quantity of molecules AB, and ions A + and B -is driven by the following equation: [START_REF] Montanari | Electrical degradation threshold of polyethylene investigated by space charge and conduction current measurements[END_REF] where k d is the dissociation coefficient, while k r is the recombination coefficient. It is to note that k r and k d do not have the same dimension, and are not equal. In a first simplified hypothesis, we will consider that a thermodynamic equilibrium has been reached prior to any voltage application. Hence, there is the same quantity of ions A + and B - everywhere inside the material. We will also consider that the dissociation and recombination rates are low compared to charge transport, so that no additional ions or additional neutral molecules AB are created during the experiment. This is a strong hypothesis, as it is certainly not the case in reality, but i/ the goal is here to understand how the presence of ions can influence the space charge pattern or the current profiles, and ii/ although a literature exists on dissociation and recombination coefficients in plasma science or even for liquids, it is not available for solid organic dielectrics. Moreover, no attempt has been made in this paper to link the molecule AB or ions A + and B -to real cases. Hence, A + and B - are present with the same initial amount C 0 , i.e. the net charge density remains zero prior to any voltage application. Ions can also move, and their mobility is thought as being constant. In a first attempt, we do not consider that these ions have a mobility dependence upon temperature, concentrations of sites/vacancies where ions can move or even electric field, even if it is the case in reality [START_REF] Kao | Electrical Transport in Solids[END_REF][START_REF] Dissado | Electrical Degradation and Breakdown in Dielectrics[END_REF]. There is no recombination between ions, i.e. no creation of additional molecules AB, and no recombination between an ion and an electronic charge.

The space and time equations describing the space charge behavior are the following, for each electronic and ionic carrier:

  ) t , x ( s x t , x j t ) t , x ( n i i i       (2) ε t) ρ(x, x t) E(x,    (3) x t) (x, n t) (x, D t) t).E(x, (x, t).μ (x, n t) (x, j i i i i i     (4)
Here, n i refers to the carrier density of specie i (C/m 3 ), j i refers to the current density (A/m 2 ), µ i the mobility, D i the diffusion coefficient, and s i refers to the source term for each specie i, and corresponds to all physical processes not linked to transport. For electrons and holes, it refers to trapping, detrapping and recombination. For ionic species, s i is equal to zero. t is the time, x is the space coordinate, ε the dielectric permittivity and E is the electric field.  is the net charge density, and is, in this case:

t) (x, n t) (x, n t) (x, n t) (x, n t) (x, n t) (x, n t) ρ(x, B et eµ A ht hµ         (5)
where n hµ , n ht , n eµ and n et are the mobile and trapped electronic charge density for holes and electrons, and n A + and n B -are the positive and negative ion charge densities. Electronic (i.e. electrons or holes) charge generation at each electrode is function of the electric field at this electrode and follows a modified Schottky law of the type [START_REF] Roy | Modelling space charge in a cable geometry[END_REF]:

                    1 4π t) , eE(x T k e exp T k ew A.T²exp t) , (x j 0/L b b inj 0/L inj  (6)
where A is the Richardson constant, T is the temperature, e the elementary charge, w inj the barrier height, k b refers to the Boltzmann constant, x 0/L is the position of the electrode (0 or L, thickness of the sample). The description of the model has already been published in [START_REF] Roy | Description of charge transport in polyethylene using a fluid model with a constant mobility: fitting model and experiments[END_REF]. It is a Finite Volume explicit model, using a third order accuracy in space and fourth order in time. The grid to discretize the insulating material is divided into 301 elements, non-uniform, having a smaller size when approaching the electrodes. A flux limiter (ULTIMATE) is used to correct the charge flux within the dielectric. The time steps follow the CFL condition. Parameters used for the simulations are presented in Table 1. All simulations are performed on a polyethylene based material, of relative permittivity 2.3, of thickness 200 µm. The experimental protocol consists of the application of a voltage of 5 kV (field of 25 kV/mm) during 3 hours, at T=25°C.

III. SIMULATION RESULTS FOR IONS TRANSPORT

Simulations have been performed when only ions are present, in order to see the impact of the presence of ions on the charge density, current and electric field. The initial ion concentration is set to 2 C/m 3 . It is difficult to anticipate the exact quantity of ion available for dissociation and/or transport. Taking cross-linked polyethylene (XLPE) as an example, when 1.3% of dicumylperoxide is used for crosslinking, the amount of residue, as cumyl alcohol (M=136.9 g/mol), in the material is around 0.5% in weight [START_REF] Guffond | Characterization and modeling of microstructure evolution of cable insulation system under high continuous electric field[END_REF]. The charge density associated with the total ionization of this molecule would be of the order of 10 6 C/m 3 . The initial ion density (2 C/m 3 ) chosen here is then really low, and has been set considering that a low amount of molecules can ionize.

A. Impact of ion mobility

Only A + is mobile, B -is considered as too big to move. The ion mobility has been changed from 10 -14 to 10 -17 m 2 /V/s. It is to note that the color scale for the space charge pattern is the same for all space charge patterns, and is given for charge densities from -2 to 2 C/m 3 (Fig. 1). The simulated space charge patterns are given in Fig. 2 for the four mobility values. All the simulations are presented with the influence charges on each electrode, linked to capacitive and to image charges. A Gaussian filter of width 10 µm was introduced to approach patterns usually obtained experimentally. For most of the space charge profiles, even if only the positive ions are mobile, only negative charges are observed. The specific behavior shows that positive ions move towards the ground electrode, leaving a net negative charge next to the top electrode first. This net negative charge, due to immobile negative ions, appears nearer to the ground electrode with time. Positive ions accumulate next to the ground electrode, but are not observed here due to the low quantity of initial ions, and the fact that influence charges at the electrode mask these positive charges. The higher the positive ions mobility, the sooner the material bulk becomes negative. The electric field after 3 hours of voltage application is 3.5 10 7 V/m at the bottom electrode, i.e. 40% higher than the applied electric field, while the electric field at the top electrode is around 1.5 10 7 V/m for a positive ion mobility of 10 -15 m 2 /V/s. This variation of the electric field at the electrodes will have a strong impact on the electronic charge generation at the electrodes. The range of mobility values explored is close to the one obtained for electronic charges when we could expect lower mobilities for ions. Increasing the mobility make it possible to observe some specific patterns within 3 hours stressing. Experimentally, such kind of pattern is observed for XLPE (see as example [START_REF] Roy | Ion generation and transport in low density polyethylene under electric stress[END_REF]) soon after voltage application. It could mean that i/ the ion mobility is effectively rather high, ii/ the charge that moves after dissociation is not an ion but an electronic charge, iii/ other processes are involved prior to voltage application, and the initial state of charge is not the one that has been taken for these calculations (e.g. mix of electronic and ionic carriers).

Fig. 3 presents the current density as a function of time for the different positive ion mobility values. When ions are present, there is a specific behavior in the current density, i.e. the presence of a 'plateau'. The higher the ion mobility is the higher the current density on this plateau. When the ion mobility is increased, the drop in the current density appears sooner in time, when most of the mobile ions have crossed the sample and accumulate at the counter electrode. This is also a kind of result observed experimentally in XLPE [START_REF] Cènes | Implementation of space charge measurement using the Pulsed Electro-Acoustic method during ageing of HVDC model cable[END_REF].

B. Impact of ion concentration

For these simulations, the mobility of positive ions is set to 10 -15 m 2 /V/s, while the negative ions are still immobile. The initial ion concentration has been varied from 0.5 to 10 C/m 3 . Fig. 4 shows the simulated space charge profiles obtained for different initial ion densities. The space charge pattern does not change to a large extent when the initial ion density is increased. The dynamic is almost the same for all initial ion densities, the net charge varying slightly. Only when the initial ion density is high (5 C/m 3 ), positive charge accumulation is observed at the counter electrode. As for the previous case, positive charges are present at the vicinity of the cathode, but are masked by the negative charges and the influence charges. When the initial ion density is high, the negative charge density is not enough to mask the presence of positive charges at the cathode. An increase by one decade of the initial ion concentration increases the value of the current density on the plateau by one decade, but does not change the transient current shape (not shown here). Looking at the electric field value at the cathode (Fig. 5), an increase of the initial ion concentration has a large impact on the electric field value. An increase of a decade of the ion density implies a value of electric field which is multiplied by two. This then has a direct impact on the electronic charge injection. Simulations have also been performed for mobile positive and negative ions. The space charge pattern and current density as a function of time for these cases do not change to a large extent as long as negative ions have a lower mobility compared to positive ones. If negative ions have a higher mobility than positive ones, then the charge pattern would be reversed, i.e. a dominant positive charge that seems to arise from the cathode.

IV. COMBINED IONIC AND ELECTRONIC TRANSPORT

For this case study, electronic and ionic species have been considered. The parameters for electronic species are the ones that have been optimized in [START_REF] Roy | Description of charge transport in polyethylene using a fluid model with a constant mobility: fitting model and experiments[END_REF], for a constant effective mobility for electrons and holes, see Table 1. For ionic species the initial ion density is 2 C/m 3 and the positive ion mobility is 10 -15 m 2 /V/s. Fig. 6 presents the net charge density as a function of time and space for this case study. A positive charge front is observed arising from the anode. These charges are injected holes that cross the dielectric. The positive charge front disappears slowly in the middle of the sample, progressively replaced by negative charges. These negative charges are the negative ions that cannot move and are 'uncovered' progressively by the positive ions leaving their initial position to accumulate at the ground electrode (cathode). These negative ions completely mask the positive charge injection and continuous flow towards the cathode. It is to note that the injection flux does not change to a large extent at the anode, as the net charge density is small. At the cathode however, the positive ion accumulation in this region implies a higher negative charge injection after 3000 s, i.e. when most of the positive ions have migrated towards the cathode. The simulated current density transient associated to this case study is given in Fig. 7. For comparison, the current density when only electronic carriers are considered is also presented. The current density shape is different without and with ions, as a small increase is observed below 10 s, and then the current density is quasi stable during 1000 s. This behavior is directly linked to the movement of positive ions. Once most of them accumulate close to the cathode, the current density decreases, to a value higher than the one calculated without ions. Although ionic and electronic carriers are not considered to interact, the total current is not just the superposition of the two contributions taken independently.

V. CONCLUSIONS

A bipolar charge transport model has been implemented to take into account the presence of positive and negative ions and their transport inside an insulating material. A small amount of ions has a large impact on the space charge density, on the electric field at the electrode, and on the current density. The simulated space charge and current transients has clear common features with what is measured in reality for insulating polymers which are known or supposed to contain ionic contributions to conduction. More work needs now to be done to take into account the neutral molecule impact and the recombination between positive and negative ions as first step.

Fig. 1 .

 1 Fig. 1. Colorbar used for all space charge patterns.

Fig. 3 .

 3 Fig. 3. Simulated current density vs. time for different positive ion mobilities.

2 .

 2 Simulated space charge patterns calculated when only ions are present, and when only positive ions are mobile, for different mobilities: a) 10 -14 m 2 /V/s, b) 10 -15 m 2 /V/s, c) 10 -16 m 2 /V/s, d) 10 -17 m 2 /V/s.

Fig. 5 .

 5 Fig. 5. Simulated electric field at the cathode as a function of time for different initial ion concentrations.

Fig. 6 .

 6 Fig. 6. Simulated space charge pattern calculated when electronic and ionic carriers are considered, for an initial ion density of 2 C/m 3 . Applied voltage 5 kV, temperature 25°C.

Fig. 7 .

 7 Fig. 7. Simulated current density function of time calculated for electronic carriers only and when electronic and ionic carriers are considered, for an initial ion density of 2 C/m 3 . Applied voltage 5 kV, temperature 25°C.

4 .

 4 Simulated space charge patterns calculated when only ions are present, and for different ion concentrations: a) 0.5 C/m 3 , b) 1 C/m 3 , c) 2 C/m 3 , d) 5 C/m 3 .

TABLE I .

 I PARAMETERS USED IN THE SIMULATIONS

	Symbol	value	units
	trapping coefficients
	Be electrons	1. 10 -1	s -1
	Bh holes	2. 10 -1	s -1
	Maximal trap densities
	noet for electrons	100	C.m -3
	noht for holes	100	C.m -3
	injection barrier heights
	wei for electrons	1.27	eV
	whi for holes	1.16	eV
	Detrapping barrier heights
	wtre for electrons	0.96	eV
	wtrh for holes	0.99	eV
	Mobility
	µe for electrons	10 -14	m 2 /V/s
	µh for holes	2 10 -13	m 2 /V/s
	µA for positive ions	10 -17 to 10 -14	m 2 /V/s
	µB for negative ions	0	m 2 /V/s
	Initial ions density	0.5 to 10	C/m 3