Jehiel Nteme
email: jehiel.nteme-mukonzo@univ-orleans.fr

Stéphane Scaillet

Pascal Brault

Laurent Tassan-Got

Atomistic Simulations of 40 Ar Diffusion in Muscovite

Keywords: 40 Ar diffusion, Muscovite, Molecular Dynamics, Transition State Theory, Thermochronology

Muscovite ranks among the most commonly dated minerals in 40 Ar/ 39 Ar geochronology. Yet, its use in thermochronological reconstructions is hampered by the lack of reliable data on its 40 Ar diffusional behavior. In this contribution, we investigate 40 Ar lattice diffusion in muscovite at the atomic scale using Molecular Dynamics (MD) simulations combined with Nudged Elastic Band (NEB) and Transition State Theory (TST). Classical MD simulations of 40 Ar recoil dynamics in 2M 1 muscovite reveal that 40 Ar initially resides predominantly in the interlayer region, close to its production site. Systematic computations of migration barriers coupling NEB with TST identify the divacancy mechanism as the more energetically favorable pathway for 40 Ar diffusion in the interlayer region, with characteristic enthalpy of motion E = 66 kcal.mol -1 and vibrational entropy term D 0 = 6 × 10 -4 cm 2 .s -1 . For typical cooling rates between 1-100 • C.Ma -1 and grain size varying from 0.1 and 1 mm, these parameters predict closure temperatures significantly higher (∼ 200 • C) than currently accepted maximum estimates (∼ 500 • C). Consistent with long-standing empirical evidence, our theoretical results downplay the role of purely thermally activated diffusion in promoting efficient 40 Ar transport in ideal (stoichiometrically stable and undefective) muscovite. Along with experimental and field-based evidence, they call for more complex physics to explain the 40 Ar retention properties of natural muscovite, most notably by considering crystal-chemical disequilibrium interactions and the reactivity of the interlayer with the external medium.

Introduction

Diffusion is the main mechanism governing atomic transport in minerals. It is fundamental to the understanding of many geological processes in terms of kinetics of metamorphic transformations and closure behavior of radio-isotopic systems [START_REF] Dodson | Closure temperature in cooling geochronological and petrological systems[END_REF]. Among these, the 40 Ar/ 39 Ar dating technique has been widely applied for studying Earth processes in the thermal range controlling the retention behavior of 40 Ar-bearing minerals [START_REF] Mcdougall | Geochronology and Thermochronology by the 40 Ar/ 39 Ar Method[END_REF]. Being non bonded, 40 Ar is expected to be much more mobile than silicateforming elements and therefore much more sensitive to thermal effects. The commonly observed discordance between for example 40 Ar/ 39 Ar and U/Pb ages has been attributed to the fact that common minerals lose their 40 Ar at crustal medium to moderate temperatures, giving apparent 40 Ar/ 39 Ar ages that do not reflect the time of crystallization but the time when the mineral was sufficiently cold to become effectively closed to 40 Ar loss by diffusion [START_REF] Hart | Mineral ages and metamorphism[END_REF][START_REF] Armstrong | K-Ar dating of plutonic and volcanic rocks in orogenic belts[END_REF]. Such a closed-system behavior has been theoretized in terms of a physically defined closure temperature (T c ) taking into account the Arrhenian dependence of intracrystalline (lattice) diffusion to temperature according to [START_REF] Dodson | Closure temperature in cooling geochronological and petrological systems[END_REF]:

T c = E R ln AR T 2 c D 0 /a 2 E dT /dt (1)
in which E is the activation energy and D 0 the pre-exponential factor. dT /dt is the cooling rate, a the diffusion radius, A a numerical constant appropriate for the chosen geometry for diffusion, and R the gas constant. This expression emphasizes the importance of the diffusion parameters E and D o in any quantitative interpretation of 40 Ar/ 39 Ar metamorphic ages in terms of thermal-kinetic crustal processes.

Among those minerals commonly dated by 40 Ar/ 39 Ar, muscovite stands apart due to its high potassium content and ubiquity in most crustal lithologies, all the way down from the epizone to subduction-zone metamorphic conditions [START_REF] Mcdougall | Geochronology and Thermochronology by the 40 Ar/ 39 Ar Method[END_REF].

However, the 40 Ar diffusion behavior in muscovite remains poorly understood. Early on, maximum retention temperatures of ca. 350 • C to as high as 400 • C (or greater) were inferred from field calibration of apparent ages vs. metamorphic grade [START_REF] Purdy | K-Ar ages on rock forming minerals from the Central Alps[END_REF] 2 and pioneering diffusion experiments [START_REF] Hart | Mineral ages and metamorphism[END_REF]. Other 40 Ar diffusivity estimates were later proposed by [START_REF] Hames | An empirical evaluation of the argon diffusion geometry in muscovite[END_REF], [START_REF] Lister | Modelling the effect of arbitrary PTt histories on argon diffusion in minerals using the MacArgon program for the Apple Macintosh[END_REF] and [START_REF] Kirschner | Staircase 40 Ar/ 39 Ar spectra of fine-grained white mica: Timing and duration of deformation and empirical constraints on argon diffusion[END_REF]. [START_REF] Hames | An empirical evaluation of the argon diffusion geometry in muscovite[END_REF] and [START_REF] Lister | Modelling the effect of arbitrary PTt histories on argon diffusion in minerals using the MacArgon program for the Apple Macintosh[END_REF] re-evaluated the early thermal-hydrothermal experimental data of [START_REF] Robbins | Radiogenic argon diffusion in muscovite under hydrothermal conditions[END_REF] to suggest the values E = 52 kcal.mol -1 and D 0 = 0.04 cm 2 .s -1 [START_REF] Hames | An empirical evaluation of the argon diffusion geometry in muscovite[END_REF], and E = 41.8 kcal.mol -1 and D 0 = 3.352 × 10 -7 cm 2 .s -1 [START_REF] Lister | Modelling the effect of arbitrary PTt histories on argon diffusion in minerals using the MacArgon program for the Apple Macintosh[END_REF]. [START_REF] Kirschner | Staircase 40 Ar/ 39 Ar spectra of fine-grained white mica: Timing and duration of deformation and empirical constraints on argon diffusion[END_REF] inferred an activation energy of E = 58 kcal.mol -1 with D 0 = 0.24 cm 2 .s -1 from an empirical thermal model binned to 40 Ar/ 39 Ar ages in muscovites grown during nappe thrusting in the Western Alps. Nowadays, commonly adopted values are those determined by [START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF] from hydrothermal experiments conducted in the range 730-600

• C and P = 10, 20 kbar (E = 63 kcal.mol -1 ; D 0 = 2.3 cm 2 .s -1 ). Some other estimates were inferred by [START_REF] Forster | 40 Ar/ 39 Ar geochronology and the diffusion of 39 Ar in phengite-muscovite intergrowths during step-heating experiments in vacuo[END_REF] and [START_REF] Lister | White mica 40 Ar/ 39 Ar age spectra and the timing of multiple episodes of high-P metamorphic mineral growth in the Cycladic eclogite-blueschist belt, Syros, Aegean Sea, Greece[END_REF] using in vacuo experiments (E = 67 kcal.mol -1 and D 0 /a 2 = 1.2 × 10 9 s -1 ).

These parameters imply theoretical 40 Ar-closure temperatures lower than 500 • C [START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF]). However, cases of effective 40 Ar retention have been consistently reported for muscovites exposed to temperatures in excess of this value [START_REF] Monié | Preservation of Hercynian 40 Ar/ 39 Ar ages through high-pressure low-temperature Alpine metamorphism in the Western Alps[END_REF][START_REF] Scaillet | 40 Ar/ 39 Ar laser-probe dating by step heating and spot fusion of phengites from the Dora Maira nappe of the Western Alps, Italy[END_REF][START_REF] Scaillet | Mg/Fe and [(Mg,Fe)Si-Al 2 ] compositional control on argon behaviour in high-pressure white micas: A 40 Ar/ 39 Ar continuous laser-probe study from the Dora-Maira nappe of the internal Western Alps, Italy[END_REF][START_REF] Monié | 40 Ar/ 39 Ar dating in coesite-bearing and associated units of the Dora Maira massif[END_REF][START_REF] Arnaud | Evidence for excess argon during high pressure metamorphism in the Dora Maira Massif (western Alps, Italy), using an ultra-violet laser ablation microprobe 40 Ar-39 Ar technique[END_REF][START_REF] Scaillet | Excess 40 Ar transport scale and mechanism in high-pressure phengites: A case study from an eclogitized metabasite of the Dora-Maira nappe, western Alps[END_REF][START_REF] Vincenzo | Evidence of a close link between petrology and isotope records: constraints from SEM, EMP, TEM and in situ 40 Ar-39 Ar laser analyses on multiple generations of white micas (Lanterman Range, Antarctica)[END_REF][START_REF] Rapaport | The art of molecular dynamics simulation[END_REF], 2006;[START_REF] Balogh | Argon and fission track dating of Alpine metamorphism and basement exhumation in the Sopron Mts.(Eastern Alps, Hungary): thermochronology or mineral growth?[END_REF][START_REF] Hames | Single-crystal 40 Ar/ 39 Ar age variation in muscovite of the Gassetts Schist and associated gneiss, Vermont Appalachians[END_REF][START_REF] Itaya | Regional-Scale Excess Ar wave in a Barrovian type metamorphic belt, eastern Tibetan Plateau[END_REF][START_REF] Allaz | The effects of retrograde reactions and of diffusion on 40 Ar-39 Ar ages of micas[END_REF][START_REF] Beltrando | Preservation of sub-microscopic structural relicts in micas from the Gran Paradiso Massif (Western Alps): Implications for 40 Ar-39 Ar geochronology[END_REF][START_REF] Halama | Effects of fluid-rock interaction on 40 Ar/ 39 Ar geochronology in high-pressure rocks (Sesia-Lanzo Zone, Western Alps)[END_REF][START_REF] Schertl | Tracking the incidence of excess argon in white mica Ar-Ar data from UHP conditions to upper crustal levels in the Dora-Maira Massif, Western Alps[END_REF][START_REF] Laurent | 40 Ar behaviour and exhumation dynamics in a subduction channel from multi-scale 40 Ar/ 39 Ar systematics in phengite[END_REF]. While most of these studies relate to high-pressure phengites from eclogitefacies rocks for which the combined effect of pressure and composition (e.g., Tschermak substitution) may have influenced the 40 Ar diffusion behavior [START_REF] Wijbrans | 40 Ar/ 39 Ar dating of white micas from an Alpine highpressure metamorphic belt on Naxos (Greece): the resetting of the argon isotopic system[END_REF][START_REF] Scaillet | Mg/Fe and [(Mg,Fe)Si-Al 2 ] compositional control on argon behaviour in high-pressure white micas: A 40 Ar/ 39 Ar continuous laser-probe study from the Dora-Maira nappe of the internal Western Alps, Italy[END_REF][START_REF] Scaillet | K-Ar ( 40 Ar/ 39 Ar) geochronology of ultrahigh pressure rocks[END_REF], a high retentivity has also been suggested in other settings for which compositional and pressure effects are not prevalent [START_REF] Vincenzo | The relationship between tectono-metamorphic evolution and argon isotope records in white mica: constraints from in situ 40 Ar-39 Ar laser analysis of the Variscan basement of Sardinia[END_REF][START_REF] Hames | Single-crystal 40 Ar/ 39 Ar age variation in muscovite of the Gassetts Schist and associated gneiss, Vermont Appalachians[END_REF][START_REF] Itaya | Regional-Scale Excess Ar wave in a Barrovian type metamorphic belt, eastern Tibetan Plateau[END_REF][START_REF] Allaz | The effects of retrograde reactions and of diffusion on 40 Ar-39 Ar ages of micas[END_REF]. The discordance between experimental predictions and geologic evidence raises the fundamental question as to the meaning of such 40 Ar/ 39 Ar ages: Do they reflect cooling or (re)crystallization? The answer requires a better understanding of the 40 Ar diffusion behavior in muscovite as a function of several parameters including temperature, chemical composition, and other intensive variables such as pressure.

temperature. Unlike previous studies that investigated the 40 Ar diffusivity in muscovite experimentally and empirically, we re-evaluate the 40 Ar diffusion dependence to temperature via atomistic simulations. These simulations, which include various methods (e.g., Density Functional Theory -DFT, Molecular Dynamics -MD, Kinetic Monte-Carlo -KMC, etc.), show many advantages over laboratory experiments. While the determination of rare-gas diffusivities in hydrothermal experiments is based on finite-state observables (fractional loss), atomistic simulations investigate the diffusion behavior at the atomic scale, allowing the identification and characterization of the diffusion mechanism (e.g., interstitial vs. vacancy mechanism) apt to explain the bulk migration of the diffusing species across the crystal. Also, atomistic simulations permit to analyze the diffusion process in a perfect lattice structure, unlike laboratory experiments that use natural crystals which commonly contain impurities and microstructural defects.

Atomistic simulations have been widely applied to the study of solid-state diffusion in the field of Material Science. Quite recently, they have been used for the theoretical characterization of low-temperature (fast-diffusing) thermochronometers like He and Ne in hematite, goethite, zircon, apatite and quartz via a multi-scale approach combining DFT with KMC [START_REF] Saadoune | He incorporation and diffusion pathways in pure and defective zircon ZrSiO 4 : A density functional theory study[END_REF]Bengston et al., 2012;[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][START_REF] Bassal | Neon diffusion in goethite, α-FeO(OH): a theoretical multi-scale study[END_REF][START_REF] Domingos | Simulations and Experiments Reveal Effect of Nanopores on Helium Diffusion in Quartz[END_REF][START_REF] Gautheron | A multi-method, multi-scale theoretical study of He and Ne diffusion in zircon[END_REF]. DFT simulations of rare-gas (He, Ne, Ar, Kr, Xe) diffusion in hydrated minerals have further been applied to quantify the efficiency of noble gas recycling in subduction zones [START_REF] Wang | Diffusion of noble gases in subduction zone hydrous minerals[END_REF]. DFT has also been successfully used to constrain He diffusion in olivine, quartz and coesite at mantle pressures [START_REF] Wang | Helium diffusion in olivine based on first principles calculations[END_REF][START_REF] Liu | A firstprinciples study of helium diffusion in quartz and coesite under high pressure up to 12GPa[END_REF]. In this paper, we investigate 40 Ar diffusion in muscovite by combining MD with TST. Our results reveal 40 Ar diffusion kinetics several orders of magnitude slower than current estimates.

Methods

Crystal structure model

The crystal structure of muscovite [KAl 2 (Si 3 Al)O 10 (OH) 2 ] consists of a stacking of tetrahedral -octahedral -tetrahedral (T-O-T) layers held together across the interlayer space by K + ions, which compensate the negative structural charge resulting from the substitution of Al for Si [START_REF] Bailey | Classification and structures of the Micas[END_REF]. Tetrahedral silicon layers and dioctahedral aluminum layers are connected through bridging oxygens. In this work, we use a 2M 1 muscovite atomic model based on the structural data of [START_REF] Catti | Thermal strain analysis in the crystal structure of muscovite at 700[END_REF]. The model consists of a periodic supercell of 6 × 3 × 1 unit cells along the crystallographic directions a, b, and c, respectively (Fig. 1a).

Tetrahedral Al atoms are distributed in an irregular pattern while obeying the short range Lowenstein rule [START_REF] Loewenstein | The distribution of aluminum in the tetrahedra of silicates and aluminates[END_REF] and the principle of local charge balance [START_REF] Herrero | Si, Al distribution in micas; analysis by high-resolution 29 Si NMR spectroscopy[END_REF][START_REF] Herrero | 29 Si MAS-NMR spectroscopy of mica-type silicates: Observed and predicted distribution of tetrahedral Al-Si[END_REF]. Accordingly, no IV Al-O-IV Al bond was allowed and each hexagonal tetrahedral ring contains either one or two Al atoms (Fig. 1b). Each K + is located between one Si 4 Al 2 ring and one Si 5 Al 1 ring to ensure an equitable charge distribution.

Computational details

Classical Molecular dynamics (MD) consists in simulating the trajectory of a system of interacting particles (atoms and/or molecules) over time by solving Newton's equations of motion (Frenkel and Smith, 2002;[START_REF] Rapaport | The art of molecular dynamics simulation[END_REF]. Interatomic interactions are computed based on force fields, which are parameterized using experimental results or ab initio simulations. We performed MD simulations using the simulation package LAMMPS [START_REF] Thompson | LAMMPS -a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales[END_REF]. The interatomic interactions were computed using ClayFF, a non-reactive force field that describes hydrated mineral systems through essentially nonbonded potentials [START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF]. Since its development, this force field has been widely used in MD simulations of phyllosilicates and other common minerals for which it has been shown to faithfully reproduce the structure and many thermodynamic properties [START_REF] Greathouse | Molecular dynamics simulation of uranyl (VI) adsorption equilibria onto an external montmorillonite surface[END_REF][START_REF] Greathouse | Water structure and aqueous uranyl (VI) adsorption equilibria onto external surfaces of beidellite, montmorillonite, and pyrophyllite: Results from[END_REF][START_REF] Kirkpatrick | Molecular dynamics modelling of hydrated mineral interlayers and surfaces: structure and dynamics[END_REF][START_REF] Kumar | Hydration, swelling, interlayer structure, and hydrogen bonding in organolayered double hydroxides: Insights from molecular dynamics simulation of citrate-intercalated hydrotalcite[END_REF][START_REF] Larentzos | An ab initio and classical molecular dynamics investigation of the structural and vibrational properties of talc and pyrophyllite[END_REF][START_REF] Liu | Interlayer structure and dynamics of alkylammonium-intercalated smectites with and without water: A molecular dynamics study[END_REF][START_REF] Liu | Effects of layer-charge distribution on the thermodynamic and microscopic properties of Cs-smectite[END_REF][START_REF] Kerisit | Molecular dynamics simulations of the orthoclase (001)-and (010)-water interfaces[END_REF][START_REF] Meleshyn | Two-dimensional ordering of water adsorbed on a mica surface at room temperature[END_REF][START_REF] Wander | Structural and dielectric properties of quartz-water interfaces[END_REF][START_REF] Cygan | Molecular models and simulations of layered materials[END_REF][START_REF] Malani | Influence of hydrophilic surface specificity on the structural properties of confined water[END_REF][START_REF] Ferrage | Hydration properties and interlayer organization of water and ions in synthetic Na-smectite with tetrahedral layer charge. Part 1. Results from X-ray diffraction profile modeling[END_REF][START_REF] Kalinichev | Molecular dynamics computer simulations of the effects of hydrogen bonding on the properties of layered double hydroxides intercalated with organic acids[END_REF][START_REF] Tao | Swelling of K + , Na + and Ca 2+ -montmorillonites and hydration of interlayer cations: a molecular dynamics simulation[END_REF]Argyris et al., 2011a,b;[START_REF] Ho | Interfacial water on crystalline silica: a comparative molecular dynamics simulation study[END_REF][START_REF] Teich-Mcgoldrick | Molecular dynamics simulations of structural and mechanical properties of muscovite: pressure and temperature effects[END_REF][START_REF] Churakov | Mobility of Na and Cs on montmorillonite surface under partially saturated conditions[END_REF][START_REF] Loganathan | On the hydrogen bonding structure at the aqueous interface of ammonium-substituted mica: A molecular dynamics simulation[END_REF][START_REF] Morrow | Structure, energetics, and dynamics of smectite clay interlayer hydration: Molecular dynamics and metadynamics investigation of Na-hectorite[END_REF]. The pairwise energy for atoms i and j separated by a distance r ij is given by

E ij = e 2 q i q j 4πϵ 0 r ij + 4ϵ ij σ ij r ij 12 - σ ij r ij 6 ( 2
)
where e is the fundamental charge of an electron, q i and q j are partial atomic charges and ϵ 0 is the permittivity of vacuum. The Lennard-Jones parameters σ ij and ϵ ij correspond to the zero-crossing distance for the energy and the depth of the potential energy well, respectively.

Interaction between hydroxyl atoms are computed with a bonded potential that describes the bond stretch energy by a simple harmonic term:

E bond ij = k 1 (r ij -r 0 ) 2 (3)
where k 1 is a force constant and r 0 the equilibrium bond length. Interaction parameters used in this study are provided in Table 1. As there are no parameters for Ar in ClayFF, we assigned it a zero charge and Lennard-Jones parameters as reported by [START_REF] White | Lennard-Jones as a model for argon and test of extended renormalization group calculations[END_REF].

The nonbonded parameters are given for interactions between identical atom types (i = j). Interaction parameters between unlike atoms (i ̸ = j) are calculated using the Lorentz-Berthelot mixing rules:

σ ij = 1 2 (σ i + σ j ) (4)
ϵ ij = √ ϵ i ϵ j (5)
The muscovite atomic model was optimized in a series of preliminary calculations consisting of energy minimization followed by 250 ps MD runs performed in the isobaric-isothermal (NPT) ensemble using the Nosé-Hoover thermostat and barostat, with damping parameters for temperature and pressure of 0.01 and 1 ps, respectively, and a time step of 0.001 ps.

Short-range interactions were evaluated with a cut-off distance of 10 Å whereas long-range interactions were treated using the Ewald summation. The 40 Ar diffusion energy barriers were quantified using the Nudged Elastic Band (NEB) method [START_REF] Jónsson | Nudged elastic band method for finding minimum energy paths of transitions[END_REF], which determines the minimum energy path (MEP) between two stable sites by solving simultaneously independent system images, each one containing one 40 Ar atom located on the path connecting the stable sites. Images are linked by a spring force to ensure the continuity of the path and the equal spacing between images along the MEP. The migration energy is obtained by taking the difference in energy between the saddle point and the initial stable site. In this study, eight images were generated and all the NEB calculations were performed with a spring constant of 23 kcal.mol -1 . Å-2 .

The transport rate based on those migration barriers was evaluated using the Transition State Theory (TST) that expresses the jump probability γ as:

γ = ν exp -E RT ( 6
)
where ν is the attempt frequency, E the migration barrier, R the gas constant, and T the temperature of the system (in Kelvin). The attempt frequency is given by:

ν = 3N i=1 ν init i 3N -1 i=1 ν * i (7)
in which ν init i corresponds to the modal vibrational frequencies at the stable site and ν * i represents the non-imaginary modal frequencies at the saddle point [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF]. For most solids, the attempt frequency is comprised between 10 12 s -1 and 10 13 s -1 . This parameter is usually approximated by a fixed value within this interval to avoid the computationally expensive determination of modal frequencies. In this work, we used ν = 5 × 10 12 s -1 , which is close to the attempt frequency of 4.61 × 10 12 s -1 computed via DFT for Ar in lizardite [START_REF] Wang | Diffusion of noble gases in subduction zone hydrous minerals[END_REF]. The diffusion coefficient is related to the transition rate via:

D = ga 2 γ = ga 2 ν exp -E RT ( 8
)
where g is a numerical constant depending on the geometry of the lattice and a is the distance between two neighboring stable sites. The frequency factor D 0 is given by:

D 0 = ga 2 ν (9)

Force field validation

To validate the application of the ClayFF force field to the study of muscovite, we compared the lattice parameters computed at ambient conditions (P = 0.001 kbar and T = 25 • C) for our muscovite atomic model with experimental data [START_REF] Catti | Thermal strain analysis in the crystal structure of muscovite at 700[END_REF]. The results (Table 2) show that ClayFF reproduces the experimental crystal structure well. The lattice parameters a, b, c and the angles α, β, γ are all within 2 % (and in most cases 0.5 %)

of the experimental values. The interlayer separation is also well reproduced (within 1 % of experimental determinations). Moreover, ClayFF also faithfully reproduces the stacking of adjacent T-O-T layers. In the interlayer region, the basal oxygens of one T-O-T layer are directly above those of the adjacent layer, in keeping with experimental determinations.

Each potassium atom is located at the center of the upper and lower tetrahedral rings. As expected for a di-octahedral mica, the hydroxyl groups point towards the vacant octahedral site and remain subparallel to the (001) plane.

40 Ar initial position from interlayer recoil dynamics

The investigation of 40 Ar diffusion at the atomic scale requires the determination of its insertion site within the muscovite structure. Since 40 Ar is produced with a kinetic energy of 660 kcal.mol -1 corresponding to a velocity of 0.11765 Å.fs -1 [START_REF] Dong | Mechanisms of argon retention in clays revealed by laser 40 Ar-39 Ar dating[END_REF][START_REF] Szczerba | Molecular modeling of the effects of 40 Ar recoil in illite particles on their K-Ar isotope dating[END_REF], the 40 Ar insertion site can be determined by modelling the natural 40 Ar recoil via MD simulations. Such simulations have been performed by [START_REF] Szczerba | Molecular modeling of the effects of 40 Ar recoil in illite particles on their K-Ar isotope dating[END_REF], who investigated the natural recoil of 40 Ar in illite, a clay mineral structurally similar to muscovite.

Using both REAXFF and ClayFF force fields, they performed several simulations in which In this work, we simulate the trajectory of the 40 Ar recoil in muscovite using the same methodology as [START_REF] Szczerba | Molecular modeling of the effects of 40 Ar recoil in illite particles on their K-Ar isotope dating[END_REF]. Since the penetration of 40 Ar into the T-O-T layer induces deformations in the crystal structure [START_REF] Szczerba | Molecular modeling of the effects of 40 Ar recoil in illite particles on their K-Ar isotope dating[END_REF], and ClayFF being a non-reactive force field designed to model undeformed structures [START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF], we investigated recoil trajectories parallel to the interlayer region only. We consider that 40 Ar penetrates into the muscovite T-O-T layer in the same proportions as in illite, due to the similarity of their structures. We define the recoil angle ϕ as the angle between the crystallographic direction a and the recoil vector. The simulations were carried out for recoil angles from 0 • to 360 • with a step of 20 • . In order to simulate a system under geological conditions, we equilibrated first the system in the NPT ensemble, at a temperature of 500 • C and a pressure of 5 kbar. The recoil simulations were then carried out in the NVT ensemble, for a duration of 250 ps with a time step of 1 fs.

The simulations reveal that, in the absence of vacancies in the immediate vicinity of the parent 40 K initial site, when the recoil vector is oriented towards a neighboring potassium atom (i.e. for ϕ = 0 • , 60 • , 120 • , 180 • and 240 • ) the radiogenic 40 Ar bounces back from that neighboring K + and returns to its initial interlayer position. For the other recoil angles, 40 Ar is first propelled into an interstitial position where it remains for about 50 ps, before either returning in its initial position or replacing a neighboring K + which is in turn pushed into the 40 Ar initial site. Along with [START_REF] Szczerba | Molecular modeling of the effects of 40 Ar recoil in illite particles on their K-Ar isotope dating[END_REF], this is the first theoretical verification that radiogenic 40 Ar initially resides predominantly close to its production site in K-bearing phyllosilicates [START_REF] Dong | Mechanisms of argon retention in clays revealed by laser 40 Ar-39 Ar dating[END_REF], and not many unit cells away as sometimes suggested.

The variation of the total energy of the system when the radiogenic 40 Ar is propelled into the interstitial site was computed relative to the mean energy of the system when 40 Ar resides in a regular interlayer site. The simulations reveal that when 40 Ar is in the interstitial site, i.e., during the first 50 ps, the system has a mean total energy of about 100 kcal.mol -1 higher than when 40 Ar occupies a regular interlayer site (Fig. 2). This energy increase is due to the distortion that 40 Ar causes to the muscovite structure because of its large size. Simulation snapshots (not shown) reveal that when 40 Ar is located in the interstitial position, it repeals its neighboring basal oxygens, increasing locally the interlayer spacing by ca. 15 %. This is a highly unfavourable configuration energetically, precluding protracted 40 Ar residence in the interstitial site and making unlikely a reincorporation of 40 Ar into the interstitial sites once it has lost its initial recoil velocity. This rules out interstitial mechanism as a viable mechanism for 40 Ar diffusion in muscovite on geologic timescales.

Pathways and 40 Ar diffusion barriers in muscovite

The foregoing simulations suggest that the easiest way for 40 Ar to move within the muscovite structure without causing major distortions is by jumping into neighboring interlayer vacancies. Vacancies are always present in minerals. They form at high temperature to reduce the Gibbs free energy of the crystal [START_REF] Mehrer | Diffusion in solids: fundamentals, methods, materials, diffusion-controlled processes[END_REF]. Beside those intrinsic vacancies, additional vacancies can be introduced in the mica interlayer region via the pyrophyllitic substitution Si + □ → IV Al + K [START_REF] Tardy | An ideal solid solution model for calculating solubility of Clay miner[END_REF][START_REF] Aagaard | Activity/composition relations among silicates and aqueous solutions: II. Chemical and thermodynamic consequences of ideal mixing of atoms on homological sites in montmorillonites, illites, and mixed-layer clays[END_REF][START_REF] Helgeson | Activity/composition relations among silicates and aqueous solutions; I, Thermodynamics of intrasite mixing and substitutional order/disorder in minerals[END_REF][START_REF] Wang | Non-stoichiometry of interlayer cations in micas from low-to middle-grade metamorphic rocks in the Ryoke and the Sanbagawa belts[END_REF][START_REF] Baldelli | Ferrimuscovite and celadonite substitutions in muscovite from Fe 3+ -rich low-grade psammitic rocks (Northern Apennines, Italy)[END_REF][START_REF] Wunder | Interlayer vacancy characterization of synthetic phlogopitic micas by IR spectroscopy[END_REF].

Their abundance estimated from the commonly measured stochiometric deficit of interlayers cations (up to ca. 10 %) suggests that they can be involved in the motion of interlayer atoms. Therefore, we performed MD simulations in which we initially positioned 40 Ar in an interlayer K + site and added interlayer vacancies in its vicinity. Simulations were performed in the NPT ensemble at temperatures of 650 to 750 • C and P = 5 kbar, for time durations up to 25 ns. No migration could be observed; neither 40 Ar nor K + atoms jumped into the vacancies. This immobility suggests high migration energy barriers, making diffusion through a vacancy mechanism an infrequent event, unobservable on the timescales of classical MD.

To overcome the timescale limitation, we studied 40 Ar diffusion using TST. This requires the migration paths and the migration energies of diffusing atom to be determined first, unlike classical MD simulations which only require initial positions and velocities of atoms. The 40 Ar migration energies from its interlayer sites to neighboring vacancies were computed using the NEB method. Since the value of the migration energy in Eqs. ( 6) and ( 8) corresponds to the migration barrier of a system at the absolute zero, we first equilibrated the system via NPT simulations at P = 0 kbar and T = -273.1 • C (i.e., 0.5 K since there is no dynamics at 0 K). The migration of 40 Ar in the muscovite interlayer region depends on the presence of a vacancy in its immediate vicinity. Since vacancies can initially form far from the 40 Ar site, a motion of potassium atoms along the interlayer region is required for a vacancy to end next to 40 Ar so the K + migration barriers must be also determined. Furthermore, because vacancies move essentially either individually or by pairs (vacancy clusters involving more than two vacancies being thermodynamically less stable; [START_REF] Mehrer | Diffusion in solids: fundamentals, methods, materials, diffusion-controlled processes[END_REF], we considered the migration barriers for these two configurations, i.e., the migration energies of 40 Ar and K + from their sites towards a simple vacancy (monovacancy mechanism, labelled "□"; Fig. 3a),

or towards a vacancy adjacent to another vacancy (divacancy mechanism, labelled "□□"; Fig. 3b). As the environment around each interlayer site is variable, in particular because of the irregular distribution of Al in the tetrahedral layers, the statistical (probabilistic) calculations must take into account the configuration of atoms in the vicinity of 40 Ar and K + .

We positioned thus the 40 Ar and K + in all the interlayer sites of the model and computed the migrations energies towards all possible neighboring sites. The charge deficit induced by the incorporation of 40 Ar and vacancies into the interlayer region was compensated by substituting Al atoms for Si in tetrahedral sites far from the 40 Ar migration path (some unit cells away). Our results show that, in the case of the divacancy mechanism, the migration barriers related to the "black" jumps (Fig. 3b) are identical to the energy barriers associated with the monovacancy mechanism. Therefore, we report below 216 migration energies for the monovacancy mechanism (all possible combinations), and 432 for the divacancy mechanism (corresponding only to the "blue" jumps).

40 Ar migration barriers

The energy barriers associated with the interlayer migration of 40 Ar towards neighboring vacancies display wide distributions (Fig. 4a,b). The migration barriers associated with the monovacancy mechanism yield a mean value of 71 ± 5 kcal.mol -1 (1σ), while those related to the divacancy mechanism show a slightly lower but statistically indistinguishable mean value of 69 ± 6 kcal.mol -1 (1σ). These values are similar to the energy barrier of 69 kcal.mol -1 computed using DFT for the migration of Ar in lizardite [START_REF] Wang | Diffusion of noble gases in subduction zone hydrous minerals[END_REF], which consists of a stacking of tetrahedral and octahedral layers like muscovite [START_REF] Mellini | The crystal structure of lizardite 1T: hydrogen bonds and polytypism[END_REF]. The similarity between the migration barriers related to both mechanisms is illustrated by the variation of the 40 Ar potential energy along its migration path between two neighboring stable sites (corresponding to images "1" and "8"; Fig. 4c). In this example, the energy barrier related to the migration of 40 Ar towards a single vacancy is 65 kcal.mol -1 while that associated with the jump of 40 Ar into a pair of vacancies is 61 kcal.mol -1 . In both cases, the forward and backward migration barriers are similar, and this has been observed for all computed paths.

The minimum energy paths followed by 40 Ar during these jumps (Fig. 4d) reveal that when jumping into a single vacancy, 40 Ar follows an almost rectilinear trajectory, passing directly below the basal oxygens crossed along the path. When it migrates into a pair of vacancies, 40

Ar follows a similar trajectory but avoids the basal oxygens, passing quite closer to the vacancy not involved in the migration path. The slight difference in trajectories shows that the difference in magnitude between the respective migration barriers is directly related to the presence or not of a K + around the pathway.

The migration barriers are anticorrelated to the total electrostatic energies calculated for each Ar-vacancy configuration in the muscovite model (Fig. 4e). Ruiz [START_REF] Pestana | Direct exchange mechanism for interlayer ions in non-swelling clays[END_REF] also report such anticorrelation between the potassium migration energies and the cohesive energies of K + at stable sites in illite. These observations suggest that the observed variability in the magnitude of the energy barriers in our model is related to the nonuniform distribution of electrostatic charges around the migration paths, which results from the irregular distribution of Al atoms in the tetrahedral layers. The highest migration barriers are related to paths with less tetrahedral Al around, in contrast to the lowest energies obtained for those surrounded by more tetrahedral Al. In addition to the irregular distribution of tetrahedral Al atoms, variations in the magnitude of migrations barriers also depend to some extent on the relative position of the octahedral vacancies with respect to the migration path. In 2M 1 muscovite, the projection of octahedral cations and interlayer atoms (K + ) on the (001) plane

shows that the octahedral vacancies are located directly above the migration paths forming angles of ± 60 • and ± 120 • with the a axis. We observe that for migration paths with the same number of surrounding tetrahedral Al, migration energies are systematically higher for paths parallel to the a axis. These observations suggest a non-negligible effect of the mica structure on the intrinsic diffusion of 40 Ar within the interlayer region, as suggested by previous observations and empirical predictions [START_REF] Scaillet | Mg/Fe and [(Mg,Fe)Si-Al 2 ] compositional control on argon behaviour in high-pressure white micas: A 40 Ar/ 39 Ar continuous laser-probe study from the Dora-Maira nappe of the internal Western Alps, Italy[END_REF][START_REF] Dahl | The crystal-chemical basis for Ar retention in micas: inferences from interlayer partitioning and implications for geochronology[END_REF].

K + migration barriers

As with 40 Ar, the K + migration barriers define wide distributions (Fig. 5a,b). Migrations The K + energy barriers associated with the divacancy mechanism yield a mean value of 66 ± 10 kcal.mol -1 (1σ). This value is significantly lower than the mean energy barrier related to the migration of K + into a single vacancy. This pronounced difference is illustrated by the variation of the potential energy along the K + migration path between two adjacent stable sites (Fig. 5c). In this example, the K + energy barrier associated with the monovacancy mechanism is 70 kcal.mol -1 and the K + energy barrier related to the divacancy mechanism 55 kcal.mol -1 . The minimum energy paths followed by K + during these migrations show that, in both cases, the potassium atom follows a curvilinear trajectory avoiding the basal oxygen atoms along the path (Fig. 5d). However, in the case of the divacancy mechanism, potassium passes further away from the oxygen atoms than in the monovacancy mechanism.

This suggests that the decrease in the migration energy is due to the ability of K + ions to minimize interactions with the oxygen atoms lying on their path. Moreover, while 40 Ar is at the same energy level in two neighboring vacancies, potassium shows an energy difference between the two sites (Fig. 5c). This is because potassium (charged) is more sensitive to the local charge distribution than 40 Ar (neutral). Despite these structural-charge differences, the K + and 40 Ar energy barriers are positively correlated (Fig. 5e), suggesting a similar influence of the local environment on the migration energy landscape.

40 Ar diffusion coefficients

In this section, we calculate the 40 Ar diffusion coefficients in muscovite for both the monovacancy and the divacancy mechanisms by TST (Eq. ( 8)). Since the investigated migrations occur along the interlayer region, we perform calculations for 2D diffusion in a hexagonal lattice. As given above, Eq. ( 8) allows the analytical derivation of the diffusion coefficient for a pristine and homogenous crystal, i.e., a crystal in which the migration barrier is the same for all equivalent jumps. In our case, because of the irregular distribution of Al atoms in the tetrahedral layers, the energy barriers of interlayer atoms vary from one direction to another, and from one interlayer site to another. To take into account the spatial variation of the migration barrier, we consider for each diffusing atom an overall effective migration energy, E eff , defined as:

exp E eff RT = n i=1 1 n exp E i RT (10)
where E i is the effective energy barrier required for an atom to leave an interlayer site i and n is the total number of stable sites in the interlayer region. The effective energy barrier for each site (E i ) is derived from:

exp

-E i RT = m j=1 1 m exp -E ij RT (11)
in which E ij is the migration energy for leaving the interlayer site i in the direction j, and m is the number of all possible directions (m = 6 in a hexagonal lattice).

Eqs ( 10) and ( 11) are justified by making an analogy with an electric circuit. When different paths from a given site are possible, the conductances of the paths add up. As the conductance is proportional to the Boltzmann factor, this leads to Eq. ( 11). A consequence is that the path having the lowest migration energy dominates. When considering a set of sites, they are travelled in sequence so that the resistances add up as expressed by Eq. ( 10).

In particular a site having a high migration energy will slow down the entire diffusion process, this is the trapping effect. Using the energy barriers calculated in section 5, we obtain effective vacancy is present in the interlayer region, it is filled by the jump of a neighboring K + which, in turn, releases a vacancy in its initial site (Fig. 3a). The displacement vectors have as components:

- → V 0 =   1 2 a √ 3 2 a   - → V 1 =   -1 2 a √ 3 2 a   - → V 2 =   -a 0   - → V 3 =   -1 2 a - √ 3 2 a   - → V 4 =   1 2 a - √ 3 2 a   - → V 5 =   a 0   (12)
and the variance of the displacement of the vacancy is:

∆x 2 = ∆y 2 = a 2 2 (13)
We denote by γ K □ the frequency of the jump of one of the neighbouring K + towards a monovacancy (noted K + ⇆ □), ν the attempt frequency, and E K □ eff the migration barrier of K + . According to TST, the jump frequency is given by:

γ K □ = ν exp - E K □ eff RT (14)
The timelife of the vacancy (i.e., the average time for a successful displacement) is:

∆t K = 1 6γ K □ ( 15
)
and the 2D diffusion coefficient of the vacancy is:

D □ = ∆x 2 + ∆y 2 4∆t K = 3 2 a 2 ν exp - E K □ eff RT ( 16
)
where a is the distance between two neighboring interlayer sites (a = 5.2 Å). At T = 700 • C, the diffusion coefficient of a single vacancy is D □ = 2.88 × 10 -19 cm 2 .s -1 .

40 Ar diffusion in a monovacancy flux

Each time a vacancy gets next to a 40 Ar atom it is surrounded by five K + ions. Interchange between the 40 Ar and vacancy positions (noted 40 Ar ⇆ □) occurs due to thermal agitation and the lower barrier potential of 40 Ar ⇆ □ relative to K + ⇆ □. The average duration of these Ar-vacancy oscillations is therefore:

∆t osc = 1 5γ K □ (17)
During this time, many oscillations will have occurred with equal probability 0.5 of ending up either in the initial configuration (no exchange) or in effective migration of 40 Ar into the vacancy. The oscillations stop when one of the five K + neighbors jumps into the vacancy.

For three of these K + , a successful jump will remove the vacancy away from 40 Ar, stopping the 40 Ar ⇆ □ oscillations. For the other two, the vacancy stays in the vicinity of 40 Ar and the oscillations can resume among the newly created Ar-vacancy pair.

Summing up, when a vacancy gets close to 40 Ar, after an average time ∆t osc , four different configurations are possible:

-The vacancy has gone away leaving 40 Ar in its initial position with a probability of 3/10.

-The vacancy has moved into one of the two neighboring positions of 40 Ar with a probability of 2/10. Oscillations resume among the newly-created 40 Ar-vacancy pair.

-The 40 Ar ⇆ □ interchange is successful and the vacancy moves further away by more than one lattice position with probability of 3/10.

-The 40 Ar ⇆ □ interchange has occurred but leaving the vacancy next to 40 Ar with probability of 2/10. Oscillations resume again among the new Ar-vacancy pair.

The average number of jumps ⟨n s ⟩ needed to move the vacancy away from 40 Ar can be calculated by combining the respective probabilities:

⟨n s ⟩ = 3 5 ∞ n=1 n 2 5 n-1 = 5 3 1 (1 -2/5) 2 = 5 3 (18)
According to Eq. ( 18), the average residence time of the vacancy in the vicinity of 40 Ar (the timelife of this configuration) is:

∆t o = 5 3 ∆t osc = 1 3γ K □ (19)
During this time interval, several trains of oscillations will have occured. The calculation of the variance of the mean displacement of 40 Ar is complicated by the fact the successive displacements are correlated. However, KMC simulations show that the positional variance of the displacement of □ along either -→ x or -→ y is completely defined by the number of jumps n s of the vacancy remaining in the neighborhood:

⟨∆x 2 ⟩ = ⟨∆y 2 ⟩ = n s + 1 8 a 2 (20)
By summing over the n s probabilities, we obtain the average variance resulting in the nearestneighbor vacancy configuration:

⟨∆x 2 ⟩ = ⟨∆y 2 ⟩ = a 2 3 (21)
Note that, despite its higher energy barrier, the variance is stochastically dominated by the bulk motion of K + via the K + ⇆ □ interchange (vacancy diffusion), not the mobility of 40 Ar with respect to the vacancy. The vibrational motion of 40 Ar just ensures that it will take (with probability 0.5) a quantum step away from its initial position when the vacancy passes close by. Note also that this is the variance induced by the proximity of a vacancy. Most of the time, no vacancy hangs around 40 Ar which cannot move then (the time goes on with the variance remaining unchanged).

We now inquire into the bulk vacancy flux that governs the overall mobility of 40 Ar through the muscovite interlayer. Consider a portion of lattice containing N regular interlayer sites with a single 40 Ar atom and a density of vacancies ρ. The number of neighboring sites around a vacancy is six and the number of vacancies N ρ. At each vacancy jump, the probability of starting from either a normal "bulk" position or next to 40 Ar is proportional to their respective populations. The probability of starting from a neighboring site must be weighted by a factor of 5/6 because only five directions allow it to be reached, instead of six for a normal site. This quantitative reduction has been tested by KMC simulation (Appendix A). Finally, during a vacancy jump the probability of starting from a bulk site is P b = (N -6)/N while that of jumping from the vicinity of 40 Ar is P v = 6/N (5/6) = 5/N .

Consider now a trajectory of m jumps of the N ρ vacancies. Among these m jumps mP b will start from a bulk site whose lifetime is t b = 1/(6γ K □ ) and mP v will come from neighboring sites with lifetime t v = 1/(5γ K □ ). The total duration of these trajectories of m jumps is therefore:

t = m 1 - 6 N 1 6γ K □ + m 5 N 1 5γ K □ = m 6γ K □ (22)
The average number of vacancy jumps ending next to 40 Ar is:

m v = m 5 N N ρ = 5ρ m (23)
We have seen that in the direct-neighbor configuration (Fig. 3a), on average 5/3 jumps move the vacancy away from 40 Ar (Eq. ( 18)). The number of such moves in a bulk vacancy flux is:

m ′ v = 3ρ m (24)
According to Eq. ( 21), each of these moves has variance a 2 /3 so that on average:

⟨∆x 2 ⟩ = ⟨∆y 2 ⟩ = ρa 2 m = 6ργ K □ a 2 t (25)
The diffusion coefficient for 40 Ar is thus:

D Ar □ = ⟨x 2 ⟩ + ⟨y 2 ⟩ 4 t = 3 a 2 γ K □ ρ = 3 a 2 ν ρ exp - E K □ eff RT (26)
The temperature dependence of this diffusion mechanism is reported in the Arrhenius diagram (Fig. 6). For example, at T = 700 • C and ρ = 0.1, the diffusion coefficient is D Ar □ = 5.76 × 10 -20 cm 2 .s -1 . The diffusion parameters associated with the monovacancy mechanism are E = 75 kcal.mol -1 and D 0 = 4 × 10 -3 cm 2 .s -1 .

6.2. Divacancy mechanism 6.2.1. Diffusion of divacancies in a hexagonal lattice of K +

In the case of a divacancy mechanism, K + atoms are adjacent either to only one vacancy or both (Fig. 3b). In the latter configuration, the jumps (blue arrows) have a lower effective energy barrier (58 kcal.mol -1 ) and therefore a higher probability of occurrence than those of the first configuration (in black), whose migration barriers are identical to those associated with the monovacancy mechanism (75 kcal.mol -1 ). When one of these "blue" migrations occurs, a K + atom fills one vacancy and a new vacancy is created in its initial position.

This new vacancy lies in one of the sites neighboring the unoccupied vacancy of the initial pair, and both form a new vacancy pair. Therefore, one can expect a correlated mobility of divacancies along the interlayer region.

The two pertinent quantities to calculate here are the diffusion coefficient of the divacancy and its lifetime. The divacancy is characterized by its position (barycenter of the vacancy pair), the direction passing through the vacancy pair, and the angle made by the direction with the x-axis which can be 0 • (Fig. 3b) or ± 60 • . Each atom jump filling a vacancy produces a virtual displacement and the apparent rotation of the divacancy direction. Let us calculate the average variance of such a displacement.

When the starting angle is 0

• the jumps - → δ 0 , - → δ 1 , - → δ 2 , - → δ 3
produce respectively the angles 60 • , 60 • , -60 • , -60 • , and the divacancy displacement vectors:

- → V 0 =   1 4 a √ 3 4 a   - → V 1 =   -1 4 a - √ 3 4 a   - → V 2 =   1 4 a - √ 3 4 a   - → V 3 =   -1 4 a √ 3 4 a   (27)
By symmetry, these jumps are equiprobable.

For a starting angle of 60 • , the divacancy displacement vectors -→ V i are just rotated by this angle :

- → V i ′ =   1 2 - √ 3 2 √ 3 2 1 2   - → V i ( 28
)
and:

- → V 0 ′ =   -1 4 a √ 3 4 a   - → V 1 ′ =   1 4 a - √ 3 4 a   - → V 2 ′ =   1 2 a 0   - → V 3 ′ =   -1 2 a 0   (29)
For the -60 • configuration, the corresponding rotation is:

- → V 0 ′′ =   1 2 a 0   - → V 1 ′′ =   -1 2 a 0   - → V 2 ′′ =   -1 4 a - √ 3 4 a   - → V 3 ′′ =   1 4 a √ 3 4 a   (30)
Whatever the starting configuration, it can be shown that ⟨x⟩ = ⟨x⟩ = 0. For a given trajectory, the orientations are equiprobable. For all jumps keeping the vacancy pair connected the displacements

- → V i , - → V i ′ and - → V i
′′ occur with equal probability of 1/12. One may verify that ⟨x • y⟩ = 0, indicating that the displacements along either x or y are uncorrelated which is a consequence of hexagonal symmetry. In terms of variance of displacements we have:

⟨x 2 ⟩ = ⟨y 2 ⟩ = 1 8 a 2 (31)
due again to the symmetry in x and y.

Denote by γ K □ the jump frequency corresponding to a "black" migration and γ K □□ for the "blue" migration (Fig. 3b). The residence time of the divacancy at any given position is:

t r = 1 4γ K □□ + 6γ K □ (32)
The average migration time after m jumps is:

⟨t⟩ = m t r = m 4γ K □□ + 6γ K □ (33)
and the variance of the displacement:

⟨x 2 ⟩ = ⟨y 2 ⟩ = m a 2 8 ( 34
)
From these, we deduce the divacancy diffusion coefficient:

D □□ = 1 4 a 2 γ K □□ + 3 2 γ K □ D □□ = 1 4 a 2 ν exp - E K □□ eff RT 1 + 3 2 exp E K □□ eff -E K □ eff RT ( 35
)
where E K □□ eff is the K + migration barrier towards a divacancy. The second exponential in parentheses (rightmost term) is much smaller than unity and the diffusion coefficient reduces to:

D □□ = 1 4 a 2 ν exp - E K □□ eff RT (36)
At T = 700 • C, the diffusion coefficient associated with the divacancy migration is D □□ = 3.89 × 10 -16 cm 2 .s -1 , three orders of magnitude faster than the monovacancy mechanism.

40 Ar diffusion in a divacancy flux

When a divacancy lies in the vicinity of 40 Ar, the latter can be adjacent to either one or both vacancies and the probability that one is filled by 40 Ar depends on the specific nearest-neighbor configuration. Denote by E Ar □ eff the 40 Ar energy barrier faced by the "black" transitions and by E Ar □□ eff that associated with "blue" migrations. When the 40 Ar is adjacent to only one of the two vacancies (Fig. 3b) the probability of jump before a move of the divacancy is:

P Ar □ = exp - E Ar □ eff RT exp - E Ar □ eff RT + 5 exp - E K □ eff RT + 2 exp - E K □□ eff RT = exp - E Ar □ eff -E K □□ eff RT exp - E Ar □ eff -E K □□ eff RT + 5 exp - E K □ eff -E K □□ eff RT + 2 (37)
The exponentials in the denominator are much smaller than one because of the difference in the energy barriers, and Eq. ( 37) simplifies to

P Ar □ = 1 2 exp - E Ar □ eff -E K □□ eff RT (38)
On the other hand, there are two positions where the 40 Ar is close to both vacancies (Fig. 3b) and the probability of jumping into one of them before the divacancy move is:

P Ar □□ = exp - E Ar □□ eff RT exp - E Ar □ eff RT + 6 exp - E K □ eff RT + exp - E K □□ eff RT (39)
As before, Eq. ( 39) simplifies to

P Ar □□ = exp - E Ar □□ eff -E K □□ eff RT (40)
The relative importance of the probabilities Eqs. ( 38) and ( 40) depends on the frequency of each configuration.

Denote by ρ ′ the density of divacancies in the interlayer region, defined as the number of divacancies divided by the number of interlayer sites in the lattice. Consider a trajectory of m jumps of the divacancies. For 6mρ ′ positions the 40 Ar will be close to one only vacancy, and for mρ ′ it will be close to both vacancies. The above expressions can be deduced by weighing the different orientations and by observing that the double-vicinity of the 40 Ar can be reached only from one side, and they have been confirmed by running dedicated KMC simulations (Appendix A.).

In the bulk divacancy flux, the residence time of a divacancy is given by Eq. ( 32) in which γ K □ can be neglected. This residence time is nearly the same when a 40 Ar is close to only one vacancy, but it doubles when it neighbors the two vacancies. Therefore the time for a trajectory of m jumps is:

t = m 4γ K □□ (1 + ρ ′ ) = m 4ν (1 + ρ ′ ) exp E K □□ eff RT (41)
For 6mρ ′ configurations, 40 Ar sees only one vacancy and the jump probability is given by Eq.

(38). For mρ ′ configurations it neighbors the two vacancies and has a jump probability given by Eq. ( 40). For any jump the variance in coordinates is given by Eq. ( 13). In the end, for m jumps the variance is:

⟨x 2 ⟩ = ⟨y 2 ⟩ = mρ ′ a 2 2 3 exp - E Ar □ eff -E K □□ eff RT + exp - E Ar □□ eff -E K □□ eff RT (42)
From Eq. ( 41) and ( 42) we can deduce the diffusion coefficient for 40 Ar:

D Ar □□ = ρ ′ 1 + ρ ′ a 2 ν 3 exp - E Ar □ eff RT + exp - E Ar □□ eff RT (43)
in which the first term inside the brackets is small compared to the last one.

It remains to determine the density of divacancies ρ ′ . It can be expressed with the density of monovacancies ρ. The divacancies are produced when two monovacancies hit each other, because their coherent motion is energetically favored. On the other hand they decay when one of its vacancies is filled by a neighboring K + (black jumps; Fig. 3). This is a dynamical equilibrium reached when the rate of production equates the rate of decay.

Consider N regular interlayer sites. The number of monocancies is N ρ. At any jump of a vacancy, the probability to make a collision (i.e. to reach the neighborhood of another vacancy) is proportional to ρ, so that the number of produced divacancies is kN ρ 2 , where k is a number to be determined. A dedicated random walk simulation shows that k = 2.5. The duration of each step is 1/(6γ K □ ) giving a production rate of 15N ρ 2 γ K □ . The decay rate is 6N ρ ′ γ K □ . By equating the two rates one gets ρ ′ = 2.5ρ 2 .

D Ar □□ = 5 2 ρ 2 a 2 ν 3 exp - E Ar □ eff RT + exp - E Ar □□ eff RT (44)
At T = 700 • C and for ρ = 0.1, the 40 Ar diffusion coefficient associated with the divacancy mechanism is D Ar □□ = 1.16 × 10 -18 cm 2 .s -1 , two orders of magnitude higher than the diffusivity of the monovacancy mechanism. The 40 Ar diffusion parameters deduced from the Arrhenius diagram (Fig. 6) are E = 66 kcal.mol -1 and D 0 = 6 × 10 -4 cm 2 .s -1 . A summary listing of the final estimates is provided in Table 3.

Discussion

Comparison with previous 40 Ar diffusion estimates

As mentioned above, only a handful of studies have attempted to quantify 40 Ar diffusion in muscovite [START_REF] Robbins | Radiogenic argon diffusion in muscovite under hydrothermal conditions[END_REF][START_REF] Hames | An empirical evaluation of the argon diffusion geometry in muscovite[END_REF][START_REF] Kirschner | Staircase 40 Ar/ 39 Ar spectra of fine-grained white mica: Timing and duration of deformation and empirical constraints on argon diffusion[END_REF]Lister andBaldwin, 1996, 2016;[START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF][START_REF] Forster | 40 Ar/ 39 Ar geochronology and the diffusion of 39 Ar in phengite-muscovite intergrowths during step-heating experiments in vacuo[END_REF]. While

40

Ar diffusivity estimates inferred by these studies are broadly consistent, they considerably differ from the present results. Our MD and TST estimates are up to five orders of magnitude lower than the corresponding Arrhenian laws (Fig. 6). Note that our results are calculated for muscovite at P = 0 kbar. Assuming a pressure effect on diffusivity of 40 Ar in muscovite [START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF] would enhance the disparity even further.

Hydrothermal experiments have been suspected to overestimate 40 Ar diffusivity in natural muscovite for quite a while now [START_REF] Monié | Preservation of Hercynian 40 Ar/ 39 Ar ages through high-pressure low-temperature Alpine metamorphism in the Western Alps[END_REF][START_REF] Scaillet | 40 Ar/ 39 Ar laser-probe dating by step heating and spot fusion of phengites from the Dora Maira nappe of the Western Alps, Italy[END_REF][START_REF] Scaillet | Mg/Fe and [(Mg,Fe)Si-Al 2 ] compositional control on argon behaviour in high-pressure white micas: A 40 Ar/ 39 Ar continuous laser-probe study from the Dora-Maira nappe of the internal Western Alps, Italy[END_REF][START_REF] Monié | 40 Ar/ 39 Ar dating in coesite-bearing and associated units of the Dora Maira massif[END_REF][START_REF] Arnaud | Evidence for excess argon during high pressure metamorphism in the Dora Maira Massif (western Alps, Italy), using an ultra-violet laser ablation microprobe 40 Ar-39 Ar technique[END_REF][START_REF] Scaillet | Excess 40 Ar transport scale and mechanism in high-pressure phengites: A case study from an eclogitized metabasite of the Dora-Maira nappe, western Alps[END_REF][START_REF] Vincenzo | Evidence of a close link between petrology and isotope records: constraints from SEM, EMP, TEM and in situ 40 Ar-39 Ar laser analyses on multiple generations of white micas (Lanterman Range, Antarctica)[END_REF][START_REF] Rapaport | The art of molecular dynamics simulation[END_REF], 2006;[START_REF] Balogh | Argon and fission track dating of Alpine metamorphism and basement exhumation in the Sopron Mts.(Eastern Alps, Hungary): thermochronology or mineral growth?[END_REF][START_REF] Hames | Single-crystal 40 Ar/ 39 Ar age variation in muscovite of the Gassetts Schist and associated gneiss, Vermont Appalachians[END_REF][START_REF] Itaya | Regional-Scale Excess Ar wave in a Barrovian type metamorphic belt, eastern Tibetan Plateau[END_REF][START_REF] Allaz | The effects of retrograde reactions and of diffusion on 40 Ar-39 Ar ages of micas[END_REF][START_REF] Beltrando | Preservation of sub-microscopic structural relicts in micas from the Gran Paradiso Massif (Western Alps): Implications for 40 Ar-39 Ar geochronology[END_REF][START_REF] Halama | Effects of fluid-rock interaction on 40 Ar/ 39 Ar geochronology in high-pressure rocks (Sesia-Lanzo Zone, Western Alps)[END_REF][START_REF] Schertl | Tracking the incidence of excess argon in white mica Ar-Ar data from UHP conditions to upper crustal levels in the Dora-Maira Massif, Western Alps[END_REF][START_REF] Laurent | 40 Ar behaviour and exhumation dynamics in a subduction channel from multi-scale 40 Ar/ 39 Ar systematics in phengite[END_REF]. As emphasized before, most of these studies relate to high-pressure (phengitic) muscovite. The high retentivity (i.e., low 40 Ar diffusivity) of muscovite in these conditions has long been traced to crystal-structural effects suspected to affect the interlayer at increasingly higher pressure [START_REF] Scaillet | Mg/Fe and [(Mg,Fe)Si-Al 2 ] compositional control on argon behaviour in high-pressure white micas: A 40 Ar/ 39 Ar continuous laser-probe study from the Dora-Maira nappe of the internal Western Alps, Italy[END_REF][START_REF] Scaillet | K-Ar ( 40 Ar/ 39 Ar) geochronology of ultrahigh pressure rocks[END_REF]. It is also supported by the diffusion data of [START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF] which predict lower diffusivities at high pressure. Evidence for a pronounced 40 Ar retentivity in muscovite has been also reported in other settings as well (e.g. Di [START_REF] Vincenzo | The relationship between tectono-metamorphic evolution and argon isotope records in white mica: constraints from in situ 40 Ar-39 Ar laser analysis of the Variscan basement of Sardinia[END_REF][START_REF] Hames | Single-crystal 40 Ar/ 39 Ar age variation in muscovite of the Gassetts Schist and associated gneiss, Vermont Appalachians[END_REF][START_REF] Itaya | Regional-Scale Excess Ar wave in a Barrovian type metamorphic belt, eastern Tibetan Plateau[END_REF][START_REF] Allaz | The effects of retrograde reactions and of diffusion on 40 Ar-39 Ar ages of micas[END_REF], suggesting that sluggish 40 Ar kinetics in natural muscovite is more the norm than currently acknowledged by existing experimental data.

Field-based calibration of 40 Ar retention kinetics in muscovite is scarce and highly contextspecific (e.g., Wijbrans andMcDougall, 1986, 1988;[START_REF] Blanckenburg | Time calibration of a PTpath from the Western Tauern Window, Eastern Alps: the problem of closure temperatures[END_REF][START_REF] Kirschner | Staircase 40 Ar/ 39 Ar spectra of fine-grained white mica: Timing and duration of deformation and empirical constraints on argon diffusion[END_REF]. So far, such efforts have been taken to indicate moderate or intermediate 40 Ar retention properties of muscovite relative to other isotope-mineral pairs (by order of decreasing retentivity: U/Pb-zircon > Rb/Sr-muscovite > K/Ar-amphibole > K/Ar-muscovite > K/Ar-biotite ≈ K/Ar-alkali feldspar). Mounting evidence, in particular from grainscale Rb/Sr and laser-based 40 Ar/ 39 Ar disequilibrium data (e.g. [START_REF] Bröcker | The geological significance of 40 Ar/ 39 Ar and Rb-Sr white mica ages from Syros and Sifnos, Greece: a record of continuous (re)crystallization during exhumation?[END_REF][START_REF] Cliff | Geochronological challenges posed by continuously developing tectonometamorphic systems: insights from Rb/Sr mica ages from the Cycladic Blueschist Belt, Syros (Greece)[END_REF][START_REF] Laurent | 40 Ar behaviour and exhumation dynamics in a subduction channel from multi-scale 40 Ar/ 39 Ar systematics in phengite[END_REF], shows that Ar-Rb-Sr retention systematics in natural muscovite are actually quite similar [START_REF] Glodny | Precise eclogitization ages deduced from Rb/Sr mineral systematics: the Maksyutov complex, Southern Urals, Russia[END_REF][START_REF] Meleshyn | Two-dimensional ordering of water adsorbed on a mica surface at room temperature[END_REF][START_REF] Bosse | Rb-Sr and 40 Ar/ 39 Ar ages in blueschists from the Ile de Groix (Armorican Massif, France): Implications for closure mechanisms in isotopic systems[END_REF][START_REF] Scibiorski | Cryptic excess argon in metamorphic biotite: Anomalously old 40 Ar/ 39 Ar plateau dates tested with Rb/Sr thermochronology and Ar diffusion modelling[END_REF], pointing to otherwise more sluggish kinetics than suggested by the above trend for elements that are either ionically bonded (Rb, Sr) and/or simply too big (Ar)

to easily jump from a lattice position to the next. This is most clearly illustrated by peakmetamorphic 40 Ar/ 39 Ar retention ages commonly preserved in muscovite despite deformation and overprinting [START_REF] Putlitz | Prograde mica 40 Ar/ 39 Ar growth ages recorded in high pressure rocks (Syros, Cyclades, Greece)[END_REF][START_REF] Beaudoin | Situ and Step-Heating 40 Ar/ 39 Ar Dating of White Mica in Low-Temperature Shear Zones (Tenda Massif[END_REF][START_REF] Laurent | 40 Ar behaviour and exhumation dynamics in a subduction channel from multi-scale 40 Ar/ 39 Ar systematics in phengite[END_REF]. This is also emphasized by congruent 40 Ar/ 39 Ar muscovite and U/Pb titanite ages metastably preserved over map-scales in excess of 450-750 km 2 in dry granulite-amphibolite crust exposed to T ⩾ 700 • C for tens of Myr [START_REF] Spencer | Campaign-style titanite U-Pb dating by laser-ablation ICP: Implications for crustal flow, phase transformations and titanite closure[END_REF].

Using the diffusion data retrieved from our modeling (E = 66 kcal.mol -1 and D 0 = 6 × 10 -4 cm 2 .s -1 ), nominal closure temperatures can be calculated via Eq. ( 1) for various cooling rates and grain sizes. For cooling rates ranging from 1 to 100 • C.Ma -1 and grain size varying between 0.1-1.0 mm, this returns closure temperatures in the range 560-800 • C (Fig. 7). These are significantly (up to 200 • C) higher than predicted by current experimental estimates [START_REF] Hames | An empirical evaluation of the argon diffusion geometry in muscovite[END_REF][START_REF] Kirschner | Staircase 40 Ar/ 39 Ar spectra of fine-grained white mica: Timing and duration of deformation and empirical constraints on argon diffusion[END_REF][START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF] and close to the upper stability limit of muscovite (∼ 700 • C; Massonne andSchreyer, 1987, 1989), which by the way is significantly lower in presence of quartz [START_REF] Althaus | An experimental re-examination of the upper stability limit of muscovite plus quartz[END_REF]. This opens up the prospect that crystallization ages may be more likely than cooling ages in natural muscovite.

However, empirical evidence also abounds to suggest that intracrystalline 40 Ar mobility may be enhanced by other factors beyond purely diffusion-driven rates such as predicted by our model. Our model calculations apply to ideal (undefective) muscovite. Natural micas are far from perfect and always contain microstructural defects such as dislocations, stacking faults, microcleavages, etc., which can develop during crystal growth, or more commonly, in response to deformation (e.g. Bell andWilson, 1977, 1981;[START_REF] Amouric | Growth and deformation defects in phyllosilicates as seen by HRTEM[END_REF][START_REF] Kronenberg | Basal slip and mechanical anisotropy of biotite[END_REF][START_REF] Shea | Rheology and deformation mechanisms of an isotropic mica schist[END_REF][START_REF] Dunlap | Neocrystallization or cooling? 40 Ar/ 39 Ar ages of white micas from low-grade mylonites[END_REF][START_REF] Beaudoin | Situ and Step-Heating 40 Ar/ 39 Ar Dating of White Mica in Low-Temperature Shear Zones (Tenda Massif[END_REF]. Such extended microstructural defects act as high diffusivity pathways that can considerably enhance the mobility of diffusing atoms [START_REF] Lee | Multipath diffusion in geochronology[END_REF]. Sizeable extents of 40 Ar loss have been documented in naturally and experimentally deformed muscovite in connection with microstructural defects [START_REF] Kramar | Heterogeneous 40 Ar* distributions in naturally deformed muscovite: in situ UV-laser ablation evidence for microstructurally controlled intragrain diffusion[END_REF][START_REF] Kramar | Stacking fault-enhanced argon diffusion in naturally deformed muscovite[END_REF][START_REF] Mulch | In-situ UV-laser 40 Ar/ 39 Ar geochronology of a micaceous mylonite: an example of defect-enhanced argon loss[END_REF][START_REF] Cosca | 40 Ar* loss in experimentally deformed muscovite and biotite with implications for 40 Ar/ 39 Ar geochronology of naturally deformed rocks[END_REF]. Ascribing the observed losses in such cases to diffusion only is bound to overestimate 40 Ar diffusion rates, but to a largely unknown extent if the crystal-defective structure of the mica is not properly quantified [START_REF] Kramar | Heterogeneous 40 Ar* distributions in naturally deformed muscovite: in situ UV-laser ablation evidence for microstructurally controlled intragrain diffusion[END_REF][START_REF] Kramar | Stacking fault-enhanced argon diffusion in naturally deformed muscovite[END_REF]. A case in point are the diffusion rates empirically derived from the muscovites studied by [START_REF] Kirschner | Staircase 40 Ar/ 39 Ar spectra of fine-grained white mica: Timing and duration of deformation and empirical constraints on argon diffusion[END_REF]. These crystallized during deformation by thrusting and most likely contain microstructural defects as well.

The same caveat applies to diffusion experiments performed under hydrothermal conditions [START_REF] Dunlap | Neocrystallization or cooling? 40 Ar/ 39 Ar ages of white micas from low-grade mylonites[END_REF]. To optimize the fractional loss of their starting material, [START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF] used fine-grained muscovite fractions (< 100 µm across) obtained by mechanical grinding of larger grains. Such grinding is likely to have produced extended intracrystalline defects in their samples, a problem already raised by [START_REF] Dunlap | Neocrystallization or cooling? 40 Ar/ 39 Ar ages of white micas from low-grade mylonites[END_REF] for the early thermalhydrothermal experiments conducted by [START_REF] Robbins | Radiogenic argon diffusion in muscovite under hydrothermal conditions[END_REF]. While these defects may have little effect on the effective loss of 40 Ar in vacuum (< 0.5 %; [START_REF] Dunlap | Argon loss during deformation of micas: constraints from laboratory deformation experiments[END_REF], their impact under hydrothermal conditions may be vastly underestimated because of the longer duration of these experiments. Incidentally, [START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF] used a multi-domain diffusion model in an attempt to fit their hydrothermal-loss data to complex in vacuo generated 40 Ar/ 39 Ar spectrum patterns. Such a non-ideal behavior may to some extent reflect the complex defective structure artificially induced during sample preparation.

An additional complication is the presence of reactive fluids, both in Nature and experiments. When out of equilibrium with the solid phase, fluid-phase dissolution-precipitation ture and width, and exchange duration [START_REF] Cole | Rates and mechanisms of isotopic exchange[END_REF][START_REF] Putnis | Mineral replacement reactions[END_REF]. Fluidphase dissolution kinetics are much faster than solid-state intracrystalline diffusion [START_REF] Scaillet | K-Ar ( 40 Ar/ 39 Ar) geochronology of ultrahigh pressure rocks[END_REF][START_REF] Putnis | Mineral replacement reactions[END_REF], potentially enhancing 40 Ar loss well beyond lattice diffusion rates in hydrothermal experiments. [START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF] report a fractional loss 10 % higher for a hydrothermal run compared to the same P -T run under dry conditions. Although the authors infer that a fluid phase is a necessary reservoir for 40 Ar to be diffusively lost from muscovite during the experiments, the difference may also reflect 40 Ar diffusion enhanced in the presence of water. [START_REF] Harrison | Diffusion of 40 Ar in muscovite[END_REF] also report an Al enrichment of their run products compared to the starting material. As discussed by [START_REF] Villa | The in vacuo release of Ar from minerals: 1. Hydrous minerals[END_REF], this may hint at dissolution-precipitation transformations during the runs, potentially compromising (i.e., overestimating) the diffusion rates derived from the observed losses. Evidence of recrystallization of micas during hydrothermal experiments is also reported by [START_REF] Hess | Interpretation of 40 Ar/ 39 Ar spectra of biotites: Evidence from hydrothermal degassing experiments and TEM studies[END_REF] and [START_REF] Villa | Geochronology of the Larderello geothermal field: new data and the "closure temperature" issue[END_REF].

Although not dealing with 40 Ar diffusion in muscovite proper, studies of Sr diffusion in biotite [START_REF] Hammouda | Diffusion of Sr in fluorphlogopite determined by Rutherford backscattering spectrometry[END_REF] have revealed four orders magnitude difference in diffusion rates between hydrothermally treated biotite (faster-diffusing) and dry F-phlogopite (slower-diffusing). These authors rule out F-substitutional effects in favor of the weakening effect of water on the mica structure to explain this behavior, possibly enhanced by dissolution-recrystallization processes during the hydrothermal treatment. Water molecules can also act as source of H + , which weaken Si-O and Al-O bonds and consequently enhance the diffusion of oxygen in minerals such as quartz [START_REF] Elphick | The effect of hydrogen on oxygen diffusion in quartz: evidence for fast proton transients?[END_REF], alkali feldspar [START_REF] Farver | The effect of hydrogen, oxygen, and water fugacity on oxygen diffusion in alkali feldspar[END_REF] and micas [START_REF] Fortier | Volume self-diffusion of oxygen in biotite, muscovite, and phlogopite micas[END_REF]. Hydrolysis of K-O bonds in muscovite, similar to that of Si-O and Al-O in other common silicates, could lower the 40 Ar migration barrier, thereby accelerating its diffusion. All these mechanisms suggests that 40 Ar diffusivity measured under hydrothermal conditions provides only an upper limit to actual 40 Ar transport rates across the mica lattice.

7.2. Applicability to natural systems and implications for 40 Ar/ 39 Ar thermochronology All the above processes work the same in Nature, but with even more free thermodynamic variables (e.g., multiphase equilibria, transient fluids) and much longer durations than in hy-drothermal diffusion cells. In fluid-assisted metamorphic transformations, the 40 Ar initially present in recrystallized (or dissolved) muscovite may be released in the free fluid phase, locally enhancing or disrupting 40 Ar concentration gradients pre-existing in the original crystals (e.g. [START_REF] Hames | On the loss of 40 Ar* from muscovite during polymetamorphism[END_REF][START_REF] Scaillet | K-Ar ( 40 Ar/ 39 Ar) geochronology of ultrahigh pressure rocks[END_REF][START_REF] Vincenzo | Evidence of a close link between petrology and isotope records: constraints from SEM, EMP, TEM and in situ 40 Ar-39 Ar laser analyses on multiple generations of white micas (Lanterman Range, Antarctica)[END_REF][START_REF] Rapaport | The art of molecular dynamics simulation[END_REF][START_REF] Hames | Single-crystal 40 Ar/ 39 Ar age variation in muscovite of the Gassetts Schist and associated gneiss, Vermont Appalachians[END_REF][START_REF] Beltrando | Preservation of sub-microscopic structural relicts in micas from the Gran Paradiso Massif (Western Alps): Implications for 40 Ar-39 Ar geochronology[END_REF][START_REF] Beaudoin | Situ and Step-Heating 40 Ar/ 39 Ar Dating of White Mica in Low-Temperature Shear Zones (Tenda Massif[END_REF][START_REF] Laurent | 40 Ar behaviour and exhumation dynamics in a subduction channel from multi-scale 40 Ar/ 39 Ar systematics in phengite[END_REF]. Partial to complete 40 Ar resetting with intracrystalline 40 Ar/ 39 Ar gradients do exist however, as documented in situ in natural muscovite (e.g., [START_REF] Scaillet | 40 Ar/ 39 Ar laser-probe dating by step heating and spot fusion of phengites from the Dora Maira nappe of the Western Alps, Italy[END_REF][START_REF] Scaillet | Mg/Fe and [(Mg,Fe)Si-Al 2 ] compositional control on argon behaviour in high-pressure white micas: A 40 Ar/ 39 Ar continuous laser-probe study from the Dora-Maira nappe of the internal Western Alps, Italy[END_REF][START_REF] Hames | Laser 40 Ar/ 39 Ar evaluation of slow cooling and episodic loss of 40 Ar from a sample of polymetamorphic muscovite[END_REF][START_REF] Hodges | 40 Ar/ 39 Ar age gradients in micas from a high-temperature-low-pressure metamorphic terrain: Evidence for very slow cooling and implications for the interpretation of age spectra[END_REF][START_REF] Hodges | 40 Ar/ 39 Ar thermochronology of isotopically zoned micas: Insights from the southwestern USA proterozoic orogen[END_REF][START_REF] Scaillet | Excess 40 Ar transport scale and mechanism in high-pressure phengites: A case study from an eclogitized metabasite of the Dora-Maira nappe, western Alps[END_REF][START_REF] Boundy | Fine-scale isotopic heterogeneities and fluids in the deep crust: a 40 Ar/ 39 Ar laser ablation and TEM study of muscovites from a granulite-eclogite transition zone[END_REF].

The mere existence of these spatial gradients is indication that bulk (volume) diffusive motion is effectively operating in Nature, but not with the rate constants and barrier potentials theoretically predicted for ideal (undefective) muscovite. The empirical Ar-Rb-Sr parallel drawn above and the highly disparate 40 Ar diffusion rates documented both experimentally and in the field clearly indicate that 40 Ar/ 39 Ar systematics in natural systems are not fully captured yet in terms of lattice-scale processes.

Our findings point to a fundamental, yet largely unexplored, crystal-structural thread connecting the Ar isotopic record with mixed-scale diffusional ± compositional ± structural disequilibrium interactions during the natural "growth → closure" cycle of muscovite.

During metamorphic P -T changes, white micas react via segregation and re-organization of Si 4+ -rich domains involving the cooperative interdiffusion of several components such as 2Na + + Al 3+ = K + + Si 4+ + □ that progressively produce discrete compositional domains during re-equilibration [START_REF] Shau | Genesis and solvus relations of submicroscopically intergrown paragonite and phengite in a blueschist from northern California[END_REF], analogous to segregation of interlayer cations (Na + , Ca 2+ , K + ) into nanodomains [START_REF] Livi | White mica domain formation: A model for paragonite, margarite, and muscovite formation during prograde metamorphism[END_REF][START_REF] Beltrando | Preservation of sub-microscopic structural relicts in micas from the Gran Paradiso Massif (Western Alps): Implications for 40 Ar-39 Ar geochronology[END_REF]. Crystal-structural evidence above suggests that such modifications are decoupled kinetically across the interlayer hosting radiogenic 40 Ar, especially in the presence of grain-scale chemical potential and stress gradients [START_REF] Scaillet | Mg/Fe and [(Mg,Fe)Si-Al 2 ] compositional control on argon behaviour in high-pressure white micas: A 40 Ar/ 39 Ar continuous laser-probe study from the Dora-Maira nappe of the internal Western Alps, Italy[END_REF][START_REF] Scaillet | Excess 40 Ar transport scale and mechanism in high-pressure phengites: A case study from an eclogitized metabasite of the Dora-Maira nappe, western Alps[END_REF][START_REF] Scaillet | K-Ar ( 40 Ar/ 39 Ar) geochronology of ultrahigh pressure rocks[END_REF], severely impairing the interpretation of microscale 40 Ar/ 39 Ar patterns in terms of either diffusion or (re)crystallization [START_REF] Beaudoin | Situ and Step-Heating 40 Ar/ 39 Ar Dating of White Mica in Low-Temperature Shear Zones (Tenda Massif[END_REF][START_REF] Laurent | 40 Ar behaviour and exhumation dynamics in a subduction channel from multi-scale 40 Ar/ 39 Ar systematics in phengite[END_REF].

Hence, the reason why empirical 40 Ar diffusion rates are so distinct from DFT and MD estimates is not failure of these approaches in describing fundamental lattice-based processes.

Rather, it is the manifestation of atomic-scale kinetics combining diffusion with (as yet) largely unmodeled physics, most notably those controlling reactivity of the interlayer with the external medium. Thermodynamic equilibrium is the common-denominator assumption of both experimental and theoretical diffusion modeling. So far, fulfilling this condition has proven elusive in the hydrothermal treatment of muscovite [START_REF] Villa | The in vacuo release of Ar from minerals: 1. Hydrous minerals[END_REF]. In contrast, equilibrium is a built-in condition in theoretical models ensuring, as here, full control on boundary conditions along with a test that the system indeed remains thermodynamically stable. As such, our MD and TST diffusion estimates provide the first reliable benchmark against which to gauge first-order 40 Ar diffusion kinetics in a stoichiometrically stable muscovite structure.

Combined with hydrothermal predictions (Fig. 6), they collectively define the probable D-T space spanned by intracrystalline 40 Ar kinetics in natural muscovite. Further empirical and theoretical work is required to clarify kinetic interactions between crystal-chemical disequilibrium and effective 40 Ar diffusion across such D-T space.

Conclusions

A lattice-based study of 40 Ar diffusion in muscovite has been attempted for the first time via Molecular Dynamics simulations coupled with Nudged Elastic Band and Transition State

Theory. The main results lead to the following conclusions:

• 40 Ar diffusion in 2M 1 muscovite predominantly occurs along the (001) interlayer with recoil energies too low to permit extensive or systematic implantation of the radiogenic 40 Ar atom more than a lattice spacing away from its production site from 40 K.

• Systematic energy-barrier mapping by NEB coupled with TST allows to identify the most probable diffusion pathway and barrier potential for 40 Ar motion along the (001) interlayer;

40 Ar transport via a divacancy mechanism is shown to be the more energetically favorable pathway for 40 Ar migration along (001) with characteristic diffusion parameters E = 66 kcal.mol -1 and D 0 = 6 × 10 -4 cm 2 .s -1 .

• These theoretical estimates suggest that purely diffusion-driven (Brownian) motion is not an efficient mechanism for 40 Ar intracrystalline transport in a stoichiometrically pure and stable muscovite structure.

• Along with hydrothermal and field-based evidence, they call for more complex (multi-scale and multi-process) physics to explain 40 Ar exchange and retention properties of natural muscovite, most notably by considering the chemical-stoichiometric reactivity of the interlayer with the external medium.

Figure captions

40

  Ar had the same initial velocity (0.11765 Å.fs -1 ) but varying recoil angles to cover all possible values within a hemisphere. They identified four possible final positions of 40 Ar. For recoil vectors subparallel to the (001) plane, 40 Ar systematically remained in the interlayer region, occupying either the initial site of the parent 40 K or a neighboring vacancy. For recoil vectors pointing towards the T-O-T sheet, 40 Ar remained mainly in the interlayer region (79.5 %) but also penetrated into the T-O-T layer (20.5 %), where it ended up either in the tetrahedral layer adjacent to the interlayer region (16.5 %), or in the octahedral layer (1.5 %) or in the opposite tetrahedral layer (2.5 %).

  energies towards monovacancies show a mean value of 79 ± 9 kcal.mol -1 (1σ). Using DFT simulations,[START_REF] Yu | First-principles study of organically modified muscovite mica with ammonium (NH + 4 ) or methylammonium (CH 3 NH + 3 )[END_REF] computed energy barriers associated with the migration of K + towards a single vacancy in muscovite, for atomic models differing in their number of tetrahedral Al surrounding the interlayer region where the migrating potassium is located.They report energy barriers of (i) 82 kcal.mol -1 for the Al-Al model which contains Al atoms in the two tetrahedral layers adjacent to the interlayer region; (ii) 146 kcal.mol -1 for the Al-Si model in which Al atoms are present only in one of the tetrahedral layers, and (iii) 227 kcal.mol -1 for the Si-Si model with tetrahedral layer adjacent to the considered interlayer layer region fully occupied by Si atoms. Noteworthy, the energy barrier reported for their Al-Al model, which is equivalent to ours, is within the range of values computed in our study with ClayFF. Ruiz Pestana et al. (2017) report a distribution of K + energy barriers in illite computed via ClayFF and DFT. They get mean values of 54 ± 12 kcal.mol -1 forClayFF energy barriers and 75 ± 8 kcal.mol -1 for those computed using DFT. They conclude that ClayFF underestimates energy barriers due to the non-bonded nature of the force field which makes structures excessively deformable. However, this apparent discordance between ClayFF and DFT is probably related to the fact that Ruiz[START_REF] Pestana | Direct exchange mechanism for interlayer ions in non-swelling clays[END_REF] performed their ClayFF simulations without prior equilibration of their illite structure. In our study, the K + energy barriers computed using ClayFF after prior equilibration of the muscovite structure are in the range of their DFT-computed barriers.

40

  Ar migration energy barriers E Ar □ eff = 69 kcal.mol -1 and E Ar □□ eff = 65 kcal.mol -1 for the monovacancy and the divacancy mechanism, respectively. For K + (or vacancy) migration, the calculated effective barriers are E K □ eff = 75 kcal.mol -1 (monovacancy mechanism) andE K □□ eff= 58 kcal.mol -1 (divacancy mechanism).6.1. Monovacancy mechanism 6.1.1. Diffusion of monovacancies in a hexagonal lattice of K +The diffusion of 40 Ar via the monovacancy mechanism requires the presence of a vacancy in one of its neighboring sites. As long as no vacancies are available, 40 Ar cannot move, requiring first to study the diffusivity of the vacancies in the interlayer region. When a

Fig. 1 .

 1 Fig. 1. The crystal structure atomic model of muscovite imaged with VESTA (Momma and Izumi, 2011). The system (n = 1512 atoms) consists of a supercell (6 × 3 × 1 unit cells) of initial dimensions L x = 31.146 Å, L y = 27.018 Å and L z = 19.966 Å. a) Perspective view down the a axis. b) Tetrahedral layers (as labelled in (a)) showing the irregular distribution of Al atoms. Atoms are drawn with their Van der Waals radii. T up and T low indicate the upper and the lower tetrahedral layers with respect to the interlayer region, respectively.

Fig. 2 .

 2 Fig.2. Variation of the energy of the system during 40 Ar recoil with respect to the mean energy of a system in which 40 Ar occupies an regular interlayer site. The red line is the running mean energy averaged over 2.5 ps. The decrease of the energy at 50 ps corresponds to the move of 40 Ar from the interstitial site to the lattice stable site.

Fig. 3 .

 3 Fig. 3. Geometry of diffusion in the muscovite interlayer plane. a) Monovacancy mechanism. b) Divacancy mechanism. Orange and green circles correspond to K + and 40 Ar atoms respectively. □ represents the vacancy, -→ V i are the displacement vectors and -→ δ i indicate atomic jumps preserving the divacancy. See text for discussion.

  Fig. 5. Results of NEB calculations of K + migration in muscovite. a) Histogram of migration barriers associated with a monovacancy mechanism. b) Histogram of migration barri-

Figures 1177 Fig.

 1177 Figures 1177

  Fig. 1

  Fig.2

Table 1 :

 1 Interaction parameters used in the simulations

	Nonbonded parameters
	Atoms	Charge (e)	σ ( Å)	ϵ (kcal.mol -1 )
	H (hydroxyl) a	0.4250	-	-
	O1 (hydroxyl) a	-0.9500	3.1655	0.1554
	O2 (bridging Oxygen) a	-1.0500	3.1655	0.1554
	O3 (Oxygen bonded to tetrahedral Al) a	-1.1688	3.1655	0.1554
	Al (octahedral) a	1.5750	4.2712	1.3298 × 10 -6
	Al2 (tetrahedral) a	1.5750	3.3020	1.8405 × 10 -6
	Si a	2.1000	3.3020	1.8405 × 10 -6
	K a	1.0000	3.3340	0.1000
	Ar b	0.0000	3.4010	0.2304
	Bonded parameters
	Atom 1	Atom 2	k 1 (kcal.mol -1 . Å-2 )	r 0 ( Å)
	O1 (hydroxyl)	H (hydroxyl)	554.1349	1.0000

a

[START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF]

b

[START_REF] White | Lennard-Jones as a model for argon and test of extended renormalization group calculations[END_REF]

Table 2 :

 2 Comparison of experimental and simulated crystallographic parameters of 2M 1 Muscovite

		Experiment MD simulations (CLAYFF) Deviation (± %)
	a ( Å)	5.191	5.182	-0.17
	b ( Å)	9.006	8.944	-0.69
	c ( Å)	20.068	19.697	-1.85
	α ( • )	90.000	89.787	-0.24
	β ( • )	95.770	95.936	+0.17
	γ ( • )	90.000	90.175	+0.19
	Interlayer sep. ( Å)	3.396	3.368	-0.88

Table 3 :

 3 Computed migration barriers and 40 Ar diffusion parameters for the monovacancy (□) and divacancy (□□) mechanisms.

		Mean energy barriers	Effective energy barriers	40	Ar Diffusion parameters
		⟨E Ar ⟩	⟨E K ⟩	E Ar eff	E K eff		E	D 0
		(kcal.mol -1 ) (kcal.mol -1 )	(kcal.mol -1 ) (kcal.mol -1 )	(kcal.mol -1 ) (cm 2 .s -1 )
	□	71 ± 5	79 ± 9	69	75		75	4 × 10 -3
	□□	69 ± 6	66 ± 10	65	58		66	6 × 10 -4

Acknowledgements

Work supported by LABEX grant VOLTAIRE (ANR-10-LABX-100-01), the Région Centre grant ARGON, and the EQUIPEX grant PLANEX (ANR-11-EQPX-0036). J. Nteme is supported by a post-doctoral grant from the project LABEX-VOLTAIRE. Jan Wijbrans (VU, Amsterdam) and two anonymous reviewers are thanked for providing insightful comments and references that helped to improve the content and clarity of the work. Associate editor F. Jourdan is thanked for efficient handling of the manuscript as well as for a flash review. Parts of the material and ideas developed in this paper form the basis of the project "DeconvAr" discarded as "unlikely to lead to a breakthrough and widespread application of the findings" under the 2021 AAPG call "Planétologie, structure et histoire de la Terre".

Appendix A. Supplementary Material

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/xx.xxx/xxxxxxxxxxxx. along its path towards a monovacancy (black) and towards a vacancy pair (blue). d) Trajectories followed by K + as it migrates from a stable site towards a single vacancy (black) and towards a vacancy pair (blue). Red circles represent basal oxygens. e)

Diagram showing the positive correlation between the 40 Ar and K + migration barriers. ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: