# Effect of heterogeneity on the diffusion of Pb in apatite for petrochronological applications: A multiscale approach to characterizing the influence of apatite chemistry and anisotropy on Pb diffusion

Matee Ullah, Urs Klötzli, Cecile Gautheron, Laurent Tassan-Got, Fakhrul Islam, Muhammad Younas, Kamran Shehzad, Muhammad Khubab, Muhammad Ibrar, Bilal Wadood

► To cite this version:

Matee Ullah, Urs Klötzli, Cecile Gautheron, Laurent Tassan-Got, Fakhrul Islam, et al. Effect of heterogeneity on the diffusion of Pb in apatite for petrochronological applications: A multiscale approach to characterizing the influence of apatite chemistry and anisotropy on Pb diffusion. Lithos, 2023, 460-461, pp.107396. 10.1016/j.lithos.2023.107396. hal-04337673

HAL Id: hal-04337673

https://hal.science/hal-04337673

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

Check for
updates

# Effect of heterogeneity on the diffusion of Pb in apatite for petrochronological applications: A multiscale approach to characterizing the influence of apatite chemistry and anisotropy on Pb diffusion

Matee Ullah<sup>a,b,*</sup>, Urs Klötzli<sup>a</sup>, Ce’cile Gautheron<sup>c</sup>, Laurent Tassan-Got<sup>d</sup>, Fakhrul Islam<sup>b</sup>, Muhammad Younas<sup>f</sup>, Kamran Shehzad<sup>b</sup>, Muhammad Khubab<sup>e</sup>, Muhammad Ibrar<sup>g</sup>, Bilal Wadood<sup>h</sup>

ARTICLE INFO

Keywords:

U-Th/Pb dating

Petrochronology

Apatite

Lead diffusion

Computation modelling

ABSTRACT

The investigation of various factors effecting the Lead (Pb) diffusion in phosphate minerals such as apatite is still challenging in the interpretation of (U—Th)/Pb geochronology. For (U—Th)/Pb system, apatite minerals have closure temperatures in the range of 375 to 600 °C and therefore can be used for the investigation of mid-temperature thermochronological and/or petrochronological questions i.e., the reconstruction of thermal events in Earth’s crust. There is still uncertainty whether Pb diffusion in apatite is characterized by thermally activated volume and/or anisotropic diffusion profiles or is instead impacted by novel growth processes and recrystallization (chemical substitutions). As the apatite structure support extensive compositional variability, including partial or total substitution of both the cationic and anionic sites and forms solid solutions therefore, it necessitates a thorough examination of these effects and anisotropy on Pb diffusivity and (U—Th)/Pb geochronometric system. For this, a multi-scale study is carried out to examine the effects of chemical composition, anisotropy, and growth structure on the diffusion of Pb in order to better understand the behaviour of Pb diffusion in apatite. This study employed computational techniques like Density Functional Theory (DFT) and Transition State Theory (TST) at the atomic level and integrates it with the Kinetic Monte Carlo (KMC) simulations at the macroscopic level. Models of this study shows that Pb diffusion is completely anisotropic along the preferred z-axis or [001] direction and Pb readily escapes faster from Na-substituted apatite when compared to pure F-apatite and Cl-substituted apatite. Because of this anisotropy and chemical substitutions, Pb diffusivity in apatite either increases by opening of diffusion channels or decreases by blocking the diffusion channels depending on the site and type of chemical substitution. Further, in case of blocking effect the Pb diffusion occurs through workaround pathways and approaches towards the isotropic diffusion. For Na-substituted apatite, the impact of Na occupation on anisotropic Pb diffusion is significantly greater while in case of Cl-substituted apatite the Cl occupation mostly leads towards isotropic diffusion by opening the diffusion paths along other directions (mostly along the in-plane direction). Furthermore, the high closure temperatures (Tc) (e.g., ~1370 °C) of the modelled apatites (except the perfect Na-substituted apatite e.g., ~500 °C) of this study when compared to the Tc of Durango apatite obtained experimentally for the effective grain size of 100 µm and cooling rates of 10 °C/Ma indicate that the effective closure temperature dominantly depends on the degree and types of chemical substitutions and play a crucial role for the closure or opening of Pb diffusion/loss in apatites.

<footnote><sup>a</sup> Department of Lithospheric Research, University of Vienna, Vienna, Austria</footnote>

<footnote><sup>b</sup> Department of Geology, Khushal Khan Khattak University, Karak, Pakistan</footnote>

<footnote><sup>c</sup> GEOPS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, Rue du Belve'de` re, Bat. 504, F-91405 Orsay, France</footnote>

<footnote><sup>d</sup> Institut de Physique Nucle'aire d'Orsay, CNRS-IN2P3, Univ. Paris-Sud, Universite' Paris-Saclay, 91406 Orsay Cedex, France</footnote>

<footnote><sup>e</sup> Department of Earth and Atmospheric Sciences, University of Houston, Houston, USA</footnote>

<footnote><sup>f</sup> Department of Geology, University of Vienna, Vienna, Austria</footnote>

<footnote><sup>g</sup> Center for Earth and Space Sciences, University of Swat, Pakistan</footnote>

<footnote><sup>h</sup> Department of Geology, University of Swabi, Anbar-Swabi, 23561, KPK, Pakistan</footnote>

<footnote>* Corresponding author at: Department of Lithospheric Research, University of Vienna, Vienna, Austria.</footnote>

<footnote>E-mail address: mateeu89@univie.ac.at (M. Ullah).</footnote>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Check for
updates

# Effect of heterogeneity on the diffusion of Pb in apatite for petrochronological applications: A multiscale approach to characterizing the influence of apatite chemistry and anisotropy on Pb diffusion

Matee Ullah<sup>a,b,*</sup>, Urs Klötzli<sup>a</sup>, Ce’cile Gautheron<sup>c</sup>, Laurent Tassan-Got<sup>d</sup>, Fakhrul Islam<sup>b</sup>, Muhammad Younas<sup>f</sup>, Kamran Shehzad<sup>b</sup>, Muhammad Khubab<sup>e</sup>, Muhammad Ibrar<sup>g</sup>, Bilal Wadood<sup>h</sup>

ARTICLE INFO

Keywords:

U-Th/Pb dating

Petrochronology

Apatite

Lead diffusion

Computation modelling

ABSTRACT

The investigation of various factors effecting the Lead (Pb) diffusion in phosphate minerals such as apatite is still challenging in the interpretation of (U—Th)/Pb geochronology. For (U—Th)/Pb system, apatite minerals have closure temperatures in the range of 375 to 600 °C and therefore can be used for the investigation of mid-temperature thermochronological and/or petrochronological questions i.e., the reconstruction of thermal events in Earth’s crust. There is still uncertainty whether Pb diffusion in apatite is characterized by thermally activated volume and/or anisotropic diffusion profiles or is instead impacted by novel growth processes and recrystallization (chemical substitutions). As the apatite structure support extensive compositional variability, including partial or total substitution of both the cationic and anionic sites and forms solid solutions therefore, it necessitates a thorough examination of these effects and anisotropy on Pb diffusivity and (U—Th)/Pb geochronometric system. For this, a multi-scale study is carried out to examine the effects of chemical composition, anisotropy, and growth structure on the diffusion of Pb in order to better understand the behaviour of Pb diffusion in apatite. This study employed computational techniques like Density Functional Theory (DFT) and Transition State Theory (TST) at the atomic level and integrates it with the Kinetic Monte Carlo (KMC) simulations at the macroscopic level. Models of this study shows that Pb diffusion is completely anisotropic along the preferred z-axis or [001] direction and Pb readily escapes faster from Na-substituted apatite when compared to pure F-apatite and Cl-substituted apatite. Because of this anisotropy and chemical substitutions, Pb diffusivity in apatite either increases by opening of diffusion channels or decreases by blocking the diffusion channels depending on the site and type of chemical substitution. Further, in case of blocking effect the Pb diffusion occurs through workaround pathways and approaches towards the isotropic diffusion. For Na-substituted apatite, the impact of Na occupation on anisotropic Pb diffusion is significantly greater while in case of Cl-substituted apatite the Cl occupation mostly leads towards isotropic diffusion by opening the diffusion paths along other directions (mostly along the in-plane direction). Furthermore, the high closure temperatures (Tc) (e.g., ~1370 °C) of the modelled apatites (except the perfect Na-substituted apatite e.g., ~500 °C) of this study when compared to the Tc of Durango apatite obtained experimentally for the effective grain size of 100 µm and cooling rates of 10 °C/Ma indicate that the effective closure temperature dominantly depends on the degree and types of chemical substitutions and play a crucial role for the closure or opening of Pb diffusion/loss in apatites.

<footnote><sup>a</sup> Department of Lithospheric Research, University of Vienna, Vienna, Austria</footnote>

<footnote><sup>b</sup> Department of Geology, Khushal Khan Khattak University, Karak, Pakistan</footnote>

<footnote><sup>c</sup> GEOPS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, Rue du Belve'de` re, Bat. 504, F-91405 Orsay, France</footnote>

<footnote><sup>d</sup> Institut de Physique Nucle'aire d'Orsay, CNRS-IN2P3, Univ. Paris-Sud, Universite' Paris-Saclay, 91406 Orsay Cedex, France</footnote>

<footnote><sup>e</sup> Department of Earth and Atmospheric Sciences, University of Houston, Houston, USA</footnote>

<footnote><sup>f</sup> Department of Geology, University of Vienna, Vienna, Austria</footnote>

<footnote><sup>g</sup> Center for Earth and Space Sciences, University of Swat, Pakistan</footnote>

<footnote><sup>h</sup> Department of Geology, University of Swabi, Anbar-Swabi, 23561, KPK, Pakistan</footnote>

<footnote>* Corresponding author at: Department of Lithospheric Research, University of Vienna, Vienna, Austria.</footnote>

<footnote>E-mail address: mateeu89@univie.ac.at (M. Ullah).</footnote>
# 1. Introduction

The thermal histories of many common and accessory minerals found on Earth, on other planets, and in meteorites can reveal important details about the timing of planetary accretion, plate tectonics, mountain creation, impact events, and sedimentary basins evolutionary histories (Kirkland et al., 2017, 2018; Spikings et al., 2010; Webster and Piccoli, 2015). Thermochronological techniques such as (U-Th-Sm)/He, $^{40}$Ar/$^{39}$Ar (alkali feldspars), fission track, and $^{238,235 U}$-206,207Pb (accessory phases, such as apatite, titanite, and rutile), can be used to compel the thermal histories of common and/or accessory minerals (e. g., Cherniak et al., 2004; Watson and Baxter, 2007; Kirkland et al., 2017; Zhang et al., 2021). These techniques are based on the fundamental concept that radiogenic isotopes (such as $^{206,207Pb}$, $^{40}$Ar, and $^{4}$He) produced are redistributed within minerals before being finally lost to some infinite reservoirs through thermally triggered volume diffusion and other diffusion mechanisms (Blackburn et al., 2011; Cherniak, 2010; Cherniak et al., 2004; Schoene and Bowring, 2007).

A variety of minerals, including zircon, monazite, apatite, and iron oxides, have been employed for (U—Th)/Pb, He thermochronology in numerous geological applications (e.g., Barnes et al., 2021; Farley and Stockli, 2002; Schoene, 2014; Watson et al., 1985; Zeitler et al., 1987). Because of its high temperature sensitivity, the U—Pb method is well suited for studying tectonic stability, exhumation, the geologic history of tectonically active margins that may have been subjected to multiple terrane collisions, extra-terrestrial phenomena such as asteroid impact and ejection events, and rates of planetary accretion (Blackburn et al., 2011; Blackburn et al., 2012; McGregor et al., 2020). Strong quantitative understanding of the diffusion behaviour in geochronometric minerals is a prerequisite for a robust U-Th-Pb age interpretation as the closure temperature $T_c$ (which denotes the temperature below which Pb are largely preserved during the process of cooling) of these minerals is a function of the diffusion behaviour (Dodson, 1973). Moreover, the diffusion behaviour of apatite and zircon has been demonstrated to be quite variable. However, theoretically, both the isotropic and anisotropic diffusion of Pb in both these minerals obey Arrhenius law (Watson et al., 1985; Cherniak et al., 1991; Cherniak and Watson, 2001; Kirkland et al., 2018).

The calcium phosphate mineral apatite, which has the general formula $Ca_5(PO_4)_3(F,Cl,OH)$, is common in magmatic as well as sedimentary rocks and can also grow and be altered during metamorphic processes. Apatite mineral crystalizes in hexagonal structure as being composed of hexagonally arranged Ca(1) and Ca(2) coordination polyhedra around a central [001] hexad (Hughes et al., 1989). The F ions are arranged in columns on the hexads along the z-axis at (0, 0, z) positions. There are three longer Ca(1)-O(1) bonds and six shorter Ca(1)-O(1) bonds that form the nine-coordinate Ca(1) polyhedron, which approximately defines a trigonal prism geometry. The Ca(1)O9 polyhedra formed chains parallel to z-axis by sharing their (001) pinacoid faces. On the other hand, the Ca(2) polyhedron has an octahedron (CaO5X) type geometry having a seventh weak O-bond of about 0.15 valence units in fluorapatite, which redefines the polyhedron as a CaO5X(O) (Hughes et al., 1989). The triangles formed by the Ca(2) atoms in mirror planes centres around [001] zone axes consists of the X anions (e.g. in fluorapatite the F anions; Hughes et al., 1989). Due to its crustal abundance and capacity to host U, apatite has found use in wide range of petrochronological and thermochronological problems (Cherniak et al., 1991; Chamberlain and Bowring, 2001; Schoene and Bowring, 2007; Yang et al., 2014; Chew and Spikings, 2015; Mark et al., 2016; Kirkland et al., 2017; Wohlgemuth-Ueberwasser et al., 2017). Unlike zircon, apatite is often found in silica-undersaturated mafic lithologies. It typically contains moderate to substantial levels of non-radiogenic/common Pb, and its U contents are frequently lower than that of zircon. Apatite U-Th/Pb geochronology faces challenges due to their large ratio of common to radiogenic Pb, and computed dates rarely match the accuracy of zircon dates (Chew et al., 2011; Chew and Spikings, 2015; Glorie et al., 2017;

Kirkland et al., 2017; Krogstad and Walker, 1994; LaFlamme et al., 2017; Lu et al., 2021). However, its typical closure temperature (which is between 375 °C and 600 °C) gives an additional time-temperature constraint on intermediate temperature prograde or retrograde processes and therefore, apatite is still considered as a powerful U-Th/Pb chronometer (Cochrane et al., 2014; Henrichs et al., 2018; Kirkland et al., 2017; Kirkland et al., 2018; Schoene and Bowring, 2007).

Radiation damage is thought to be a major contribution to the variation in diffusion, although the effects of damage type (recoil, alpha, and fission damage), topology (large damage, Frenkel-pair, and amorphous zone), and even variation in crystal's chemistry on Pb diffusion in apatite are not fully understood (Bengtson et al., 2012; Djimbi et al., 2015; Farley, 2000; Flowers et al., 2007). Other complexities including inclusions, voids, and defects could also cause comparable issues (Bassal et al., 2022; Danisík et al., 2017; Tang et al., 2017; Zeitler et al., 2017). The amount of progress made on poorly understood natural apatites and even other natural minerals by using only empirical observations and experiments is constrained due to these complexities (e.g., Cherniak et al., 1991; Reiners et al., 2004; Cherniak et al., 2004; Cherniak, 2010; Chew et al., 2011; Chew and Spikings, 2015). So, a thorough understanding of the physics behind diffusion process and how defects and change in chemistry of the minerals affect diffusion behaviour is crucial addition to the various and comprehensive data collection currently being available. However, no study has been conducted about the effect of the variable chemistry on the diffusion of Pb in apatite. The actual chemistry of these chronometric minerals is highly variable from sample to sample, although a general insight to such chemical changes and its effect on the transport/diffusion behaviour of Pb in apatite can be investigated via theoretical and computational approaches. There have been several studies that have used the DFT to estimate the quantitative diffusion of atoms in minerals (Balout et al., 2017a, 2017b; Bassal et al., 2022; Bengtson et al., 2012; Djimbi et al., 2015; Gautheron et al., 2020; Reich et al., 2007; Saadoune et al., 2009; Ullah et al., 2023).

This study also contributes to these sorts of studies to investigate the Pb diffusion in apatite at the atomistic level using a package of methods adopted in recent studies of Djimbi et al. (2015), Balout et al. (2017a, 2017b), Gerin et al. (2017), Gautheron et al. (2020), Bassal et al. (2022), and Ullah et al. (2023) to explore some of the vital factors that can affect the diffusion rates. Because Pb diffuses much faster than U and Th, the closure temperature of Pb, which is largely dependent on Pb diffusivity, is of greater interest in a U-Th/Pb geochronometric system than U and Th (Cherniak, 2005; Cherniak, 2010). To investigate the effect of apatite chemistry on Pb diffusion, the Durango apatite was chosen for comparison due to the availability of experimental data, and the chemical substitutions of Na (almost negligible with 0.08 atomic proportion per formula unit (apfu) of Ca) and Cl were made accordingly in the perfect apatite supercell used in this study. At first, the DFT is implemented at the atomic level to optimize the ideal apatite super cell and exhibiting for quantum characteristics with inventory insertion and interstitial sites. Following the identification of insertion and interstitial sites, migratory energies between the two adjacent symmetrically equivalent sites are calculated using a combined approach of Nudged Elastic Band (NEB) and DFT method. The jump rates between the two sites are then estimated using Transition State Theory (TST). To account for heterogeneities at higher scale, the results of both DFT and TST are employed in a Kinetic Monte Carlo (KMC) simulation to investigate the 3D diffusion behaviour and estimate the 3D effective diffusion coefficients. The KMC simulations model the 3D obstruction and/or opening of the diffusion pathways caused by change in chemistry of the supercell at mineral scales and estimate their effect on the closure temperature. At the end all the results were compared with the published experiments so that to get the idea about the heterogeneity that can modify Pb diffusion behaviour in apatite which is a prime criterion for a robust (U—Th)/Pb age interpretation, thermal histories, and the closure temperature $T_c$ of apatite as $T_c$ is a function of the diffusion behaviour (Smye et al., 2018).
# 1. Introduction

The thermal histories of many common and accessory minerals found on Earth, on other planets, and in meteorites can reveal important details about the timing of planetary accretion, plate tectonics, mountain creation, impact events, and sedimentary basins evolutionary histories (Kirkland et al., 2017, 2018; Spikings et al., 2010; Webster and Piccoli, 2015). Thermochronological techniques such as (U-Th-Sm)/He, $^{40}$Ar/$^{39}$Ar (alkali feldspars), fission track, and $^{238,235 U-206,207}$Pb (accessory phases, such as apatite, titanite, and rutile), can be used to compel the thermal histories of common and/or accessory minerals (e. g., Cherniak et al., 2004; Watson and Baxter, 2007; Kirkland et al., 2017; Zhang et al., 2021). These techniques are based on the fundamental concept that radiogenic isotopes (such as $^{206,207}$Pb, $^{40}$Ar, and $^{4}$He) produced are redistributed within minerals before being finally lost to some infinite reservoirs through thermally triggered volume diffusion and other diffusion mechanisms (Blackburn et al., 2011; Cherniak, 2010; Cherniak et al., 2004; Schoene and Bowring, 2007).

A variety of minerals, including zircon, monazite, apatite, and iron oxides, have been employed for (U—Th)/Pb, He thermochronology in numerous geological applications (e.g., Barnes et al., 2021; Farley and Stockli, 2002; Schoene, 2014; Watson et al., 1985; Zeitler et al., 1987). Because of its high temperature sensitivity, the U—Pb method is well suited for studying tectonic stability, exhumation, the geologic history of tectonically active margins that may have been subjected to multiple terrane collisions, extra-terrestrial phenomena such as asteroid impact and ejection events, and rates of planetary accretion (Blackburn et al., 2011; Blackburn et al., 2012; McGregor et al., 2020). Strong quantitative understanding of the diffusion behaviour in geochronometric minerals is a prerequisite for a robust U-Th-Pb age interpretation as the closure temperature $T_c$ (which denotes the temperature below which Pb are largely preserved during the process of cooling) of these minerals is a function of the diffusion behaviour (Dodson, 1973). Moreover, the diffusion behaviour of apatite and zircon has been demonstrated to be quite variable. However, theoretically, both the isotropic and anisotropic diffusion of Pb in both these minerals obey Arrhenius law (Watson et al., 1985; Cherniak et al., 1991; Cherniak and Watson, 2001; Kirkland et al., 2018).

The calcium phosphate mineral apatite, which has the general formula $Ca_5(PO_4)_3(F,Cl,OH)$, is common in magmatic as well as sedimentary rocks and can also grow and be altered during metamorphic processes. Apatite mineral crystalizes in hexagonal structure as being composed of hexagonally arranged Ca(1) and Ca(2) coordination polyhedra around a central [001] hexad (Hughes et al., 1989). The F ions are arranged in columns on the hexads along the z-axis at (0, 0, z) positions. There are three longer Ca(1)-O(1) bonds and six shorter Ca(1)-O(1) bonds that form the nine-coordinate Ca(1) polyhedron, which approximately defines a trigonal prism geometry. The Ca(1)O9 polyhedra formed chains parallel to z-axis by sharing their (001) pinacoid faces. On the other hand, the Ca(2) polyhedron has an octahedron (CaO5X) type geometry having a seventh weak O-bond of about 0.15 valence units in fluorapatite, which redefines the polyhedron as a CaO5X(O) (Hughes et al., 1989). The triangles formed by the Ca(2) atoms in mirror planes centres around [001] zone axes consists of the X anions (e.g. in fluorapatite the F anions; Hughes et al., 1989). Due to its crustal abundance and capacity to host U, apatite has found use in wide range of petrochronological and thermochronological problems (Cherniak et al., 1991; Chamberlain and Bowring, 2001; Schoene and Bowring, 2007; Yang et al., 2014; Chew and Spikings, 2015; Mark et al., 2016; Kirkland et al., 2017; Wohlgemuth-Ueberwasser et al., 2017). Unlike zircon, apatite is often found in silica-undersaturated mafic lithologies. It typically contains moderate to substantial levels of non-radiogenic/common Pb, and its U contents are frequently lower than that of zircon. Apatite U-Th/Pb geochronology faces challenges due to their large ratio of common to radiogenic Pb, and computed dates rarely match the accuracy of zircon dates (Chew et al., 2011; Chew and Spikings, 2015; Glorie et al., 2017;

Kirkland et al., 2017; Krogstad and Walker, 1994; LaFlamme et al., 2017; Lu et al., 2021). However, its typical closure temperature (which is between 375 °C and 600 °C) gives an additional time-temperature constraint on intermediate temperature prograde or retrograde processes and therefore, apatite is still considered as a powerful U-Th/Pb chronometer (Cochrane et al., 2014; Henrichs et al., 2018; Kirkland et al., 2017; Kirkland et al., 2018; Schoene and Bowring, 2007).

Radiation damage is thought to be a major contribution to the variation in diffusion, although the effects of damage type (recoil, alpha, and fission damage), topology (large damage, Frenkel-pair, and amorphous zone), and even variation in crystal's chemistry on Pb diffusion in apatite are not fully understood (Bengtson et al., 2012; Djimbi et al., 2015; Farley, 2000; Flowers et al., 2007). Other complexities including inclusions, voids, and defects could also cause comparable issues (Bassal et al., 2022; Danisík et al., 2017; Tang et al., 2017; Zeitler et al., 2017). The amount of progress made on poorly understood natural apatites and even other natural minerals by using only empirical observations and experiments is constrained due to these complexities (e.g., Cherniak et al., 1991; Reiners et al., 2004; Cherniak et al., 2004; Cherniak, 2010; Chew et al., 2011; Chew and Spikings, 2015). So, a thorough understanding of the physics behind diffusion process and how defects and change in chemistry of the minerals affect diffusion behaviour is crucial addition to the various and comprehensive data collection currently being available. However, no study has been conducted about the effect of the variable chemistry on the diffusion of Pb in apatite. The actual chemistry of these chronometric minerals is highly variable from sample to sample, although a general insight to such chemical changes and its effect on the transport/diffusion behaviour of Pb in apatite can be investigated via theoretical and computational approaches. There have been several studies that have used the DFT to estimate the quantitative diffusion of atoms in minerals (Balout et al., 2017a, 2017b; Bassal et al., 2022; Bengtson et al., 2012; Djimbi et al., 2015; Gautheron et al., 2020; Reich et al., 2007; Saadoune et al., 2009; Ullah et al., 2023).

This study also contributes to these sorts of studies to investigate the Pb diffusion in apatite at the atomistic level using a package of methods adopted in recent studies of Djimbi et al. (2015), Balout et al. (2017a, 2017b), Gerin et al. (2017), Gautheron et al. (2020), Bassal et al. (2022), and Ullah et al. (2023) to explore some of the vital factors that can affect the diffusion rates. Because Pb diffuses much faster than U and Th, the closure temperature of Pb, which is largely dependent on Pb diffusivity, is of greater interest in a U-Th/Pb geochronometric system than U and Th (Cherniak, 2005; Cherniak, 2010). To investigate the effect of apatite chemistry on Pb diffusion, the Durango apatite was chosen for comparison due to the availability of experimental data, and the chemical substitutions of Na (almost negligible with 0.08 atomic proportion per formula unit (apfu) of Ca) and Cl were made accordingly in the perfect apatite supercell used in this study. At first, the DFT is implemented at the atomic level to optimize the ideal apatite super cell and exhibiting for quantum characteristics with inventory insertion and interstitial sites. Following the identification of insertion and interstitial sites, migratory energies between the two adjacent symmetrically equivalent sites are calculated using a combined approach of Nudged Elastic Band (NEB) and DFT method. The jump rates between the two sites are then estimated using Transition State Theory (TST). To account for heterogeneities at higher scale, the results of both DFT and TST are employed in a Kinetic Monte Carlo (KMC) simulation to investigate the 3D diffusion behaviour and estimate the 3D effective diffusion coefficients. The KMC simulations model the 3D obstruction and/or opening of the diffusion pathways caused by change in chemistry of the supercell at mineral scales and estimate their effect on the closure temperature. At the end all the results were compared with the published experiments so that to get the idea about the heterogeneity that can modify Pb diffusion behaviour in apatite which is a prime criterion for a robust (U—Th)/Pb age interpretation, thermal histories, and the closure temperature $T_c$ of apatite as $T_c$ is a function of the diffusion behaviour (Smye et al., 2018).
## 2. Methods and approaches

The stack workflow of our study is given in Fig. 1. A computational modelling technique called density functional theory (DFT) determines the energy of interacting atoms by determining their exterior electronic configuration (Kohn and Sham, 1965). In contrast to molecular dynamics, the DFT only addresses static problems involving atoms in certain places and does not address atoms dynamics. However, the quantum effects for the electrons are included in DFT and this allows for the most precise estimation of energy profiles/distribution for an ensemble of interacting atoms. In comparison to a real crystal, a supercell with a small number of atoms (usually <150 atoms) is used, and an infinite repeat of this supercell is believed to be representative of

an infinite lattice. The reasons for not selecting a supercell with large number of atoms is that it is computationally very expensive because of its high computation time and required cluster computing. An insertion site for Pb is identified by introducing it into a fully relaxed perfect lattice and the sites with the lowest energy are investigated after running another full geometry optimization. The stable/minimum energy positions of the inserted atom(s) can be determined by adjusting or nudging its position till the total energy is minimized for that inserted atom(s). At finite temperature, Pb atom start vibration and statistically moves across insertion sites. However, this is a dynamic process that falls outside the purview of DFT, which solely deals with static configurations. Therefore, another approach called the Transition State Theory (TST) is used which can statistically estimate the jumps between two adjacent energy

<figure><img src="image_3.png" /><figcaption>Fig. 1. workflow showing the stepwise procedure of the computational methods used in this work. The black arrow indicates the output of the method to be used as an input in the next step.</figcaption></figure>
## 2. Methods and approaches

The stack workflow of our study is given in Fig. 1. A computational modelling technique called density functional theory (DFT) determines the energy of interacting atoms by determining their exterior electronic configuration (Kohn and Sham, 1965). In contrast to molecular dynamics, the DFT only addresses static problems involving atoms in certain places and does not address atoms dynamics. However, the quantum effects for the electrons are included in DFT and this allows for the most precise estimation of energy profiles/distribution for an ensemble of interacting atoms. In comparison to a real crystal, a supercell with a small number of atoms (usually <150 atoms) is used, and an infinite repeat of this supercell is believed to be representative of

an infinite lattice. The reasons for not selecting a supercell with large number of atoms is that it is computationally very expensive because of its high computation time and required cluster computing. An insertion site for Pb is identified by introducing it into a fully relaxed perfect lattice and the sites with the lowest energy are investigated after running another full geometry optimization. The stable/minimum energy positions of the inserted atom(s) can be determined by adjusting or nudging its position till the total energy is minimized for that inserted atom(s). At finite temperature, Pb atom start vibration and statistically moves across insertion sites. However, this is a dynamic process that falls outside the purview of DFT, which solely deals with static configurations. Therefore, another approach called the Transition State Theory (TST) is used which can statistically estimate the jumps between two adjacent energy

<figure><img src="image_4.png" /><figcaption>Fig. 1. workflow showing the stepwise procedure of the computational methods used in this work. The black arrow indicates the output of the method to be used as an input in the next step.</figcaption></figure>
minima. This method is based on a comprehensive mapping of the energy field (provided by DFT) that an atom would face as it travels to its neighbouring site. The perfect crystal depicted by DFT and TST is not a natural crystal, which will always have defects, dislocations, distortions, substitutions, and usually fluid or mineral inclusions. So, at last of this study, we will determine how the perfect crystal modifies Pb diffusion and its closure temperature when additions to the pristine lattice are made i.e., the periodicity of the lattice is broken, and heterogeneity is introduced on length scales far greater than the unit cell dimension, alteration of the unit cell by interstitial atoms and chemical replacements/substitutions.

## 2.1. Density functional theory and nudged elastic band approach

In this study, the detailed stepwise procedure discussed in Bengtson et al. (2012), Djimbi et al. (2015), Gautheron et al. (2020), Bassal et al. (2022) and Ullah et al. (2023) is used for structural optimization, and insertion site analyses of Pb doped apatite lattices. It is based on the periodic DFT formulation of Hohenberg and Kohn (1964) and Kohn and Sham (1965) as implemented in the Quantum Espresso (QE) software package (Giannozzi et al., 2009). The apatite unit cell is projected in each spatial direction, creating a supercell with 84 atoms. This is done to reduce the volume relaxation effect caused by the inserted Pb atom and to prevent Pb from interacting with nearby cells. The insertion energies of Pb in each interstitial site are determined after relaxing the apatite structure using a gradient approach at constant volume. The energy needed to introduce a Pb atom from outside or inside (radioactive decay) the crystal into the interstitial site located inside that crystal as defined by the insertion energy which can be calculated as:

$$E_i = E_{apatite+Pb} - (E_{apatite} + E_{Pb}) \quad (1)$$

where $E_{Pb}$, $E_{apatite}$, and $E_{apatite+Pb}$ denotes the energy of the isolated Pb atom, the absolute energy of the perfect apatite supercell, and absolute energy of the Pb-dopped apatite host supercell respectively.

To calculate the migration energy between insertion sites, this study follows the Nudged Elastic Band (NEB) approach of Mills et al. (1995), and Jo'nnsson et al. (1998). NEB is a useful method for calculating the migratory energy or barrier energy profile (defined as: the energy difference between the energy of an atom at the saddle point and that when the same atom is at local minimum or insertion position) and minimal energy path between two interstitial sites. The NEB is a chain-of-states method executed by building a sequence of system images between states (Mills et al., 1995; Pratt, 1986). To maintain or ensure the path's continuity, adjacent images are connected by a spring force, imitating an elastic band. The spring constant of -5 eV/A° is used in this work for all the NEB calculation.

## 2.2. Transition state theory

An atom moves randomly from one site in the crystal lattice to a nearby symmetrically equivalent free site causing an interstitial diffusion. Diffusion, which is a thermally activated process because the diffusing atom requires enough energy to move and attempt to overcome an existing energy barrier. As was previously stated, DFT, which deals with static configurations, cannot directly calculate this dynamical (diffusion) process. However, the canonical ensemble is presumable in phase space since the motion takes place in thermodynamic equilibrium. This gives rise to the transition state hypothesis (Voter, 1986; Wert and Zener, 1949), in which the atomic motion is governed by jump probabilities, Γ, that depends exponentially on temperature.

$$\Gamma = \frac{\prod_{i=1}^{3} v_i}{\prod_{i=1}^{2} v_i^*} \exp\left(\frac{-E_{mig}}{kT}\right), \quad (2)$$

where, $v_i$ and $v_i^*$ indicates the normal mode frequencies at the initial local minimum and transition state positions respectively. $E_{mig}$ is the migration energy (kJ/mol) of the diffusing Pb atom. For a non-linear molecule there are 3 N-6 normal modes. As the energy is minimized at the insertion location via DFT structural optimization process therefore an upward curvature of energy is obtained for every displacement of the inserted atom (which in this case is Pb). Three eigenvalues, which can directly be translated to the three common frequencies, can be used to express this 3D-curvature in same manner as for any harmonic oscillator (e.g., pendulum). At the saddle point, the same characteristics apply, but there are only two frequencies because only transverse motion in relation to the transition route is involved. As per Vineyard (1957) approximation, the energetics of other atoms remain unaffected by the position of the migratory atom being taken into account for one element (i.e., Pb). The TST thus offers a clear theoretical method for determining the frequency factor.

## 2.3. Kinetic Monte Carlo (KMC) simulation method

KMC simulates a random walk over the interstitial sites at a specific temperature by using microscopic jump probabilities determined from TST using DFT results (Bortz et al., 1975; Hoffmann et al., 2014; Andersen et al., 2019; Fig. 1). Further, it also makes the computation easier for all the subjected heterogeneities like chemical changes, obstructions, traps, and anisotropy in the lattice. To start the simulation, a Pb atom is randomly inserted into one available interstitial site for a chosen temperature. Then, using the estimated jumping frequencies provided by the TST, a random value is chosen to mimic the jumping duration for each of the potential diffusion directions. The simulation extracts the jump direction with the lowest possible time, updates the time and coordinates, and adds the time to the total travel time at each step. This is done m times for each atom, and n atoms or paths (atoms trajectories) are simulated (Djimbi et al., 2015). The total time is calculated for each trajectory as the sum of the residency times at each site between jumps. Let's $D_x$, $D_y$, $D_z$, denotes the diffusion-coefficients along the x, y, and z axes, then in accordance with the Einstein relation as:

$$D_x = \frac{\langle x_{MSD}^2 \rangle}{2 \langle t \rangle}, \quad (3)$$

$$D_y = \frac{\langle y_{MSD}^2 \rangle}{2 \langle t \rangle}, \quad (4)$$

$$D_z = \frac{\langle z_{MSD}^2 \rangle}{2 \langle t \rangle}. \quad (5)$$

Here, $\langle x_{MSD}^2 \rangle$, $\langle y_{MSD}^2 \rangle$, and $\langle z_{MSD}^2 \rangle$ are the averages of the mean square displacements (MSD) over their respective trajectories along the x, y, and z directions and $\langle t \rangle$ represents the average time taken by the diffusing atom along the trajectories. From the time dependence of the mean square displacement of the Pb atom, $\langle r^2(t) \rangle$ at temperature T, the effective (isotropic) diffusion coefficient at that temperature, $D_t$, was calculated (Finger and Hazen, 1980; Pozun and Henkelman, 2011).

$$\langle r^2(t) \rangle = 2d\langle t \rangle D_t \quad (6)$$

where d indicates the dimension, which is equal to 3 in case of 3D random walk diffusion of Pb in apatite crystal. The diffusion coefficient at various temperatures was calculated using Eqs. (4)–(7), and diffusion coefficient (D) equations are modelled using Arrhenius law:

$$D = D_o \exp\left(-E_a/k_b T\right) \quad (7)$$
minima. This method is based on a comprehensive mapping of the energy field (provided by DFT) that an atom would face as it travels to its neighbouring site. The perfect crystal depicted by DFT and TST is not a natural crystal, which will always have defects, dislocations, distortions, substitutions, and usually fluid or mineral inclusions. So, at last of this study, we will determine how the perfect crystal modifies Pb diffusion and its closure temperature when additions to the pristine lattice are made i.e., the periodicity of the lattice is broken, and heterogeneity is introduced on length scales far greater than the unit cell dimension, alteration of the unit cell by interstitial atoms and chemical replacements/substitutions.

## 2.1. Density functional theory and nudged elastic band approach

In this study, the detailed stepwise procedure discussed in Bengtson et al. (2012), Djimbi et al. (2015), Gautheron et al. (2020), Bassal et al. (2022) and Ullah et al. (2023) is used for structural optimization, and insertion site analyses of Pb dopped apatite lattices. It is based on the periodic DFT formulation of Hohenberg and Kohn (1964) and Kohn and Sham (1965) as implemented in the Quantum Espresso (QE) software package (Giannozzi et al., 2009). The apatite unit cell is projected in each spatial direction, creating a supercell with 84 atoms. This is done to reduce the volume relaxation effect caused by the inserted Pb atom and to prevent Pb from interacting with nearby cells. The insertion energies of Pb in each interstitial site are determined after relaxing the apatite structure using a gradient approach at constant volume. The energy needed to introduce a Pb atom from outside or inside (radioactive decay) the crystal into the interstitial site located inside that crystal as defined by the insertion energy which can be calculated as:

$$E_i = E_{apatite+Pb} - (E_{apatite} + E_{Pb}) \quad (1)$$

where $E_{Pb}$, $E_{apatite}$, and $E_{apatite+Pb}$ denotes the energy of the isolated Pb atom, the absolute energy of the perfect apatite supercell, and absolute energy of the Pb-dopped apatite host supercell respectively.

To calculate the migration energy between insertion sites, this study follows the Nudged Elastic Band (NEB) approach of Mills et al. (1995), and Jo'nnsson et al. (1998). NEB is a useful method for calculating the migratory energy or barrier energy profile (defined as: the energy difference between the energy of an atom at the saddle point and that when the same atom is at local minimum or insertion position) and minimal energy path between two interstitial sites. The NEB is a chain-of-states method executed by building a sequence of system images between states (Mills et al., 1995; Pratt, 1986). To maintain or ensure the path's continuity, adjacent images are connected by a spring force, imitating an elastic band. The spring constant of -5 eV/A° is used in this work for all the NEB calculation.

## 2.2. Transition state theory

An atom moves randomly from one site in the crystal lattice to a nearby symmetrically equivalent free site causing an interstitial diffusion. Diffusion, which is a thermally activated process because the diffusing atom requires enough energy to move and attempt to overcome an existing energy barrier. As was previously stated, DFT, which deals with static configurations, cannot directly calculate this dynamical (diffusion) process. However, the canonical ensemble is presumable in phase space since the motion takes place in thermodynamic equilibrium. This gives rise to the transition state hypothesis (Voter, 1986; Wert and Zener, 1949), in which the atomic motion is governed by jump probabilities, Γ, that depends exponentially on temperature.

$$\Gamma = \frac{\prod_{i=1}^{3} v_i}{\prod_{i=1}^{2} v_i^*} \exp\left(\frac{-E_{mig}}{kT}\right), \quad (2)$$

where, $v_i$ and $v_i^*$ indicates the normal mode frequencies at the initial local minimum and transition state positions respectively. $E_{mig}$ is the migration energy (kJ/mol) of the diffusing Pb atom. For a non-linear molecule there are 3 N-6 normal modes. As the energy is minimized at the insertion location via DFT structural optimization process therefore an upward curvature of energy is obtained for every displacement of the inserted atom (which in this case is Pb). Three eigenvalues, which can directly be translated to the three common frequencies, can be used to express this 3D-curvature in same manner as for any harmonic oscillator (e.g., pendulum). At the saddle point, the same characteristics apply, but there are only two frequencies because only transverse motion in relation to the transition route is involved. As per Vineyard (1957) approximation, the energetics of other atoms remain unaffected by the position of the migratory atom being taken into account for one element (i.e., Pb). The TST thus offers a clear theoretical method for determining the frequency factor.

## 2.3. Kinetic Monte Carlo (KMC) simulation method

KMC simulates a random walk over the interstitial sites at a specific temperature by using microscopic jump probabilities determined from TST using DFT results (Bortz et al., 1975; Hoffmann et al., 2014; Andersen et al., 2019; Fig. 1). Further, it also makes the computation easier for all the subjected heterogeneities like chemical changes, obstructions, traps, and anisotropy in the lattice. To start the simulation, a Pb atom is randomly inserted into one available interstitial site for a chosen temperature. Then, using the estimated jumping frequencies provided by the TST, a random value is chosen to mimic the jumping duration for each of the potential diffusion directions. The simulation extracts the jump direction with the lowest possible time, updates the time and coordinates, and adds the time to the total travel time at each step. This is done m times for each atom, and n atoms or paths (atoms trajectories) are simulated (Djimbi et al., 2015). The total time is calculated for each trajectory as the sum of the residency times at each site between jumps. Let's $D_x$, $D_y$, $D_z$, denotes the diffusion-coefficients along the x, y, and z axes, then in accordance with the Einstein relation as:

$$D_x = \frac{\langle x_{MSD}^2 \rangle}{2\langle t \rangle}, \quad (3)$$

$$D_y = \frac{\langle y_{MSD}^2 \rangle}{2\langle t \rangle}, \quad (4)$$

$$D_z = \frac{\langle z_{MSD}^2 \rangle}{2\langle t \rangle}. \quad (5)$$

Here, $\langle x_{MSD}^2 \rangle$, $\langle y_{MSD}^2 \rangle$, and $\langle z_{MSD}^2 \rangle$ are the averages of the mean square displacements (MSD) over their respective trajectories along the x, y, and z directions and $\langle t \rangle$ represents the average time taken by the diffusing atom along the trajectories. From the time dependence of the mean square displacement of the Pb atom, $\langle r^2(t) \rangle$ at temperature T, the effective (isotropic) diffusion coefficient at that temperature, $D_t$, was calculated (Finger and Hazen, 1980; Pozun and Henkelman, 2011).

$$\langle r^2(t) \rangle = 2d\langle t \rangle D_t \quad (6)$$

where d indicates the dimension, which is equal to 3 in case of 3D random walk diffusion of Pb in apatite crystal. The diffusion coefficient at various temperatures was calculated using Eqs. (4)–(7), and diffusion coefficient (D) equations are modelled using Arrhenius law:

$$D = D_o \exp\left(-\frac{E_a}{k_b T}\right) \quad (7)$$
## 3. Results

### 3.1. Pb interstitial sites delineation in apatite

A full geometric optimization on apatite supercell is applied using periodic DFT and its structure are fully relaxed. The model is then utilized to pinpoint all potential interstitial insertion locations for Pb atoms along with diffusion pathways. In order to replicate and investigate the Pb insertion in the pure apatite crystal lattice free of boundary effects, a 1 × 1 × 2 supercells are used as the host lattices (Fig. 2-3). In the fully optimized supercell of apatite, the two adjacent Ca(2) and Ca(1) ployhedra connect through O atoms shared with PO₄ tetrahedra. Such connections of both polyhedral and PO₄ tetrahedra generates two preferential channels for the diffusion of Pb in apatite minerals identical to those discussed in Bengtson et al. (2012) and Djimbi et al. (2015) for the diffusion of helium (He) in apatite. The first one is situated along the F atoms that run parallel to the z-axis. This channel is made up of six structurally and symmetrically equivalent sub-channels according to the apatite's hexagonal structure (Fig. 3b-c). The second channel, which is more difficult to identify, runs along a direction inside the plane approximately orthogonal to z-axis and contains the spatial coordinates x and y. For this in-plane diffusion, there are three preferential directions, each at a 60-degree angle. Two distinct Pb interstitial sites; interstitial site along z-axis (ISZ) and in-plane interstitial site (ISIP), are examined during DFT calculations (Fig. 3a). Two symmetrically equivalent ISZ sites are separated by an energy barrier along the z-axis of the apatite supercell model, and this pattern is periodically repeated throughout its crystal structure. The Pb atom is situated at ISZ sites at a mean distance of 2.89 Å from six Ca1 atoms, 2.95 Å from six O atoms, and 1.75 Å from the two nearby fluorine atoms. For the in-plane diffusion the insertion site (ISIP) is demarcated inside the channel that is nearly orthogonal to the fluorine axis. As per the hexagonal geometry of apatite lattice, this channel exists along the in-between gaps of Ca1 and Ca2 polyhedra in [110] directions (i.e. Pb is placed in-between oppositely equivalent Ca(1) and Ca(2) coordination polyhedra) as shown in Fig. 3b. In this second type (ISIP) of Pb insertion sites, the Pb atom is

placed in relation to the surrounding two Ca1, two Ca2, two P, and six O atoms at mean distances of 2.97, 2.74, 2.06 Å, and 2.15 Å respectively. The local relaxation before and after geometry optimization at both interstitial sites are shown in Fig. 4 and Table 1.

Similarly, for apatite with 10% Na concentration i.e., Na₀.₁-apatite and for apatite with 10% Na and 25% Cl concentration i.e., Cl₀.₂₅-Na₀.₁-apatite, the interstitial sites were also calculated to further investigate the effect of chemical substitution on the insertion as well as activation energies. This Na and Cl concentrations are used to bring the chemistry of the apatite closer to that of natural Durango apatite (with chemical composition Ca₉.₃₃ Na₀.₀₈ Sr₀.₀₁ REE₀.₀₉ (PO₄)₅.₈₇ (SO₄)₀.₀₅ (SiO₄)₀.₀₆ (AsO)₀.₀₁ (CO₃)₀.₀₁ F₁.₉₀ Cl₀.₁₂ (OH)₀.₀₁) for comparison (which is approximately close ideal F₃.₇₅-Cl₀.₂₅-apatite as other substituents is almost negligible). The overall heterogeneity cannot be incorporated in single computation as doing so would entail expanding the supercell and therefore will increase the number of atoms which is specifically limited to or even can't be handle by DFT. The Pb atom's insertion energies in F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite are -2.03 eV, -2.21 eV, and -1.94 eV respectively. Here, the negative energies imply that those places are energetically favoured. We can infer two diffusion paths in pure F-apatite from the aforementioned findings: (1) Diffusion from an ISZ site to the next equivalent ISZ site along the fluorine or z-axis passing through a barrier via three potential sub-channels (Figs. 2 and 4) and (2) Diffusion from one ISZ (or ISZ barrier) site to the next symmetrically equivalent ISZ (or ISZ barrier) site along the in-plane diraction while passing through ISIP interstitial site (Fig. 3c).

### 3.2. Estimation of the activation energies

In order to simulate the migration of a Pb atoms along the z-axis in pure F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite, a climbing image NEB calculation was carried out in one of the six equivalent sub-channels as shown in Figs. 5-7. This climbing technique is more precise than the standard NEB one in that it yields a strict convergence to a saddle position. The results obtained through this method shows that the migration of Pb from one energy minima to the next equivalent one

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** (a): crystal structure of pure F-apatite with z-axis onto the plane and x- and y-axes on the plane of paper. A network of hexagonally arranged Ca(1) and Ca(2) coordination polyhedral (brown colour) are labelled on the lattice. The Ca(1)O9 is characterized by a trigonal prism while the Ca(2) polyhedron show an octahedron (CaO5X) type geometry (b): A (1 × 1 × 2) supercell of perfect F-apatite structure comprises of 84 atoms showing a 3D projected image. The periodic distribution of the phosphorus tetrahedra showing a zigzag pattern along the z-axis can be seen in orange colour. (All the colour visuals can be seen in the web version of this manuscript). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
# 3. Results

## 3.1. Pb interstitial sites delineation in apatite

A full geometric optimization on apatite supercell is applied using periodic DFT and its structure are fully relaxed. The model is then utilized to pinpoint all potential interstitial insertion locations for Pb atoms along with diffusion pathways. In order to replicate and investigate the Pb insertion in the pure apatite crystal lattice free of boundary effects, a 1 × 1 × 2 supercells are used as the host lattices (Fig. 2-3). In the fully optimized supercell of apatite, the two adjacent Ca(2) and Ca(1) ployhedra connect through O atoms shared with PO₄ tetrahedra. Such connections of both polyhedral and PO₄ tetrahedra generates two preferential channels for the diffusion of Pb in apatite minerals identical to those discussed in Bengtson et al. (2012) and Djimbi et al. (2015) for the diffusion of helium (He) in apatite. The first one is situated along the F atoms that run parallel to the z-axis. This channel is made up of six structurally and symmetrically equivalent sub-channels according to the apatite's hexagonal structure (Fig. 3b-c). The second channel, which is more difficult to identify, runs along a direction inside the plane approximately orthogonal to z-axis and contains the spatial coordinates x and y. For this in-plane diffusion, there are three preferential directions, each at a 60-degree angle. Two distinct Pb interstitial sites; interstitial site along z-axis (ISZ) and in-plane interstitial site (ISIP), are examined during DFT calculations (Fig. 3a). Two symmetrically equivalent ISZ sites are separated by an energy barrier along the z-axis of the apatite supercell model, and this pattern is periodically repeated throughout its crystal structure. The Pb atom is situated at ISZ sites at a mean distance of 2.89 Å from six Ca1 atoms, 2.95 Å from six O atoms, and 1.75 Å from the two nearby fluorine atoms. For the in-plane diffusion the insertion site (ISIP) is demarcated inside the channel that is nearly orthogonal to the fluorine axis. As per the hexagonal geometry of apatite lattice, this channel exists along the in-between gaps of Ca1 and Ca2 polyhedra in [110] directions (i.e. Pb is placed in-between oppositely equivalent Ca(1) and Ca(2) coordination polyhedra) as shown in Fig. 3b. In this second type (ISIP) of Pb insertion sites, the Pb atom is

placed in relation to the surrounding two Ca1, two Ca2, two P, and six O atoms at mean distances of 2.97, 2.74, 2.06 Å, and 2.15 Å respectively. The local relaxation before and after geometry optimization at both interstitial sites are shown in Fig. 4 and Table 1.

Similarly, for apatite with 10% Na concentration i.e., Na₀.₁-apatite and for apatite with 10% Na and 25% Cl concentration i.e., Cl₀.₂₅-Na₀.₁-apatite, the interstitial sites were also calculated to further investigate the effect of chemical substitution on the insertion as well as activation energies. This Na and Cl concentrations are used to bring the chemistry of the apatite closer to that of natural Durango apatite (with chemical composition Ca₉.₃₃ Na₀.₀₈ Sr₀.₀₁ REE₀.₀₉ (PO₄)₅.₈₇ (SO₄)₀.₀₅ (SiO₄)₀.₀₆ (AsO)₀.₀₁ (CO₃)₀.₀₁ F₁.₉₀ Cl₀.₁₂ (OH)₀.₀₁) for comparison (which is approximately close ideal F₃.₇₅-Cl₀.₂₅-apatite as other substituents is almost negligible). The overall heterogeneity cannot be incorporated in single computation as doing so would entail expanding the supercell and therefore will increase the number of atoms which is specifically limited to or even can't be handle by DFT. The Pb atom's insertion energies in F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite are -2.03 eV, -2.21 eV, and -1.94 eV respectively. Here, the negative energies imply that those places are energetically favoured. We can infer two diffusion paths in pure F-apatite from the aforementioned findings: (1) Diffusion from an ISZ site to the next equivalent ISZ site along the fluorine or z-axis passing through a barrier via three potential sub-channels (Figs. 2 and 4) and (2) Diffusion from one ISZ (or ISZ barrier) site to the next symmetrically equivalent ISZ (or ISZ barrier) site along the in-plane dirtection while passing through ISIP interstitial site (Fig. 3c).

## 3.2. Estimation of the activation energies

In order to simulate the migration of a Pb atoms along the z-axis in pure F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite, a climbing image NEB calculation was carried out in one of the six equivalent sub-channels as shown in Figs. 5-7. This climbing technique is more precise than the standard NEB one in that it yields a strict convergence to a saddle position. The results obtained through this method shows that the migration of Pb from one energy minima to the next equivalent one

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** (a): crystal structure of pure F-apatite with z-axis onto the plane and x- and y-axes on the plane of paper. A network of hexagonally arranged Ca(1) and Ca(2) coordination polyhedral (brown colour) are labelled on the lattice. The Ca(1)O9 is characterized by a trigonal prism while the Ca(2) polyhedron show an octahedron (CaO5X) type geometry (b): A (1 × 1 × 2) supercell of perfect F-apatite structure comprises of 84 atoms showing a 3D projected image. The periodic distribution of the phosphorus tetrahedra showing a zigzag pattern along the z-axis can be seen in orange colour. (All the colour visuals can be seen in the web version of this manuscript). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** A supercell of F-apatite comprises of 85 atoms used in the DFT computations. (a) showing the two possible sites for the incorporation of Pb (labelled as ISZ (interstitial site along z-axis), and in-plane interstitial site (ISIP) along the plane orthogonal to z-axis) (b) Pathway 1 is located along the z-axis is along in the [001] direction. The migration of Pb along this path encounter barriers when moving from one energy minima site to the next symmetrically equivalent energy minima site. (c) Pathway 2 exist along the plane orthogonal to the z-axis called here the in-plane diffusion path. The migration of Pb along this path also encounter barriers before reaching the z-axis sites. In (b) and (c), black spheres represent Pb atoms located at equidistance points (images) running along their respective pathways. At each of these points the energy is computed via NEB method. Note that, the computations were performed with just one Pb atom per supercell, the cells are depicted with numerous Pb atoms just for the purpose to highlight probable Pb interstitial locations along each diffusion pathway (All the colour visuals can be seen in the web version of this manuscript).</figcaption></figure>

encounter a barrier of about 5.18 eV (499.7 kJ/mol), 2.8 eV (270.3 kJ/mol), and 5.73 eV (552.8 kJ/mol), for pure F-apatite, Na<sub>0.1</sub>-apatite, and Cl<sub>0.25</sub>- Na<sub>0.1</sub>-apatite respectively along the z-axis (Figs. 5a, 6a, 7a). For a Pb atom to diffuse along this z-direction in all these mentioned apatites, there are six equal channels (i.e., three in the positive z-direction and three in negative z-direction) and therefore the barriers in all these channels are the same for Pb atom. In the orthogonal plane one NEB calculation has been performed, with a pathway that connect one ISZ point in xy-direction to the opposite equivalent ISZ point travelling through ISIP as shown in Fig. 3c and Figs. 5b-7b. The migration barrier along this orthogonal plane pathway is 8.05 eV (776.5 kJ/mol), 4.58 eV (441.5 kJ/mol), and 7.79 eV (751.6 kJ/mol), for pure F-apatite, Na<sub>0.1</sub>-apatite, and Cl<sub>0.25</sub>- Na<sub>0.1</sub>-apatite respectively. Additionally, the activation energies for the F-apatite and Cl<sub>0.25</sub>- Na<sub>0.1</sub>-apatite is significantly higher than the activation energies deduced from the experiments of natural apatite, while that for Na<sub>0.1</sub>-apatite is almost within the range of these published experimental results of natural apatite of Watson *et al.* (1985) and Cherniak *et al.* (1991) with activation energies 292.8 kJ/mol and 228.4 kJ/mol respectively as shown in Table 2.

The energy profile is altered by the substitution of sodium (in place Ca(2)) and chlorine. In fact, the Pb atom must overcome three obstacles to go to the next equivalent lattice site ISZ from initial ISZ along the z-axis. The first barrier along z-axis in Na<sub>0.1</sub>-apatite is the lowest one and is equal to 207.4 kJ/mol, which represents the jump of an atom up the sodium atom from one ISZ to the next adjacent ISZ. The second and third jumps of Pb is the same as that in F-apatite and is correspond to a barrier value of 499.7 kJ/mol. As two of the Ca(2) atoms are replace with Na in Na<sub>0.1</sub>-apatite so for in-plane diffusion the Pb encounter a barrier of 441.5 kJ/mol to reach the next equivalent site till it reach the ISZ site. So, in case of 10%Na the Pb diffusion in apatite will be mostly along the in-plane direction as compared to z-direction. This substitution of sodium in place of Ca(2) in F-apatite makes the diffusion of Pb much easier along the z-direction as compared to F-apatite. For much greater concentration of Na both the in-plane as well as the z-axis diffusion of Pb become faster than the z-axis diffusion in pure F-apatite and this is discussed latter in much detail.

Similarly, if one of the F atoms (the one that is attached to the sodium atom) in Na<sub>0.1</sub>-apatite is replaced by Cl such that in Cl<sub>0.25</sub>- Na<sub>0.1</sub>-apatite, then one of the three barrier has a value of 5.73 eV (552.8 kJ/mol) while the remaining two have same values as that in F-apatite (i.e. 5.18 eV (499.7 kJ/mol)). A chlorine atom is placed every two meshes in this structure, although with other chlorine concentrations, the greater

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** indicating the local environment of lead (Pb) guest atom in the in-plane interstitial site. Red balls represents the nearest neighbouring O atoms, the brown colour balls are nearest neighbouring Ca atoms, green are the F atoms, and the orange one are the nearest neighbouring phosphorus atoms (a) shows the pre-optimized F-apatite structure with Pb occupying the ISIP interstitial site, (b) is the DFT-optimized configuration of F-apatite characterized by a strong local structural relaxation around the guest Pb atom. (All the colour visuals can be seen in the web version of this manuscript). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** A supercell of F-apatite comprises of 85 atoms used in the DFT computations. (a) showing the two possible sites for the incorporation of Pb (labelled as ISZ (interstitial site along z-axis), and in-plane interstitial site (ISIP) along the plane orthogonal to z-axis) (b) Pathway 1 is located along the z-axis is along in the [001] direction. The migration of Pb along this path encounter barriers when moving from one energy minima site to the next symmetrically equivalent energy minima site. (c) Pathway 2 exist along the plane orthogonal to the z-axis called here the in-plane diffusion path. The migration of Pb along this path also encounter barriers before reaching the z-axis sites. In (b) and (c), black spheres represent Pb atoms located at equidistance points (images) running along their respective pathways. At each of these points the energy is computed via NEB method. Note that, the computations were performed with just one Pb atom per supercell, the cells are depicted with numerous Pb atoms just for the purpose to highlight probable Pb interstitial locations along each diffusion pathway (All the colour visuals can be seen in the web version of this manuscript).</figcaption></figure>

encounter a barrier of about 5.18 eV (499.7 kJ/mol), 2.8 eV (270.3 kJ/mol), and 5.73 eV (552.8 kJ/mol), for pure F-apatite, Na₀.₁-apatite, and Cl₀.₂₅- Na₀.₁-apatite respectively along the z-axis (Figs. 5a, 6a, 7a). For a Pb atom to diffuse along this z-direction in all these mentioned apatites, there are six equal channels (i.e., three in the positive z-direction and three in negative z-direction) and therefore the barriers in all these channels are the same for Pb atom. In the orthogonal plane one NEB calculation has been performed, with a pathway that connect one ISZ point in the xy-direction to the opposite equivalent ISZ point travelling through ISIP as shown in Fig. 3c and Figs. 5b-7b. The migration barrier along this orthogonal plane pathway is 8.05 eV (776.5 kJ/mol), 4.58 eV (441.5 kJ/mol), and 7.79 eV (751.6 kJ/mol), for pure F-apatite, Na₀.₁-apatite, and Cl₀.₂₅- Na₀.₁-apatite respectively. Additionally, the activation energies for the F-apatite and Cl₀.₂₅- Na₀.₁-apatite is significantly higher than the activation energies deduced from the experiments of natural apatite, while that for Na₀.₁-apatite is almost within the range of these published experimental results of natural apatite of Watson et al. (1985) and Cherniak et al. (1991) with activation energies 292.8 kJ/mol and 228.4 kJ/mol respectively as shown in Table 2.

The energy profile is altered by the substitution of sodium (in place Ca(2)) and chlorine. In fact, the Pb atom must overcome three obstacles to go to the next equivalent lattice site ISZ from initial ISZ along the z-axis. The first barrier along z-axis in Na₀.₁-apatite is the lowest one and is equal to 207.4 kJ/mol, which represents the jump of an atom up the sodium atom from one ISZ to the next adjacent ISZ. The second and third jumps of Pb is the same as that in F-apatite and is correspond to a barrier value of 499.7 kJ/mol. As two of the Ca(2) atoms are replace with Na in Na₀.₁-apatite so for in-plane diffusion the Pb encounter a barrier of 441.5 kJ/mol to reach the next equivalent site till it reach the ISZ site. So, in case of 10%Na the Pb diffusion in apatite will be mostly along the in-plane direction as compared to z-direction. This substitution of sodium in place of Ca(2) in F-apatite makes the diffusion of Pb much easier along the z-direction as compared to F-apatite. For much greater concentration of Na both the in-plane as well as the z-axis diffusion of Pb become faster than the z-axis diffusion in pure F-apatite and this is discussed latter in much detail.

Similarly, if one of the F atoms (the one that is attached to the sodium atom) in Na₀.₁-apatite is replaced by Cl such that in Cl₀.₂₅- Na₀.₁-apatite, then one of the three barrier has a value of 5.73 eV (552.8 kJ/mol) while the remaining two have same values as that in F-apatite (i.e. 5.18 eV (499.7 kJ/mol)). A chlorine atom is placed every two meshes in this structure, although with other chlorine concentrations, the greater

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** indicating the local environment of lead (Pb) guest atom in the in-plane interstitial site. Red balls represents the nearest neighbouring O atoms, the brown colour balls are nearest neighbouring Ca atoms, green are the F atoms, and the orange one are the nearest neighbouring phosphorus atoms (a) shows the pre-optimized F-apatite structure with Pb occupying the ISIP interstitial site, (b) is the DFT-optimized configuration of F-apatite characterized by a strong local structural relaxation around the guest Pb atom. (All the colour visuals can be seen in the web version of this manuscript). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<table><caption>**Table 1**
The local environment of the guest Pb atom at ISZ site and ISIP interstitial sites before and after complete geometry optimization by using DFT method.</caption><thead><tr><th></th><th>Pb-Ca1 (Å)</th><th>% Change</th><th>Pb-Ca2 (Å)</th><th>% Change</th><th>Pb-F (Å)</th><th>% Change</th><th>Pb-O (Å)</th><th>% Change</th></tr></thead><tbody><tr><td>Along F-atom site (ISZ)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Before</td><td>2.55</td><td></td><td>–</td><td>–</td><td>1.69</td><td></td><td>2.42</td><td></td></tr><tr><td>Optimized</td><td>3.17</td><td>24.3</td><td>–</td><td>–</td><td>2.25</td><td>33.1</td><td>2.69</td><td>11.2</td></tr><tr><td>Along in-plane (ISIP)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Before</td><td>2.97</td><td></td><td>2.74</td><td></td><td>2.06</td><td></td><td>2.15</td><td></td></tr><tr><td>Optimized</td><td>3.13</td><td>5.4</td><td>3.27</td><td>19.3</td><td>3.04</td><td>47.6</td><td>2.75</td><td>27.9</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** The activation-energy profiles of the Pb energetics in perfect F-apatite lattice. Each profile indicates the migration energy required for a Pb atom to move between neighbouring symmetrically equivalent minima. Three activation energy profiles are plotted here i.e. (a): is along z-axis or (001) direction, where Pb needs to overcome the energy barrier encountered at sites along z-axis, (b) is along (110) or in-plane directions, where Pb needs to cross the energy barrier encountered at sites along the orthogonal plane direction, and (c) shows the activation energy profile for the whole path that connect (110) direction with (001) indicating that if the Pb start diffusion first in the in-plane direction and when it encounter a site along the z-axis during the travel, the diffusion will shift again along the z-axis as can be seen in the energy profile trend.</figcaption></figure>

barrier of 5.73 eV (552.8 kJ/mol) would still apply. Given those barriers, one would anticipate quite slower 3D diffusion in the chlorine-substituted apatite crystal. Such high barrier offered by Cl along the [001] direction prevents diffusion of Pb atom. However, there are still calcium and fluorine-type barriers present along the z-axis in both Na<sub>0.1</sub>-apatite, and Cl<sub>0.25</sub>-Na<sub>0.1</sub>-apatite, indicating that the Pb atoms must work their way around and essentially jumps towards the lowest barriers in other directions if present. As a result, diffusion along the z-axis may take place in segments liked by workaround of Na, F, and Cl substitutions. Therefore, dealing with such complex problems requires the development of a suitable method i.e., 3D Random Walk KMC code for determining the effective diffusion coefficient (which is discussed in detail in later sections and in Djimbi et al., 2015).

### 3.3. Estimating the hopping rates

By adjusting the Pb atom's location at the stable site and at the saddle

point, which allows the relaxation of the surrounding atoms, the frequencies involved in eq. 2 have been calculated from the DFT based derived energy landscape. The increase in energy depicted by the displacements' quadratic form, yields the eigen-frequencies after diagonalization. In accordance with the TST, the jump probability Γ, with respect to temperature can be expressed as in Eqs. (8)–(13) after

determining the attempt frequency $v$ (in accordance with $v = \frac{\prod_{i=1}^{3} v_i}{\prod_{i=1}^{3} v_i^*}$).

$$ \text{For F-apatite along z-axis: } \Gamma = 1.98 \times 10^{12} \exp\left(-\frac{5.18}{kT}\right) (s^{-1}), \quad (8) $$

$$ \text{For F-apatite along in-plane: } \Gamma = 1.08 \times 10^{12} \exp\left(-\frac{8.05}{kT}\right) (s^{-1}), \quad (9) $$
<table><caption>**Table 1**
The local environment of the guest Pb atom at ISZ site and ISIP interstitial sites before and after complete geometry optimization by using DFT method.</caption><thead><tr><th></th><th>Pb-Ca1 (Å)</th><th>% Change</th><th>Pb-Ca2 (Å)</th><th>% Change</th><th>Pb-F (Å)</th><th>% Change</th><th>Pb-O (Å)</th><th>% Change</th></tr></thead><tbody><tr><td>Along F-atom site (ISZ)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Before</td><td>2.55</td><td></td><td>-</td><td>-</td><td>1.69</td><td></td><td>2.42</td><td></td></tr><tr><td>Optimized</td><td>3.17</td><td>24.3</td><td>-</td><td>-</td><td>2.25</td><td>33.1</td><td>2.69</td><td>11.2</td></tr><tr><td>Along in-plane (ISIP)</td><td></td><td></td><td></td><td></td><td>Pb-P (Å)</td><td></td><td></td><td></td></tr><tr><td>Before</td><td>2.97</td><td></td><td>2.74</td><td></td><td>2.06</td><td></td><td>2.15</td><td></td></tr><tr><td>Optimized</td><td>3.13</td><td>5.4</td><td>3.27</td><td>19.3</td><td>3.04</td><td>47.6</td><td>2.75</td><td>27.9</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** The activation-energy profiles of the Pb energetics in perfect F-apatite lattice. Each profile indicates the migration energy required for a Pb atom to move between neighbouring symmetrically equivalent minima. Three activation energy profiles are plotted here i.e. (a): is along z-axis or (001) direction, where Pb needs to overcome the energy barrier encountered at sites along z-axis, (b) is along (110) or in-plane directions, where Pb needs to cross the energy barrier encountered at sites along the orthogonal plane direction, and (c) shows the activation energy profile for the whole path that connect (110) direction with (001) indicating that if the Pb start diffusion first in the in-plane direction and when it encounter a site along the z-axis during the travel, the diffusion will shift again along the z-axis as can be seen in the energy profile trend.</figcaption></figure>

barrier of 5.73 eV (552.8 kJ/mol) would still apply. Given those barriers, one would anticipate quite slower 3D diffusion in the chlorine-substituted apatite crystal. Such high barrier offered by Cl along the [001] direction prevents diffusion of Pb atom. However, there are still calcium and fluorine-type barriers present along the z-axis in both Na<sub>0.1</sub>-apatite, and Cl<sub>0.25</sub>-Na<sub>0.1</sub>-apatite, indicating that the Pb atoms must work their way around and essentially jumps towards the lowest barriers in other directions if present. As a result, diffusion along the z-axis may take place in segments liked by workaround of Na, F, and Cl substitutions. Therefore, dealing with such complex problems requires the development of a suitable method i.e., 3D Random Walk KMC code for determining the effective diffusion coefficient (which is discussed in detail in later sections and in Djimbi et al., 2015).

### 3.3. Estimating the hopping rates

By adjusting the Pb atom's location at the stable site and at the saddle

point, which allows the relaxation of the surrounding atoms, the frequencies involved in eq. 2 have been calculated from the DFT based derived energy landscape. The increase in energy depicted by the displacements' quadratic form, yields the eigen-frequencies after diagonalization. In accordance with the TST, the jump probability Γ, with respect to temperature can be expressed as in Eqs. (8)–(13) after

determining the attempt frequency $v$ (in accordance with $v = \frac{\prod_{i=1}^{3} v_i}{\prod_{i=1}^{3} v_i^*}$).

$$ \text{For F-apatite along z-axis: } \Gamma = 1.98 \times 10^{12} \exp\left(-\frac{5.18}{kT}\right) (s^{-1}), \quad (8) $$

$$ \text{For F-apatite along in-plane: } \Gamma = 1.08 \times 10^{12} \exp\left(-\frac{8.05}{kT}\right) (s^{-1}), \quad (9) $$
<figure><img src="image_8.png" /><figcaption>Fig. 6. The activation-energy profiles of the Pb energetics in perfect Na₀.₁-apatite lattice. Each profile indicates the migration energy required for a Pb atom to move between neighbouring symmetrically equivalent minima. Three activation energy profiles are plotted here i.e. (a): is along z-axis or (001) direction, where Pb needs to overcome the energy barrier encountered at sites along z-axis, (b) is along (110) or in-plane directions, where Pb needs to cross the energy barrier encountered at sites along the orthogonal plane direction, and (c) shows the activation energy profile for the whole path that connect (110) direction with (001) indicating that if the Pb start diffusion first in the in-plane direction and when it encounter a site along the z-axis during the travel, the diffusion will shift again along the z-axis as can be seen in the energy profile trend.</figcaption></figure>

$$ \text{For Na}_{0.1} \text{ - apatite along z-axis : } \Gamma = 4.36 \times 10^{12} \exp\left(\frac{-2.8}{kT}\right) (s^{-1}), \quad (10) $$

$$ \text{For Na}_{0.1} \text{ - apatite along in-plane : } \Gamma = 3.51 \times 10^{12} \exp\left(\frac{-4.58}{kT}\right) (s^{-1}), \quad (11) $$

$$ \text{For Cl}_{0.25} - \text{Na}_{0.1} \text{ - apatite along z-axis : } \Gamma \\ = 3.51 \times 10^{12} \exp\left(\frac{-5.73}{kT}\right) (s^{-1}), \quad (12) $$

$$ \text{For Cl}_{0.25} - \text{Na}_{0.1} \text{ - apatite along in-plane : } \Gamma \\ = 0.62 \times 10^{12} \exp\left(\frac{-7.79}{kT}\right) (s^{-1}). \quad (13) $$

Several factors should be considered while estimating the hopping rates e.g., the transition rate will be the sum of the rates along the pathways that connect numerous paths (sub-channels) with adjacent nearby equivalent sites. For instance, the hopping rate between the ISZ sites can be determined via Eqs. ((8), (10), and (12)) with a multiplying factor of 6 in the case of diffusion along the z-axis, because six equivalent sub-channels connect the ISZ sites (Figs. 3).

There are three diffusion directions at 60 degrees along the in-plane diffusion pathway therefore, representing an isotropic nature in this plane. Eqs. ((9), (11), and (13)) will require a multiplying factor of 1.5 when estimating the diffusion coefficient along any in-plane coordinate (such as x or y). By keeping the above points in mind, one can directly

investigate the diffusion co-efficient for F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite via random walk KMC simulations.

### 3.4. Pb diffusivity (Kinetic Monte Carlo simulations)

The KMC technique is adopted in this study to simulate the three dimensional (3D) random-walk of Pb in F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite in order to calculate its effective activation energies and diffusion coefficients. Input data for such simulations includes co-ordinates of the insertion sites, the migration energies, and the attempt frequencies (Djimbi et al., 2015). KMC computes the averaged total travel time and mean square displacements (MSDs) of the Pb after one million simulation steps and at temperatures ranging from 600 K to 1500 K. The reliability of the obtained results is reflected by the linear relationship of time versus MSDs (Fig. 8) and show that Pb diffuses faster in apatite at high temperatures as compared to low temperatures. The averaged total travel time and the corresponding MSDs values are then used in Eqs. (3)–(6) and the diffusion co-efficients Dₓ, Dᵧ, D_z, and D_t were obtained at temperatures 600-1500 K. After simulating the diffusion coefficients (Dₓ, Dᵧ, D_z, and D_t), the data are modelled by using Arrhenius’ law, resulting in a linear relationship between ln(D_t) and 1/T with an intercept ln(D₀) and a slope Ea (Figs. 9a-c). These linear trends between ln(D_t) and 1/T gives effective activation energies (Ea) and pre-exponential factors (D₀) for F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite which can be represented either in the following way by using eq. 7 or in the form of plots as shown in Figs. 10-12:
<figure><img src="image_9.png" /><figcaption>Fig. 6. The activation-energy profiles of the Pb energetics in perfect Na₀.₁-apatite lattice. Each profile indicates the migration energy required for a Pb atom to move between neighbouring symmetrically equivalent minima. Three activation energy profiles are plotted here i.e. (a): is along z-axis or (001) direction, where Pb needs to overcome the energy barrier encountered at sites along z-axis, (b) is along (110) or in-plane directions, where Pb needs to cross the energy barrier encountered at sites along the orthogonal plane direction, and (c) shows the activation energy profile for the whole path that connect (110) direction with (001) indicating that if the Pb start diffusion first in the in-plane direction and when it encounter a site along the z-axis during the travel, the diffusion will shift again along the z-axis as can be seen in the energy profile trend.</figcaption></figure>

$$ \text{For Na}_{0.1} \text{ - apatite along z-axis : } \Gamma = 4.36 \times 10^{12} \exp\left(\frac{-2.8}{kT}\right) (s^{-1}), \quad (10) $$

$$ \text{For Na}_{0.1} \text{ - apatite along in-plane : } \Gamma = 3.51 \times 10^{12} \exp\left(\frac{-4.58}{kT}\right) (s^{-1}), \quad (11) $$

$$ \text{For Cl}_{0.25} - \text{Na}_{0.1} \text{ - apatite along z-axis : } \Gamma \\ = 3.51 \times 10^{12} \exp\left(\frac{-5.73}{kT}\right) (s^{-1}), \quad (12) $$

$$ \text{For Cl}_{0.25} - \text{Na}_{0.1} \text{ - apatite along in-plane : } \Gamma \\ = 0.62 \times 10^{12} \exp\left(\frac{-7.79}{kT}\right) (s^{-1}). \quad (13) $$

Several factors should be considered while estimating the hopping rates e.g., the transition rate will be the sum of the rates along the pathways that connect numerous paths (sub-channels) with adjacent nearby equivalent sites. For instance, the hopping rate between the ISZ sites can be determined via Eqs. ((8), (10), and (12)) with a multiplying factor of 6 in the case of diffusion along the z-axis, because six equivalent sub-channels connect the ISZ sites (Figs. 3).

There are three diffusion directions at 60 degrees along the in-plane diffusion pathway therefore, representing an isotropic nature in this plane. Eqs. ((9), (11), and (13)) will require a multiplying factor of 1.5 when estimating the diffusion coefficient along any in-plane coordinate (such as x or y). By keeping the above points in mind, one can directly

investigate the diffusion co-efficient for F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite via random walk KMC simulations.

### 3.4. Pb diffusivity (Kinetic Monte Carlo simulations)

The KMC technique is adopted in this study to simulate the three dimensional (3D) random-walk of Pb in F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite in order to calculate its effective activation energies and diffusion coefficients. Input data for such simulations includes co-ordinates of the insertion sites, the migration energies, and the attempt frequencies (Djimbi et al., 2015). KMC computes the averaged total travel time and mean square displacements (MSDs) of the Pb after one million simulation steps and at temperatures ranging from 600 K to 1500 K. The reliability of the obtained results is reflected by the linear relationship of time versus MSDs (Fig. 8) and show that Pb diffuses faster in apatite at high temperatures as compared to low temperatures. The averaged total travel time and the corresponding MSDs values are then used in Eqs. (3)–(6) and the diffusion co-efficients Dₓ, Dᵧ, D_z, and D_t were obtained at temperatures 600-1500 K. After simulating the diffusion coefficients (Dₓ, Dᵧ, D_z, and D_t), the data are modelled by using Arrhenius’ law, resulting in a linear relationship between ln(D_t) and 1/T with an intercept ln(D₀) and a slope Ea (Figs. 9a-c). These linear trends between ln(D_t) and 1/T gives effective activation energies (Ea) and pre-exponential factors (D₀) for F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite which can be represented either in the following way by using eq. 7 or in the form of plots as shown in Figs. 10-12:
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** The activation-energy profiles of the Pb energetics in perfect $Cl_{0.25}$-Na$_{0.1}$-apatite lattice. Each profile indicates the migration energy required for a Pb atom to move between neighbouring symmetrically equivalent minima. Three activation energy profiles are plotted here i.e. (a): is along z-axis or (001) direction, where Pb needs to overcome the energy barrier encountered at sites along z-axis, (b) is along (110) or in-plane directions, where Pb needs to cross the energy barrier encountered at sites along the orthogonal plane direction, and (c) shows the activation energy profile for the whole path that connect (110) direction with (001) indicating that if the Pb start diffusion first in the in-plane direction and when it encounter a site along the z-axis during the travel, the diffusion will shift again along the z-axis as can be seen in the energy profile trend.</figcaption></figure>

<caption>**Table 2**</caption>

computed parameters for F-apatite, Na$_{0.1}$-apatite, and $Cl_{0.25}$-Na$_{0.1}$-apatite. %Δ Ch_1991 and %Δ W_1985 represents the percentage difference between the computed Ea of this study and that of Cherniak et al. (1991) and Watson et al. (1985) respectively.

<table><thead><tr><th></th><th>Pb-Ca1<br>(Å)</th><th>Emig<br>(eV)</th><th>Ea (kJ/mol)</th><th>%Δ<br>Ch_1991</th><th>%Δ<br>W_1985</th><th>ν<br>(THz)</th></tr></thead><tbody><tr><td rowspan="2">F-apatite</td><td>z-axis</td><td>5.18</td><td>499.8</td><td>118</td><td>71</td><td>1.98</td></tr><tr><td>inplane</td><td>8.05</td><td>776.7</td><td>240</td><td>165</td><td>1.08</td></tr><tr><td rowspan="2">Na<sub>0.1</sub>-apatite</td><td>z-axis</td><td>2.80</td><td>270.2</td><td>18</td><td>8</td><td>4.36</td></tr><tr><td>inplane</td><td>4.58</td><td>441.9</td><td>93</td><td>51</td><td>3.51</td></tr><tr><td rowspan="2">Na<sub>0.1</sub>Cl<sub>0.25</sub>-<br>apatite</td><td>z-axis</td><td>5.73</td><td>552.9</td><td>142</td><td>89</td><td>3.51</td></tr><tr><td>inplane</td><td>7.79</td><td>751.6</td><td>229</td><td>157</td><td>0.62</td></tr><tr><td>Cherniak<br>et al.<br>(1991)</td><td></td><td>2.37</td><td>228.4</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Watson et al.,<br>1985</td><td></td><td>3.03</td><td>292.9</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

$$ \text{For F-apatite along z-axis: } D = 1.9 \times 10^{-1} \exp\left(-\frac{5.18}{kT}\right) (cm2s^{-1}), \quad (14) $$

$$ \text{For F-apatite along in-plane: } D = 2.8 \times 10^{-3} \exp\left(-\frac{8.04}{kT}\right) (cm2s^{-1}), \quad (15) $$

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** shows the mean square displacement (MSDs) of the diffusing Pb atom in Na-substituted apatite at 600 K and 1500 K obtained at 106 steps by using KMC simulation method.</figcaption></figure>

$$ \text{For Na}_{0.1} \text{-apatite along z-axis: } D = 3.9 \times 10^{-1} \exp\left(-\frac{2.8}{kT}\right) (cm2s^{-1}), \quad (16) $$
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** The activation-energy profiles of the Pb energetics in perfect $Cl_{0.25}$-Na$_{0.1}$-apatite lattice. Each profile indicates the migration energy required for a Pb atom to move between neighbouring symmetrically equivalent minima. Three activation energy profiles are plotted here i.e. (a): is along z-axis or (001) direction, where Pb needs to overcome the energy barrier encountered at sites along z-axis, (b) is along (110) or in-plane directions, where Pb needs to cross the energy barrier encountered at sites along the orthogonal plane direction, and (c) shows the activation energy profile for the whole path that connect (110) direction with (001) indicating that if the Pb start diffusion first in the in-plane direction and when it encounter a site along the z-axis during the travel, the diffusion will shift again along the z-axis as can be seen in the energy profile trend.</figcaption></figure>

<table><caption>**Table 2**

computed parameters for F-apatite, Na$_{0.1}$-apatite, and $Cl_{0.25}$-Na$_{0.1}$-apatite. %Δ Ch_1991 and %Δ W_1985 represents the percentage difference between the computed Ea of this study and that of Cherniak et al. (1991) and Watson et al. (1985) respectively.</caption><thead><tr><th></th><th>Pb-Ca1<br>(Å)</th><th>Emig<br>(eV)</th><th>Ea (kJ/mol)</th><th>%Δ<br>Ch_1991</th><th>%Δ<br>W_1985</th><th>ν<br>(THz)</th></tr></thead><tbody><tr><td rowspan="2">F-apatite</td><td>z-axis</td><td>5.18</td><td>499.8</td><td>118</td><td>71</td><td>1.98</td></tr><tr><td>inplane</td><td>8.05</td><td>776.7</td><td>240</td><td>165</td><td>1.08</td></tr><tr><td rowspan="2">Na<sub>0.1</sub>-apatite</td><td>z-axis</td><td>2.80</td><td>270.2</td><td>18</td><td>8</td><td>4.36</td></tr><tr><td>inplane</td><td>4.58</td><td>441.9</td><td>93</td><td>51</td><td>3.51</td></tr><tr><td rowspan="2">Na<sub>0.1</sub>Cl<sub>0.25</sub>-<br>apatite</td><td>z-axis</td><td>5.73</td><td>552.9</td><td>142</td><td>89</td><td>3.51</td></tr><tr><td>inplane</td><td>7.79</td><td>751.6</td><td>229</td><td>157</td><td>0.62</td></tr><tr><td>Cherniak<br>et al.<br>(1991)</td><td></td><td>2.37</td><td>228.4</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Watson et al.,<br>1985</td><td></td><td>3.03</td><td>292.9</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

$$ \text{For F-apatite along z-axis: } D = 1.9 \times 10^{-1} \exp\left(\frac{-5.18}{kT}\right) (cm2s^{-1}), \quad (14) $$

$$ \text{For F-apatite along in-plane: } D = 2.8 \times 10^{-3} \exp\left(\frac{-8.04}{kT}\right) (cm2s^{-1}), \quad (15) $$

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** shows the mean square displacement (MSDs) of the diffusing Pb atom in Na-substituted apatite at 600 K and 1500 K obtained at 106 steps by using KMC simulation method.</figcaption></figure>

$$ \text{For Na}_{0.1} \text{-apatite along z-axis: } D = 3.9 \times 10^{-1} \exp\left(\frac{-2.8}{kT}\right) (cm2s^{-1}), \quad (16) $$
<figure><img src="image_11.png" /><figcaption>(a) For Na<sub>0.1</sub>-apatite</figcaption></figure>

Ea = 4.57eV
D₀ = 3.02*10⁻³ cm² s⁻¹

<figure><img src="image_12.png" /><figcaption>(b) For F-apatite</figcaption></figure>

<figure><img src="image_13.png" /></figure>

Ea = 8.04eV
D₀ = 2.84*10⁻³ cm² s⁻¹

<figure><img src="image_14.png" /><figcaption>(c) For Cl<sub>0.25</sub>-Na<sub>0.1</sub>-apatite</figcaption></figure>

<figure><img src="image_15.png" /></figure>

Ea = 7.78
D₀ = 3.01*10⁻³ cm² s⁻¹

<figure><img src="image_16.png" /></figure>

Fig. 9. Estimation of the Pb diffusion coefficients in (a) Na<sub>0.1</sub>-apatite, (b) F-apatite, and (c) Cl<sub>0.25</sub>-Na<sub>0.1</sub>-apatite with respect to temperatures. The plot shows the relationship between ln(D) and 1/kBT according to the Arrhenius law.
<figure><img src="image_12.png" /><figcaption>(a) For Na<sub>0.1</sub>-apatite</figcaption></figure>

Ea = 4.57eV
D₀ = 3.02*10⁻³ cm² s⁻¹

<figure><img src="image_13.png" /><figcaption>(b) For F-apatite</figcaption></figure>

<figure><img src="image_14.png" /></figure>

Ea = 8.04eV
D₀ = 2.84*10⁻³ cm² s⁻¹

<figure><img src="image_15.png" /><figcaption>(c) For Cl<sub>0.25</sub>-Na<sub>0.1</sub>-apatite</figcaption></figure>

<figure><img src="image_16.png" /></figure>

Ea = 7.78
D₀ = 3.01*10⁻³ cm² s⁻¹

<figure><img src="image_17.png" /></figure>

Fig. 9. Estimation of the Pb diffusion coefficients in (a) Na<sub>0.1</sub>-apatite, (b) F-apatite, and (c) Cl<sub>0.25</sub>-Na<sub>0.1</sub>-apatite with respect to temperatures. The plot shows the relationship between ln(D) and 1/kBT according to the Arrhenius law.
<figure><img src="image_17.png" /><figcaption>Fig. 10. represents graphs showing the Arrhenius trends of diffusivity (cm2/s) against temperature (K) for a pure F-apatite supercell. The gray colour solid line indicates the diffusion profile along the z-axis direction while the red colour solid line represents the diffusion profile along (110) or in-plane direction computed in this work. The black and blue colours represent the diffusion profiles for the experiments of Watson et al. (1985) and Cherniak et al. (1991) respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 11. represents graphs showing the Arrhenius trends of diffusivity (cm2/s) against temperature (K) for a pure Na0.1-apatite supercell. The gray colour solid line indicates the diffusion profile along the z-axis direction while the red colour solid line represents the diffusion profile along (110) or in-plane direction computed in this work. The black and blue colours represent the diffusion profiles for the experiments of Watson et al. (1985) and Cherniak et al. (1991) respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 12. represents graphs showing the Arrhenius trends of diffusivity (cm2/s) against temperature (K) for a pure Na0.1Cl0.25-apatite supercell. The gray colour solid line indicates the diffusion profile along the z-axis direction while the red colour solid line represents the diffusion profile along (110) or in-plane direction computed in this work. The black and blue colours represent the diffusion profiles for the experiments of Watson et al. (1985) and Cherniak et al. (1991) respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$$ \text{For Na}_{0.1} - \text{apatite along in-plane : D} = 3.0 \times 10^{-3} \exp\left(\frac{-4.57}{kT}\right) (\text{cm}2\text{s}^{-1}), \qquad (17) $$

$$ \text{For Cl}_{0.25} - \text{Na}_{0.1} - \text{apatite along z-axis : D} \\ = 1.77 \times 10^{-1} \exp\left(\frac{-5.72}{kT}\right) (\text{cm}2\text{s}^{-1}), \qquad (18) $$

$$ \text{For Cl}_{0.25} - \text{Na}_{0.1} - \text{apatite along in-plane : D} \\ = 3.0 \times 10^{-3} \exp\left(\frac{-7.78}{kT}\right) (\text{cm}2\text{s}^{-1}), \qquad (19) $$

<figure><img src="image_20.png" /><figcaption>Fig. 13. shows the relative comparisons of the calculated closure temperatures (Tc) as a function of the effective grain size or effective-diffusion radius with the experimentally obtained Tc results of Watson et al. (1985) and Cherniak et al. (1991). Tc in these graphs is computed by using a spherical geometric factor (with A = 55) for all the perfect apatites considered in this work.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 10. represents graphs showing the Arrhenius trends of diffusivity (cm2/s) against temperature (K) for a pure F-apatite supercell. The gray colour solid line indicates the diffusion profile along the z-axis direction while the red colour solid line represents the diffusion profile along (110) or in-plane direction computed in this work. The black and blue colours represent the diffusion profiles for the experiments of Watson et al. (1985) and Cherniak et al. (1991) respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 11. represents graphs showing the Arrhenius trends of diffusivity (cm2/s) against temperature (K) for a pure Na<sub>0.1</sub>-apatite supercell. The gray colour solid line indicates the diffusion profile along the z-axis direction while the red colour solid line represents the diffusion profile along (110) or in-plane direction computed in this work. The black and blue colours represent the diffusion profiles for the experiments of Watson et al. (1985) and Cherniak et al. (1991) respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 12. represents graphs showing the Arrhenius trends of diffusivity (cm2/s) against temperature (K) for a pure Na<sub>0.1</sub>Cl<sub>0.25</sub>-apatite supercell. The gray colour solid line indicates the diffusion profile along the z-axis direction while the red colour solid line represents the diffusion profile along (110) or in-plane direction computed in this work. The black and blue colours represent the diffusion profiles for the experiments of Watson et al. (1985) and Cherniak et al. (1991) respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$$ \text{For Na}_{0.1} \text{ - apatite along in-plane : } D = 3.0 \times 10^{-3} \exp\left(\frac{-4.57}{kT}\right) (\text{cm2s}^{-1}), \qquad (17) $$

$$ \text{For Cl}_{0.25} \text{ - Na}_{0.1} \text{ - apatite along z-axis : } D \\ = 1.77 \times 10^{-1} \exp\left(\frac{-5.72}{kT}\right) (\text{cm2s}^{-1}), \qquad (18) $$

$$ \text{For Cl}_{0.25} \text{ - Na}_{0.1} \text{ - apatite along in-plane : } D \\ = 3.0 \times 10^{-3} \exp\left(\frac{-7.78}{kT}\right) (\text{cm2s}^{-1}), \qquad (19) $$

<figure><img src="image_21.png" /><figcaption>Fig. 13. shows the relative comparisons of the calculated closure temperatures (Tc) as a function of the effective grain size or effective-diffusion radius with the experimentally obtained Tc results of Watson et al. (1985) and Cherniak et al. (1991). Tc in these graphs is computed by using a spherical geometric factor (with A = 55) for all the perfect apatites considered in this work.</figcaption></figure>
Next step to this is the estimation of the closure temperature, Tc, at a cooling rate of 10 °C/Ma for F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite by using the parameters obtained from the KMC simulation that are modelled in Eqs. ((14)–(19)) and Figs. 10-12 and a spherical geometry with A = 55 (Dodson, 1973; Ganguly and Tirone, 1999). The estimated values of the closure temperature for F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite for different effective grain size is shown in Fig. 13. These results demonstrate that Pb retention in F-apatite, and Cl₀.₂₅-Na₀.₁-apatite is significantly higher than the Pb retention in Na₀.₁-apatite.

# 4. Discussions

## 4.1. Pb diffusion in defect-free F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite

As per results of Pb diffusion computed quantum mechanically for the super lattices of F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite, Pb diffuses primarily along the main z-axis in each of these lattices characterizing an anisotropic migration pathway. Due to the extremely large migration energy required to move Pb along the in-plane, the diffusion is very slow or minor along this direction at low temperatures, and therefore we are only left with the diffusion of Pb along the z-axes (also called the Florine-axis: Fig. 3). Along the in-plane (plane orthogonal to the z-axis) direction the Pb diffuses again to the z-axis path on its way when it moves towards the next cell. This is because, the direct jump from ISZ to ISIP has a high migration energy of 8.05 eV, making it energetically unfavourable as compared to the jump along z-axis (having migration energy of 5.18 eV). In this regard, rather than overcoming the high energy barrier of the direct in-plane direction (8.05 eV), the Pb atom favours to jump along the z-axis with an energetic effort of 5.18 eV for each jump. This study's KMC simulations produce an effective activation energy of 5.18 eV, which is the same as the migration energy along z-axis showing that the Pb diffusion along other directions is unfavourable. As a result, migrations between insertion sites along these directions could not be anticipated at temperatures below the closure temperatures. Therefore, the only direction in which the Pb diffusion may occur in each of the F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite is the z-axis. Also, several minerals, including zircon (Cherniak and Watson, 2001; Gautheron et al., 2020), exhibit such anisotropic diffusion too. Pb is retained in the perfect apatite structure of pure F-apatite, and Cl₀.₂₅-Na₀.₁-apatite even at higher temperatures in accordance with the computational results of this study except Na₀.₁-apatite for which the activation energy is 2.8 eV. Such Pb retention in both F-apatite, and Cl₀.₂₅-Na₀.₁-apatite demonstrate that a (U-Th)/Pb dating in such pure geochronometric systems can be used for any type of thermal history in the range of below ~1200 °C except Na₀.₁-apatite for which the closure temperature falls in the range of ~450-500 °C for an effective grain size of 100 µm. By increasing Na concentration results an increase in the Pb diffusion, but it still has enough closure temperature for a grain size of 100 µm (about 450-500 °C; Fig. 13) and is therefore can be used for intermediate range thermochronological applications. On the other hand, increasing Cl concentration (substitution in place of fluorine atom) decrease diffusion of Pb and enhance the closure temperature, make it more suitable for high temperature petrochronological problems. However, these findings show that the Na substitution has strong effect on the Pb diffusivity because even under low Na substitution constraints, the activation energy of Pb diffusion and closure temperature is considerably reduced (Fig. 14). For all the under-investigation apatite (F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite) in this study the Pb diffusion is strongly anisotropic and is characterized by low effective activation energy along the z-axis as compared to other directions.

These results are not in good agreement with the Pb diffusion results extracted experimentally in natural apatite (Watson et al., 1985; Cherniak et al., 1991; Figs. 9 and 11). However, the results obtained in Na

<figure><img src="image_21.png" /><figcaption>**Fig. 14.** represents graphs showing the Arrhenius trends of diffusivity (cm²/s) against temperature (K) for all the pure apatites supercells considered in this study. The gray colour solid line indicates the diffusion profile along the z-axis direction of apatite with 5%Na and no Cl, the red colour solid line represents the isotropic diffusion profile for 10%Na and 25%Cl substituted apatite, the black colour solid line denotes the anisotropic diffusion of Pb in 20% aligned Na substituted apatite, while green colour line is for 5%Na and 25%Cl substituted apatite, and at last the blue one is the isotropic diffusion line for the 20%Na and 25%Cl randomly substituted apatite. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

substitution rich apatite (i.e., along z-axis of Na₀.₁-apatite for instance) is quite close to the experimental results (Fig. 11) of the Durango apatite, the composition of which is Ca₉.₈₃ Na₀.₀₈ Sr₀.₀₁ REE₀.₀₉ (PO₄)₅.₈₇ (SO₄)₀.₀₅ (SiO₄)₀.₀₆ (AsO)₀.₀₁ (CO₃)₀.₀₁ F₁.₉₀ Cl₀.₁₂ (OH)₀.₀₁ (Cherniak et al., 1991). This is understandable given that the Pb diffusion data in natural apatite display complex behaviour, where differences in chemical compositions are evident. To model such highly variable crystal computationally is out of the scope of currently available computers therefore, only two elements (Na and Cl) with comparatively high concentration have been considered in this study. However, this study suggest that these Na-substituted regions of the apatite provides work-around pathways which probably enhance the diffusion in natural apatites. Additionally, it has already been noted that there is a significant discrepancy between results obtained using theoretical and experimental approaches (Djimbi et al., 2015; Gautheron et al., 2020; Reich et al., 2007; Saadoune et al., 2009).

## 4.2. Modelling the effect of chemical substitution on Pb diffusion at macroscopic scales

Pb diffusion data from naturally occurring apatite have been investigated experimentally by Watson et al. (1985) and Cherniak et al. (1991) but the effect of varying composition has not yet been studied in detail. Apatite's structure permits significant variations in composition, such as the partial or complete replacement of both the cationic and anionic sites, producing solid solutions and nonstoichiometric forms (Ptáček, 2016). As a result, the apatite structure can significantly accommodate almost half or more than half of the naturally occurring elements (Ptáček, 2016; Tamm and Peld, 2006). Due to their significant chemical variety and ion-exchange capability, apatite may be regarded as a geochronometer with variable diffusion behaviour. The theoretical method in this study provides a thorough understanding of the phenomenon of diffusion as a function of chemical substitutions or variation. According to previous discussions, diffusion along the z-axis and
Next step to this is the estimation of the closure temperature, Tc, at a cooling rate of 10 °C/Ma for F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite by using the parameters obtained from the KMC simulation that are modelled in Eqs. ((14)–(19)) and Figs. 10-12 and a spherical geometry with A = 55 (Dodson, 1973; Ganguly and Tirone, 1999). The estimated values of the closure temperature for F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite for different effective grain size is shown in Fig. 13. These results demonstrate that Pb retention in F-apatite, and Cl₀.₂₅-Na₀.₁-apatite is significantly higher than the Pb retention in Na₀.₁-apatite.

# 4. Discussions

## 4.1. Pb diffusion in defect-free F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite

As per results of Pb diffusion computed quantum mechanically for the super lattices of F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite, Pb diffuses primarily along the main z-axis in each of these lattices characterizing an anisotropic migration pathway. Due to the extremely large migration energy required to move Pb along the in-plane, the diffusion is very slow or minor along this direction at low temperatures, and therefore we are only left with the diffusion of Pb along the z-axes (also called the Florine-axis: Fig. 3). Along the in-plane (plane orthogonal to the z-axis) direction the Pb diffuses again to the z-axis path on its way when it moves towards the next cell. This is because, the direct jump from ISZ to ISIP has a high migration energy of 8.05 eV, making it energetically unfavourable as compared to the jump along z-axis (having migration energy of 5.18 eV). In this regard, rather than overcoming the high energy barrier of the direct in-plane direction (8.05 eV), the Pb atom favours to jump along the z-axis with an energetic effort of 5.18 eV for each jump. This study's KMC simulations produce an effective activation energy of 5.18 eV, which is the same as the migration energy along z-axis showing that the Pb diffusion along other directions is unfavourable. As a result, migrations between insertion sites along these directions could not be anticipated at temperatures below the closure temperatures. Therefore, the only direction in which the Pb diffusion may occur in each of the F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite is the z-axis. Also, several minerals, including zircon (Cherniak and Watson, 2001; Gautheron et al., 2020), exhibit such anisotropic diffusion too. Pb is retained in the perfect apatite structure of pure F-apatite, and Cl₀.₂₅-Na₀.₁-apatite even at higher temperatures in accordance with the computational results of this study except Na₀.₁-apatite for which the activation energy is 2.8 eV. Such Pb retention in both F-apatite, and Cl₀.₂₅-Na₀.₁-apatite demonstrate that a (U-Th)/Pb dating in such pure geochronometric systems can be used for any type of thermal history in the range of below ~1200 °C except Na₀.₁-apatite for which the closure temperature falls in the range of ~450-500 °C for an effective grain size of 100 µm. By increasing Na concentration results an increase in the Pb diffusion, but it still has enough closure temperature for a grain size of 100 µm (about 450-500 °C; Fig. 13) and is therefore can be used for intermediate range thermochronological applications. On the other hand, increasing Cl concentration (substitution in place of fluorine atom) decrease diffusion of Pb and enhance the closure temperature, make it more suitable for high temperature petrochronological problems. However, these findings show that the Na substitution has strong effect on the Pb diffusivity because even under low Na substitution constraints, the activation energy of Pb diffusion and closure temperature is considerably reduced (Fig. 14). For all the under-investigation apatite (F-apatite, Na₀.₁-apatite, and Cl₀.₂₅-Na₀.₁-apatite) in this study the Pb diffusion is strongly anisotropic and is characterized by low effective activation energy along the z-axis as compared to other directions.

These results are not in good agreement with the Pb diffusion results extracted experimentally in natural apatite (Watson et al., 1985; Cherniak et al., 1991; Figs. 9 and 11). However, the results obtained in Na

<figure><img src="image_22.png" /><figcaption>**Fig. 14.** represents graphs showing the Arrhenius trends of diffusivity (cm²/s) against temperature (K) for all the pure apatites supercells considered in this study. The gray colour solid line indicates the diffusion profile along the z-axis direction of apatite with 5%Na and no Cl, the red colour solid line represents the isotropic diffusion profile for 10%Na and 25%Cl substituted apatite, the black colour solid line denotes the anisotropic diffusion of Pb in 20% aligned Na substituted apatite, while green colour line is for 5%Na and 25%Cl substituted apatite, and at last the blue one is the isotropic diffusion line for the 20%Na and 25%Cl randomly substituted apatite. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

substitution rich apatite (i.e., along z-axis of Na₀.₁-apatite for instance) is quite close to the experimental results (Fig. 11) of the Durango apatite, the composition of which is Ca₉.₃₈ Na₀.₀₈ Sr₀.₀₁ REE₀.₀₉ (PO₄)₅.₈₇ (SO₄)₀.₀₅ (SiO₄)₀.₀₆ (AsO)₀.₀₁ (CO₃)₀.₀₁ F₁.₉₀ Cl₀.₁₂ (OH)₀.₀₁ (Cherniak et al., 1991). This is understandable given that the Pb diffusion data in natural apatite display complex behaviour, where differences in chemical compositions are evident. To model such highly variable crystal computationally is out of the scope of currently available computers therefore, only two elements (Na and Cl) with comparatively high concentration have been considered in this study. However, this study suggest that these Na-substituted regions of the apatite provides work-around pathways which probably enhance the diffusion in natural apatites. Additionally, it has already been noted that there is a significant discrepancy between results obtained using theoretical and experimental approaches (Djimbi et al., 2015; Gautheron et al., 2020; Reich et al., 2007; Saadoune et al., 2009).

## 4.2. Modelling the effect of chemical substitution on Pb diffusion at macroscopic scales

Pb diffusion data from naturally occurring apatite have been investigated experimentally by Watson et al. (1985) and Cherniak et al. (1991) but the effect of varying composition has not yet been studied in detail. Apatite's structure permits significant variations in composition, such as the partial or complete replacement of both the cationic and anionic sites, producing solid solutions and nonstoichiometric forms (Ptáček, 2016). As a result, the apatite structure can significantly accommodate almost half or more than half of the naturally occurring elements (Ptáček, 2016; Tamm and Peld, 2006). Due to their significant chemical variety and ion-exchange capability, apatite may be regarded as a geochronometer with variable diffusion behaviour. The theoretical method in this study provides a thorough understanding of the phenomenon of diffusion as a function of chemical substitutions or variation. According to previous discussions, diffusion along the z-axis and
along the orthogonal plane is characterized by high energy barriers and cannot occur at low to intermediate temperatures, suggesting the strong Pb retentivity in pure F-apatite. However, variety of chemical substitution or impurity may affect the Pb diffusion in apatite. Chlorine substitution along the fluorine axis (z-axis) slow down diffusion by blocking Pb. On the other hand the substitution of Na in place of Ca(2) enhance the diffusion both along the z-axis and along the orthogonal plane. Therefore, in all the three cases the Cl-substitution can be a valid candidate that increase the closure temperature (Fig. 13) and so can be used to study high temperature petrochronological problems.

Following is a summary of the key findings of the theoretical analysis, validated by the available data: (1) Pb diffusion is anisotropic and is dominantly occurring along the z-axis. The diffusion is slow in a perfect crystal that has no defects and chemical impurities, and so Pb is hardly released. (2) Variety of chemical substitutions (I.e., Na and Cl in this particular case) in the crystal effect the diffusion. The diffusion along both the in-plane and z-axis become faster in case of Na substitution while Cl blocks the diffusion since it would be too energy-intensive to hop along both these directions. So, the widening and opening of pathways by some chemical means (which in this case is Na) at high temperatures play a critical role in enhancement of Pb diffusion along the z-axis direction. However, this process is not true for every substituent as replacement of Cl for F in apatite causes a Pb atom to hops aside in the orthogonal plane after being colliding with a Cl atom during its way along the z-axis. (3) Any chemical replacement discussed in Tamm and Peld (2006) and Ptáček (2016) can either help in releasing or in the blocking of Pb atom but at any case studied here Pb will at least remain there at low to intermediate temperatures ~400 °C and above. However, the model's intent is to measure this diffusion impact as shown in Fig. 13. The activation energy and Pb loss rates of this study correspond with some of the chemical changes of natural samples (Cherniak et al., 1991; Ptáček, 2016), in accordance to the experimental data shown in Figs. 10-12. Any of the six channels present along the z-axis in pure F-apatite releases the Pb it receives from the decays at very high temperature (Fig. 13), according to (1). This is changed due to the introduction of impurities in the crystal lattice creation, which contribute to some diffusion by mechanisms 2) and 3). (4) In case, if the paths for the workaround are open, the specifics of the Cl distribution are not important, unless Cl atoms are specifically aligned along the z-axis. The results for both aligned or periodic (green) and randomly (blue) distributed Cl together with Na are shown in Fig. 14. Here, the shifts between the lines serve as an illustration of this. (5) The presence of Cl-atoms practically makes diffusion isotropic. This phenomenon can be explained by the fact that in normal cases of pure F- and Na-apatite the diffusion along the z-axis is preferred over the in-plane diffusion and happens along the F-atoms. But, when a Cl-atom is encountered in the path of diffusion along z-axis, the diffusion of Pb prefer the plane orthogonal to z-axis i.e., the in-plane diffusion. Once the blocking point is covered via the in-plane diffusion, it again restricts the Pb diffusion through the F-atoms along the z-axis and therefore imposes a diffusion coefficient for this direction that is nearly equal to that of in-plane diffusion (combination of in-plane and F-atoms direction) as shown in Fig. 14.

The scenario of the Cl blocking effect is comparable to the study of He diffusion in apatite by Djimbi et al. (2015), however is different when it comes to Na substitution where it significantly overcome the energy barriers and make it much easier for Pb to diffuse as compared to pure F-apatite, and $Cl_{0.25}-Na_{0.1}$-apatite. The aforementioned discussion show that Na-atoms have a significant impact on diffusion both along the z-axis and along the in-plane direction, depending on the site of substitution whether it have replaced the Ca(1) or the Ca(2) sites. In case of Ca (1) it may only effect the in-plane diffusion due to its location in the apatite's crystal lattice however, it is not the subject of this work. Strength of the Pb diffusion in both Na and Cl substituted apatite expected to be depends on their respective concentrations. In this study we were tried to build the crystal that nearly resemble the natural one that

is studied in Cherniak et al. (1991). Nonetheless, the chlorine atom's influence is only anticipated to have an impact on the energy landscape in its close vicinities inside an apatite crystal. Once the energy profiles or barriers and normal-mode frequencies are computed at the initial sites and saddle-points of F-apatite, 10% and/or 25% Na- and/or Cl-substituted apatites respectively. The KMC simulations are then used to calculate the diffusivity for lower/higher Na and/or Cl-concentrations with barriers and frequencies based on the local presence of respective Na-, Cl- or F-atoms concentrations. As has already been noted, the concentration of sodium and chlorine considerably effect both the in-plane as well z-axis diffusion and also our aim in this work is to project this effect close to the chemistry/heterogeneity of the natural apatite sample of Cherniak et al. (1991). Therefore, the effective Pb diffusions has been determined for Na- and Cl-concentration values range from 5 to 20% and 25-50% respectively as well as for aligned and random distribution of both Na- and Cl-atoms as shown in Fig. 14. F-apatite with 5% (gray colour line in Fig. 14) and 20% Na-concentration (aligned; black colour line in Fig. 14) and $Cl_{0.25}-Na_{0.05}$-apatite (25% Cl and 5% Na; green colour line in Fig. 14) exhibit an anisotropic diffusion along the z-axis while apatite with 10-20% Na and 25-50% Cl (red and blue colour lines in Fig. 14) randomly distributed experience a Pb diffusion of isotropic nature. As with both F-apatite and $Cl_{0.25}$-apatite, the diffusivity asymptotically approaches the in-plane diffusivity with increasing the Na-concentrations and so the Pb diffusion turns more and more isotropic when the level of Na-concentration is increased. This implies that the Pb-diffusion is strongly heterogenous depending on the apatite's chemistry.

### 4.3. Investigating the effect on the closure temperatures

For geological applications, the closure temperature (Tc) is a synthetic quantity that quantifies the temperature dependence of diffusion in a more understandable manner than $D_0$ and $E_a$ (e.g., Djimbi et al., 2015; Reiners et al., 2004; Reiners and Brandon, 2006). For a specific temperature evolution for instance, a cooling rate of 10 °C/Ma, it is defined as the temperature related to the known thermochronologic age of a system (Dodson, 1973; Ganguly and Tirone, 1999). The closure temperatures for different composition apatites are estimated using the methodology of Dodson (1973) and Cherniak (2010) and considering the geometric parameter equal to 55 (for sphere; Dodson, 1973; Ganguly and Tirone, 1999). To solve the diffusion equation requires the diffusivities that are estimated in the KMC simulations of this study. After acquiring the results of Tc, the outputs are plotted with respect to the effective grain size of the apatite and the obtained results were compared to the available experimental results of Cherniak et al. (1991) and Watson et al. (1985) for Durango apatite (Fig. 13). All the diffusivities have been calculated for crystals with only defects introduced due to the incorporation of Pb. Apart from that, no additional defects are taken into account, and the results are obtained for perfect Na- and Cl-substituted crystals. In comparison to pure F-apatite the closure temperature of $Na_{0.1}$-substituted apatite is minimal while that of $Na_{0.1}-Cl_{0.25}$-substituted apatite is higher than that of pure F-apatite for different effective grain sizes (Fig. 13). The decrease in diffusivity along the z-axis is caused by the blocking action of Cl-atoms in Na-substituted apatite. This blocking effect force the Pb atom to diffuse through workaround paths along the in-plane direction and as a result Tc are increased as compared to $Na_{0.1}$-apatite and F-apatite. Our results differ from the Tc predicted with the empirically derived diffusion coefficient from Cherniak et al. (1991) and Watson et al. (1985). However, when heterogeneity is introduced (specially Na), then they fall within the close range of the natural Durango apatite of both Cherniak et al. (1991) and Watson et al. (1985). Fig. 10-14 shows how the Na- and Cl-content has a considerable impact on the diffusion and closure temperature of apatite. However, no study is available that figure out the impact of the chemical composition of apatite on the Pb-diffusion co-efficient, but some of the work like Miro et al. (2006), Djimbi et al. (2015), and Bassal et al. (2022) have explained the effect of chemical composition on the diffusion co-
along the orthogonal plane is characterized by high energy barriers and cannot occur at low to intermediate temperatures, suggesting the strong Pb retentivity in pure F-apatite. However, variety of chemical substitution or impurity may affect the Pb diffusion in apatite. Chlorine substitution along the fluorine axis (z-axis) slow down diffusion by blocking Pb. On the other hand the substitution of Na in place of Ca(2) enhance the diffusion both along the z-axis and along the orthogonal plane. Therefore, in all the three cases the Cl-substitution can be a valid candidate that increase the closure temperature (Fig. 13) and so can be used to study high temperature petrochronological problems.

Following is a summary of the key findings of the theoretical analysis, validated by the available data: (1) Pb diffusion is anisotropic and is dominantly occurring along the z-axis. The diffusion is slow in a perfect crystal that has no defects and chemical impurities, and so Pb is hardly released. (2) Variety of chemical substitutions (i.e., Na and Cl in this particular case) in the crystal effect the diffusion. The diffusion along both the in-plane and z-axis become faster in case of Na substitution while Cl blocks the diffusion since it would be too energy-intensive to hop along both these directions. So, the widening and opening of pathways by some chemical means (which in this case is Na) at high temperatures play a critical role in enhancement of Pb diffusion along the z-axis direction. However, this process is not true for every substituent as replacement of Cl for F in apatite causes a Pb atom to hops aside in the orthogonal plane after being colliding with a Cl atom during its way along the z-axis. (3) Any chemical replacement discussed in Tamm and Peld (2006) and Ptáček (2016) can either help in releasing or in the blocking of Pb atom but at any case studied here Pb will at least remain there at low to intermediate temperatures ~400 °C and above. However, the model's intent is to measure this diffusion impact as shown in Fig. 13. The activation energy and Pb loss rates of this study correspond with some of the chemical changes of natural samples (Cherniak et al., 1991; Ptáček, 2016), in accordance to the experimental data shown in Figs. 10-12. Any of the six channels present along the z-axis in pure F-apatite releases the Pb it receives from the decays at very high temperature (Fig. 13), according to (1). This is changed due to the introduction of impurities in the crystal lattice creation, which contribute to some diffusion by mechanisms 2) and 3). (4) In case, if the paths for the workaround are open, the specifics of the Cl distribution are not important, unless Cl atoms are specifically aligned along the z-axis. The results for both aligned or periodic (green) and randomly (blue) distributed Cl together with Na are shown in Fig. 14. Here, the shifts between the lines serve as an illustration of this. (5) The presence of Cl-atoms practically makes diffusion isotropic. This phenomenon can be explained by the fact that in normal cases of pure F- and Na-apatite the diffusion along the z-axis is preferred over the in-plane diffusion and happens along the F-atoms. But, when a Cl-atom is encountered in the path of diffusion along z-axis, the diffusion of Pb prefer the plane orthogonal to z-axis i.e., the in-plane diffusion. Once the blocking point is covered via the in-plane diffusion, it again restricts the Pb diffusion through the F-atoms along the z-axis and therefore imposes a diffusion coefficient for this direction that is nearly equal to that of in-plane diffusion (combination of in-plane and F-atoms direction) as shown in Fig. 14.

The scenario of the Cl blocking effect is comparable to the study of He diffusion in apatite by Djimbi et al. (2015), however is different when it comes to Na substitution where it significantly overcome the energy barriers and make it much easier for Pb to diffuse as compared to pure F-apatite, and $Cl_{0.25}-Na_{0.1}$-apatite. The aforementioned discussion show that Na-atoms have a significant impact on diffusion both along the z-axis and along the in-plane direction, depending on the site of substitution whether it have replaced the Ca(1) or the Ca(2) sites. In case of Ca(1) it may only effect the in-plane diffusion due to its location in the apatite's crystal lattice however, it is not the subject of this work. Strength of the Pb diffusion in both Na and Cl substituted apatite expected to be depends on their respective concentrations. In this study we were tried to build the crystal that nearly resemble the natural one that

is studied in Cherniak et al. (1991). Nonetheless, the chlorine atom's influence is only anticipated to have an impact on the energy landscape in its close vicinities inside an apatite crystal. Once the energy profiles or barriers and normal-mode frequencies are computed at the initial sites and saddle-points of F-apatite, 10% and/or 25% Na- and/or Cl-substituted apatites respectively. The KMC simulations are then used to calculate the diffusivity for lower/higher Na and/or Cl-concentrations with barriers and frequencies based on the local presence of respective Na-, Cl- or F-atoms concentrations. As has already been noted, the concentration of sodium and chlorine considerably effect both the in-plane as well z-axis diffusion and also our aim in this work is to project this effect close to the chemistry/heterogeneity of the natural apatite sample of Cherniak et al. (1991). Therefore, the effective Pb diffusions has been determined for Na- and Cl-concentration values range from 5 to 20% and 25-50% respectively as well as for aligned and random distribution of both Na- and Cl-atoms as shown in Fig. 14. F-apatite with 5% (gray colour line in Fig. 14) and 20% Na-concentration (aligned; black colour line in Fig. 14) and $Cl_{0.25}-Na_{0.05}$-apatite (25% Cl and 5% Na; green colour line in Fig. 14) exhibit an anisotropic diffusion along the z-axis while apatite with 10-20% Na and 25-50% Cl (red and blue colour lines in Fig. 14) randomly distributed experience a Pb diffusion of isotropic nature. As with both F-apatite and $Cl_{0.25}$-apatite, the diffusivity asymptotically approaches the in-plane diffusivity with increasing the Na-concentrations and so the Pb diffusion turns more and more isotropic when the level of Na-concentration is increased. This implies that the Pb-diffusion is strongly heterogenous depending on the apatite's chemistry.

### 4.3. Investigating the effect on the closure temperatures

For geological applications, the closure temperature (Tc) is a synthetic quantity that quantifies the temperature dependence of diffusion in a more understandable manner than $D_0$ and $E_a$ (e.g., Djimbi et al., 2015; Reiners et al., 2004; Reiners and Brandon, 2006). For a specific temperature evolution for instance, a cooling rate of 10 °C/Ma, it is defined as the temperature related to the known thermochronologic age of a system (Dodson, 1973; Ganguly and Tirone, 1999). The closure temperatures for different composition apatites are estimated using the methodology of Dodson (1973) and Cherniak (2010) and considering the geometric parameter equal to 55 (for sphere; Dodson, 1973; Ganguly and Tirone, 1999). To solve the diffusion equation requires the diffusivities that are estimated in the KMC simulations of this study. After acquiring the results of Tc, the outputs are plotted with respect to the effective grain size of the apatite and the obtained results were compared to the available experimental results of Cherniak et al. (1991) and Watson et al. (1985) for Durango apatite (Fig. 13). All the diffusivities have been calculated for crystals with only defects introduced due to the incorporation of Pb. Apart from that, no additional defects are taken into account, and the results are obtained for perfect Na- and Cl-substituted crystals. In comparison to pure F-apatite the closure temperature of $Na_{0.1}$-substituted apatite is minimal while that of $Na_{0.1}-Cl_{0.25}$-substituted apatite is higher than that of pure F-apatite for different effective grain sizes (Fig. 13). The decrease in diffusivity along the z-axis is caused by the blocking action of Cl-atoms in Na-substituted apatite. This blocking effect force the Pb atom to diffuse through workaround paths along the in-plane direction and as a result Tc are increased as compared to $Na_{0.1}$-apatite and F-apatite. Our results differ from the Tc predicted with the empirically derived diffusion coefficient from Cherniak et al. (1991) and Watson et al. (1985). However, when heterogeneity is introduced (specially Na), then they fall within the close range of the natural Durango apatite of both Cherniak et al. (1991) and Watson et al. (1985). Fig. 10-14 shows how the Na- and Cl-content has a considerable impact on the diffusion and closure temperature of apatite. However, no study is available that figure out the impact of the chemical composition of apatite on the Pb-diffusion co-efficient, but some of the work like Miro et al. (2006), Djimbi et al. (2015), and Bassal et al. (2022) have explained the effect of chemical composition on the diffusion co-
efficient of He in various minerals. The study of Miro et al. (2006) and Djimbi et al. (2015) reports the same concept about increasing the Cl-concentration in apatite leads to block the diffusion and as a result give high closure temperature. Even though this study does not measure the He-diffusion coefficient but instead gives Pb-diffusion coefficients for varying apatite chemistry, however, still the concept of high closure temperature explained in their work appears to be consistent with our observations. In addition to the effect of composition, this work also examined the impact of anisotropy on the closure temperature as shown in Figs. 10-14.

Diffusion is important as it has the most vital role in the closure temperature of Pb in apatite and other minerals. Also, to date a particular mineral specially using U-Th/Pb geochronometric system, it is typically precise enough to know that the mineral had not experienced temperatures above its Tc as it defines the temperature range/zone within which the diffusing atoms are retained inside the crystal of that mineral (Cherniak et al., 1991; Cherniak and Watson, 2001). The closure temperatures computed using a geometric factor of 55 (for sphere) and based on the anisotropy is ~1300 °C for Cl₀.₂₅⁻ Na₀.₁-apatite, ~1100 °C for F-apatite and 450-500 °C for Na₀.₁-apatite at an effective grain size of 60 µm. Fig. 13 demonstrates how the effective radius affects the closure

temperature in different compositional apatites. The closure temperature (Fig. 13) reflects very high to intermediate (1300-500 °C) values depending on the chemical composition of apatite and therefore shows that the (U—Th)/Pb dating system of apatite can be used as a geochronometer to investigate intermediate to high temperature geological processes. In addition to the above, even a small crystal sizes (10 µm; Fig. 13) of Na₀.₁-apatite, F-apatite, and Cl₀.₂₅⁻ Na₀.₁-apatite remained at temperatures below 500 °C, 1100 °C, and 1300 °C (i.e., conditions in lower crust and mantle), indicates Pb retention, which is an important entity in the accuracy of (U—Th)/Pb geochronology.

### 4.4. Petrochronological implications to the natural apatite samples

To better understand the significance of U–Th–Pb isotopic dating for apatite and other minerals in deciphering crustal events, it is essential to examine the kinetic processes governing these elements within mineral structures and the factors influencing them. Such examination provides valuable insights into thermal histories (Cherniak and Watson, 2001; Zhang et al., 2021). Apatite is widely distributed in igneous, metamorphic, and sedimentary rocks, playing a notable role in the field of geochronology. Understanding the atomic structure and varying

<figure><img src="image_22.png" /><figcaption>**Fig. 15.** Na concentrations in natural apatites (a) from the Ilimaussaq intrusion in South Greenland (Rønsbo, 1989) and Clinohydroxylapatite from altered leucogabbro in northwestern Ontario, Canada (Chakhmouradian and Medici, 2006; inside blue rectangle), (b) from carbonatite melt inclusions in apatites given in GEOROC dataset, https://georoc.eu; DIGIS Team (2022), and (c) from Apatite-group minerals in Pilansberg alkaline complex nepheline syenite, South Africa (Liferovich and Mitchell, 2006 and references therein). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
efficient of He in various minerals. The study of Miro et al. (2006) and Djimbi et al. (2015) reports the same concept about increasing the Cl-concentration in apatite leads to block the diffusion and as a result give high closure temperature. Even though this study does not measure the He-diffusion coefficient but instead gives Pb-diffusion coefficients for varying apatite chemistry, however, still the concept of high closure temperature explained in their work appears to be consistent with our observations. In addition to the effect of composition, this work also examined the impact of anisotropy on the closure temperature as shown in Figs. 10-14.

Diffusion is important as it has the most vital role in the closure temperature of Pb in apatite and other minerals. Also, to date a particular mineral specially using U-Th/Pb geochronometric system, it is typically precise enough to know that the mineral had not experienced temperatures above its Tc as it defines the temperature range/zone within which the diffusing atoms are retained inside the crystal of that mineral (Cherniak et al., 1991; Cherniak and Watson, 2001). The closure temperatures computed using a geometric factor of 55 (for sphere) and based on the anisotropy is ~1300 °C for Cl₀.₂₅⁻ Na₀.₁-apatite, ~1100 °C for F-apatite and 450-500 °C for Na₀.₁-apatite at an effective grain size of 60 µm. Fig. 13 demonstrates how the effective radius affects the closure

temperature in different compositional apatites. The closure temperature (Fig. 13) reflects very high to intermediate (1300-500 °C) values depending on the chemical composition of apatite and therefore shows that the (U—Th)/Pb dating system of apatite can be used as a geochronometer to investigate intermediate to high temperature geological processes. In addition to the above, even a small crystal sizes (10 µm; Fig. 13) of Na₀.₁-apatite, F-apatite, and Cl₀.₂₅⁻ Na₀.₁-apatite remained at temperatures below 500 °C, 1100 °C, and 1300 °C (i.e., conditions in lower crust and mantle), indicates Pb retention, which is an important entity in the accuracy of (U—Th)/Pb geochronology.

### 4.4. Petrochronological implications to the natural apatite samples

To better understand the significance of U–Th–Pb isotopic dating for apatite and other minerals in deciphering crustal events, it is essential to examine the kinetic processes governing these elements within mineral structures and the factors influencing them. Such examination provides valuable insights into thermal histories (Cherniak and Watson, 2001; Zhang et al., 2021). Apatite is widely distributed in igneous, metamorphic, and sedimentary rocks, playing a notable role in the field of geochronology. Understanding the atomic structure and varying

<figure><img src="image_23.png" /><figcaption>**Fig. 15.** Na concentrations in natural apatites (a) from the Ilimaussaq intrusion in South Greenland (Rønsbo, 1989) and Clinohydroxylapatite from altered leucogabbro in northwestern Ontario, Canada (Chakhmouradian and Medici, 2006; inside blue rectangle), (b) from carbonatite melt inclusions in apatites given in GEOROC dataset, https://georoc.eu; DIGIS Team (2022), and (c) from Apatite-group minerals in Pilansberg alkaline complex nepheline syenite, South Africa (Liferovich and Mitchell, 2006 and references therein). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
chemical composition of apatite is essential to grasp the necessity of
knowing the impact of its structure, composition, and crystal chemistry
on petrochronology. To assess the relevance of this work to the
geochronology of natural apatites, it begins by modelling (using com-
bined DFT + TST + KMC computations) the composition close to
Durango apatite, which is used as a standard in some LA-ICP-MS analysis
of apatites (Chew et al., 2016). The composition is further simulated to
the Na (0%-20% apfu of Ca) and Cl (0-25% apfu of F) concentrations
observed in some of the natural apatite samples discussed in Rønsbo
(1989); Fig. 15a), Chakhmouradian and Medici (2006); including the
references therein; Fig. 15a), Liferovich and Mitchell (2006); including
the references therein; Fig. 15), and those with carbonatite melt in-
clusions as listed in the datasets of Geochemistry of Rocks of the Oceans
and Continents datasets (GEOROC, https://georoc.eu; DIGIS Team,
2022; Fig. 15b), which are predominantly rich in Na₂O as compared to
other melts inclusions. This study's wide range of Tc variability in these
modelled apatites suggest that chemical composition is an important
consideration when using the (U—Th)/Pb geochronometric system.
These results illustrate that (U—Th)/Pb geochronology of apatite may
be an effective petrochronological technique to predict the age and
thermal history of intermediate to high temperature rocks on Earth and
other terrestrial planets as also explained by the experimental studies of
Schoene and Bowring (2007) and Blackburn et al. (2011, 2012). Rønsbo
(1989) reported Na-rich apatites from quartz-bearing peralkaline
pegmatitic rocks and based on these Na-substituted apatites and asso-
ciated mineral assemblages it was concluded that such apatites crys-
tallize in highly peralkaline geologic environments. Chakhmouradian
and Medici (2006) also reported Na-rich apatites with excessive sulfer in
their crystal structure formed by metasomatic reaction of sulfate-
bearing fluid and plagioclase, and that Na solubility occurs even at
low crystallization temperatures (further details are in Chakhmouradian
and Medici, 2006). Similarly, the natural Sr-Na-REE-rich apatites
described by Liferovich and Mitchell (2006; and references therein)
originate in the post-metasomatic (natrolite-rich hydrothermal) alter-
ation of the host nepheline syenite of the Pilansberg alkaline complex,
South Africa, and are associated with natrolite and formed at the
reduced alkalinity stage. Furthermore, carbonatitic melt may inject Na-
oxides into apatites (calciocarbonatite co-precipitated) crystallized from
melellite-nephelenite-derived carbonatitic melt (Guzmics et al., 2011,
2015). As a result, the chemistry of the spots used for U-Th/Pb analysis
must be thoroughly investigated because it affects Pb diffusion and thus
petrochronology.

Based on the Na concentrations reported by Rønsbo (1989), Chakhmouradian and Medici (2006), and Liferovich and Mitchell (2006; and references therein) and that modelled in Fig. 14, the dotted black and cyan colour lines in Fig. 15a&c marks the threshold of Na-concentration > 10% apfu of Ca and > 5% apfu of Ca respectively that has a significant impact on the apatite's U-Th/Pb geochronology of high temperature crustal events. As the closure temperature varies with chemical composition (i.e., Na and Cl) therefore, the classification of these apatites based on chemical composition is significant before using it for U-Th/Pb petrochronology. It should be noted that this study only models the effect of Na and Cl, and because natural apatites vary greatly in composition, the effect of other elements on U-Th/Pb geochronology hosted by these natural apatites needs to be investigated further using the same methodology. Because Pb diffusion occurs in apatites over micron length scales at temperatures relevant to the formation of the middle and lower crusts, thermal history information of the middle to lower crust and even at the crust-mantle boundary can be derived from U-Pb date profiles preserved within these apatites at Tc (500-1300 °C) depending on the chemical composition of the apatite. Also, apatites with high Tc provide more reliable ages than those with low Tc. (Krogstad and Walker, 1994). Further, based on the chemistry of the apatite, it can be determined whether it is useful for obtaining intermediate or high temperature petrochronology information.

5. Conclusions

A multi-scale approach was adopted in this study to demonstrate the effect of the chemical composition on the diffusion of Pb in apatite. The computed Pb insertion energies shows that it can be inserted with an energy of -2.03 eV, -2.21 eV, and -1.94 eV in pure F-apatite, Na₀.₁-apatite, and Cl₀.₂₅- Na₀.₁-apatite respectively. The diffusion of Pb in all these investigated apatites take places along the two preferential pathways i.e., one along the z-axis and the other along the plane orthogonal to z-axis which is termed as the in-plane direction. By using the KMC approach and the obtained quantum mechanically computed data, a 3D random walk simulation is performed to see the effect of variable chemical composition on the diffusion of Pb in apatite at the macro-scale/mineral scale level. The diffusion coefficients for the Pb atom at various temperatures are determined along each path and pre-exponential factor and Ea were extracted from MSDs using Arrhenius plot. Based on the Pb diffusion profiles estimated in different apatite lattices with different composition, one can conclude that Pb diffusion changes with chemical substitutions. These substitutions contribute towards the isotropic nature of Pb diffusion in apatite and conclude that the Pb-diffusion is strongly heterogenous depending on the apatite's chemistry. This work studies the Pb diffusion in terms of geochronological application and so the integration of all the findings leads us to a conclusion that the (U-Th)/Pb geochronology of apatite may be useful for predicting the age and geological processes of intermediate to high temperature thermal events in upper and lower crust of Earth and other terrestrial planets depending on the chemical composition of these apatites. Moreover, in order to improve the accuracy of apatite U-Th/Pb geochronology, a thorough understanding of apatite chemistry is required, i.e., the effect of other elements in the various natural groups of apatites must be investigated in the context of the Pb closure temperature. It is necessary because different elements in apatites have different effects on the apatite closure temperature, limiting its use in the investigation of crustal events below their respective Pb retentivity (or closure) temperature. This study also emphasizes that the findings regarding the diffusivity of partly Na- and/or Cl-substituted apatites infers about the importance of the adopted methodology in the characterization and computational modelling of special geological systems with spatial atomic arrangements, incorporating heterogeneities in crystals, and mobility of inserted atoms and their respective energetics profiles via integrated approach of KMC, DFT, NEB, and TST simulations.

Declaration of Competing Interest

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.

Acknowledgments

This research work is published from the PhD work of the first
author, funded by the Higher Education Commission of Pakistan under
the Faculty Development Program for the Khushal Khan Khattak Uni-
versity, Karak, Pakistan with grant No. 386-388/PC/PAB/ KKKUK/21.
Partial funds are provided by the University of Vienna, Austria. The
Department of Lithospheric Research, University of Vienna is also
acknowledged for providing computational facilities and access to the
online library. Thanks are extended to the editors and anonymous re-
viewers for providing feedback to improve the manuscript.

References

Andersen, M., Panosetti, C., Reuter, K., 2019. A practical guide to surface kinetic Monte Carlo simulations. *Front. Chem.* 7, 202.
chemical composition of apatite is essential to grasp the necessity of
knowing the impact of its structure, composition, and crystal chemistry
on petrochronology. To assess the relevance of this work to the
geochronology of natural apatites, it begins by modelling (using com-
bined DFT + TST + KMC computations) the composition close to
Durango apatite, which is used as a standard in some LA-ICP-MS analysis
of apatites (Chew et al., 2016). The composition is further simulated to
the Na (0%-20% apfu of Ca) and Cl (0-25% apfu of F) concentrations
observed in some of the natural apatite samples discussed in Rønsbo
(1989); Fig. 15a), Chakhmouradian and Medici (2006); including the
references therein; Fig. 15a), Liferovich and Mitchell (2006); including
the references therein; Fig. 15), and those with carbonatite melt in-
clusions as listed in the datasets of Geochemistry of Rocks of the Oceans
and Continents datasets (GEOROC, https://georoc.eu; DIGIS Team,
2022; Fig. 15b), which are predominantly rich in Na₂O as compared to
other melts inclusions. This study's wide range of Tc variability in these
modelled apatites suggest that chemical composition is an important
consideration when using the (U—Th)/Pb geochronometric system.
These results illustrate that (U—Th)/Pb geochronology of apatite may
be an effective petrochronological technique to predict the age and
thermal history of intermediate to high temperature rocks on Earth and
other terrestrial planets as also explained by the experimental studies of
Schoene and Bowring (2007) and Blackburn et al. (2011, 2012). Rønsbo
(1989) reported Na-rich apatites from quartz-bearing peralkaline
pegmatitic rocks and based on these Na-substituted apatites and asso-
ciated mineral assemblages it was concluded that such apatites crys-
tallize in highly peralkaline geologic environments. Chakhmouradian
and Medici (2006) also reported Na-rich apatites with excessive sulfer in
their crystal structure formed by metasomatic reaction of sulfate-
bearing fluid and plagioclase, and that Na solubility occurs even at
low crystallization temperatures (further details are in Chakhmouradian
and Medici, 2006). Similarly, the natural Sr-Na-REE-rich apatites
described by Liferovich and Mitchell (2006; and references therein)
originate in the post-metasomatic (natrolite-rich hydrothermal) alter-
ation of the host nepheline syenite of the Pilansberg alkaline complex,
South Africa, and are associated with natrolite and formed at the
reduced alkalinity stage. Furthermore, carbonatitic melt may inject Na-
oxides into apatites (calciocarbonatite co-precipitated) crystallized from
melellite-nephelenite-derived carbonatitic melt (Guzmics et al., 2011,
2015). As a result, the chemistry of the spots used for U-Th/Pb analysis
must be thoroughly investigated because it affects Pb diffusion and thus
petrochronology.

Based on the Na concentrations reported by Rønsbo (1989), Chakhmouradian and Medici (2006), and Liferovich and Mitchell (2006; and references therein) and that modelled in Fig. 14, the dotted black and cyan colour lines in Fig. 15a&c marks the threshold of Na-concentration > 10% apfu of Ca and > 5% apfu of Ca respectively that has a significant impact on the apatite's U-Th/Pb geochronology of high temperature crustal events. As the closure temperature varies with chemical composition (i.e., Na and Cl) therefore, the classification of these apatites based on chemical composition is significant before using it for U-Th/Pb petrochronology. It should be noted that this study only models the effect of Na and Cl, and because natural apatites vary greatly in composition, the effect of other elements on U-Th/Pb geochronology hosted by these natural apatites needs to be investigated further using the same methodology. Because Pb diffusion occurs in apatites over micron length scales at temperatures relevant to the formation of the middle and lower crusts, thermal history information of the middle to lower crust and even at the crust-mantle boundary can be derived from U-Pb date profiles preserved within these apatites at Tc (500-1300 °C) depending on the chemical composition of the apatite. Also, apatites with high Tc provide more reliable ages than those with low Tc. (Krogstad and Walker, 1994). Further, based on the chemistry of the apatite, it can be determined whether it is useful for obtaining intermediate or high temperature petrochronology information.

5. Conclusions

A multi-scale approach was adopted in this study to demonstrate the effect of the chemical composition on the diffusion of Pb in apatite. The computed Pb insertion energies shows that it can be inserted with an energy of -2.03 eV, -2.21 eV, and -1.94 eV in pure F-apatite, Na₀.₁-apatite, and Cl₀.₂₅- Na₀.₁-apatite respectively. The diffusion of Pb in all these investigated apatites take places along the two preferential pathways i.e., one along the z-axis and the other along the plane orthogonal to z-axis which is termed as the in-plane direction. By using the KMC approach and the obtained quantum mechanically computed data, a 3D random walk simulation is performed to see the effect of variable chemical composition on the diffusion of Pb in apatite at the macro-scale/mineral scale level. The diffusion coefficients for the Pb atom at various temperatures are determined along each path and pre-exponential factor and Ea were extracted from MSDs using Arrhenius plot. Based on the Pb diffusion profiles estimated in different apatite lattices with different composition, one can conclude that Pb diffusion changes with chemical substitutions. These substitutions contribute towards the isotropic nature of Pb diffusion in apatite and conclude that the Pb-diffusion is strongly heterogenous depending on the apatite's chemistry. This work studies the Pb diffusion in terms of geochronological application and so the integration of all the findings leads us to a conclusion that the (U—Th)/Pb geochronology of apatite may be useful for predicting the age and geological processes of intermediate to high temperature thermal events in upper and lower crust of Earth and other terrestrial planets depending on the chemical composition of these apatites. Moreover, in order to improve the accuracy of apatite U-Th/Pb geochronology, a thorough understanding of apatite chemistry is required, i.e., the effect of other elements in the various natural groups of apatites must be investigated in the context of the Pb closure temperature. It is necessary because different elements in apatites have different effects on the apatite closure temperature, limiting its use in the investigation of crustal events below their respective Pb retentivity (or closure) temperature. This study also emphasizes that the findings regarding the diffusivity of partly Na- and/or Cl-substituted apatites infers about the importance of the adopted methodology in the characterization and computational modelling of special geological systems with spatial atomic arrangements, incorporating heterogeneities in crystals, and mobility of inserted atoms and their respective energetics profiles via integrated approach of KMC, DFT, NEB, and TST simulations.

Declaration of Competing Interest

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.

Acknowledgments

This research work is published from the PhD work of the first
author, funded by the Higher Education Commission of Pakistan under
the Faculty Development Program for the Khushal Khan Khattak Uni-
versity, Karak, Pakistan with grant No. 386-388/PC/PAB/ KKKUK/21.
Partial funds are provided by the University of Vienna, Austria. The
Department of Lithospheric Research, University of Vienna is also
acknowledged for providing computational facilities and access to the
online library. Thanks are extended to the editors and anonymous re-
viewers for providing feedback to improve the manuscript.

References

Andersen, M., Panosetti, C., Reuter, K., 2019. A practical guide to surface kinetic Monte Carlo simulations. *Front. Chem.* 7, 202.
Balout, H., Roques, J., Gautheron, C., Tassan-Got, L., Mbongo-Djimbi, D., 2017a. Helium diffusion in pure hematite (a-Fe2O3) for thermochronometric applications: a theoretical multi-scale study. Comput. Theor. Sci. 1099, 21–28.

Balout, H., Roques, J., Gautheron, C., Tassan-Got, L., 2017b. Computational investigation of the interstitial neon diffusion in pure hematite, a-Fe2O3. Comput. Mater. Sci. 128, 67–74.

Barnes, C.J., Majka, J., Jeanneret, P., Ziemiernik, G., Kooijman, E., Kośmińska, K., Schneider, D.A., 2021. Using Th-U-Pb geochronology to extract crystallization ages of Paleozoic metamorphic monazite contaminated by initial Pb. Chem. Geol. 582, 120450.

Bassal, F., Heller, B., Roques, J., Balout, H., Tassan-Got, L., Allard, T., Gautheron, C., 2022. Revealing the radiation damage and Al-content impacts on the diffusion in goethite. Chem. Geol. 611, 121118.

Bengtson, A., Ewing, R.C., Becker, U., 2012. He diffusion and closure temperatures in apatite and zircon: a density functional theory investigation. Geochim. Cosmochim. Acta 86, 228–238.

Blackburn, T., Bowring, S., Schoene, B., Mahan, K., Dudas, F., 2011. U-Pb thermochronology: creating a temporal record of lithosphere thermal evolution. Contrib. Mineral. Petrol. 162, 479–500.

Blackburn, T.J., Bowring, S.A., Perron, J.T., Mahan, K.H., Dudas, F.O., Barnhart, K.R., 2012. An exhumation history of continents over billion-year time scales. Science 335, 73–76. https://doi.org/10.1126/science.1213496.

Bortz, A.B., Kalos, M.H., Lebowitz, J.L., 1975. A new algorithm for Monte Carlo simulation of Ising spin systems. J. Comput. Phys. 17 (1), 10–18.

Chakhmouradian, A.R., Medici, L., 2006. Clinohydroxylapatite: a new apatite-group mineral from northwestern Ontario (Canada), and new data on the extent of Na-S substitution in natural apatites. Eur. J. Mineral. 18 (1), 105–112.

Chamberlain, K.R., Bowring, S.A., 2001. Apatite-feldspar U-Pb thermochronometer: a reliable, mid-range (~450°C), diffusion-controlled system. Chem. Geol. 172, 173–200.

Cherniak, D.J., 2005. Uranium and manganese diffusion in apatite. Chem. Geol. 219 (1-4), 297-308.

Cherniak, D.J., 2010. Diffusion in accessory minerals: zircon, titanite, apatite, monazite and xenotime. Rev. Mineral. Geochim. Acta 72 (1), 827–869.

Cherniak, D.J., Watson, E.B., 2001. Pb diffusion in zircon. Chem. Geol. 172 (1-2), 5-24.

Cherniak, D.J., Lanford, W.A., Ryerson, F.J., 1991. Lead diffusion in apatite and zircon using ion implantation and Rutherford backscattering techniques. Geochim. Cosmochim. Acta 55 (6), 1663–1673.

Cherniak, D.J., Watson, E.B., Grove, T.M., Harrison, T.M., 2004. Pb diffusion in monazite: a combined RBS/SIMS study. Geochim. Cosmochim. Acta 68 (4), 829–840.

Chew, D.M., Spikings, R.A., 2015. Geochronology and thermochronology using apatite: time and temperature, lower crust to surface. Elements 11 (3), 189–194.

Chew, D.M., Sylvester, P.J., Tubrett, M.N., 2011. U-Pb and Th-Pb dating of apatite by LA-ICPMS. Chem. Geol. 280, 200–216.

Chew, D.M., Babechuk, M.G., Cogne, N., Mark, C., O'Sullivan, G.J., Henrichs, I.A., Doepke, D., McKenna, C.A., 2016. (LA, Q)-ICPMS trace-element analyses of Durango and McClure Mountain apatite and implications for making natural LA-ICPMS mineral standards. Chem. Geol. 435, 35–48.

Cocrane, R., Spikings, R.A., Chew, D., Wotzlaw, J.F., Chiaradia, M., Tyrrell, S., Van der Lelij, R., 2014. High temperature (> 350 C) thermochronology and mechanisms of Pb loss in apatite. Geochim. Cosmochim. Acta 127, 39–56.

Danišk, M., McInnes, B.I., Kirkland, C.L., McDonald, B.J., Evans, N.J., Becker, T., 2017. Seeing is believing: Visualization of He distribution in zircon and implications for thermal history reconstruction on single crystals. Sci. Adv. 3 (2), e1601121.

**Digis Team, 2022. GEOROC Compilation: Melt Inclusions. GRO.data, Göttingen, Germany. https://doi.org/10.25625/7JW6XU.**

Djimbi, D.M., Gautheron, C., Roques, J., Tassan-Got, L., Gerin, C., Simoni, E., 2015. Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view. Geochim. Cosmochim. Acta 167, 162–176.

Dodson, M.H., 1973. Closure temperature in cooling geochronological and petrological systems. Contrib. Mineral. Petrol. 40 (3), 259–274.

Farley, K.A., 2000. Helium diffusion from apatite: General behavior as illustrated by Durango fluorapatite. J. Geophys. Res. Solid Earth 105 (B2), 2903–2914.

Farley, K.A., Stockli, D.F., 2002. (U-Th)/He dating of phosphates: apatite, monazite, and xenotime. Rev. Mineral. Geochim. Acta 48 (1), 559–577.

Finger, L.W., Hazen, R.M., 1980. Crystal structure and isothermal compression of Fe2O3, Cr2O3, and V2O3 to 50 kbars. J. Appl. Phys. 51 (10), 5362–5367.

Flowers, R.M., Shuster, D.L., Wernicke, B.P., Farley, K.A., 2007. Radiation damage control on apatite (U-Th)/He dates from the Grand Canyon region, Colorado Plateau. Geology 35 (5), 447–450.

Ganguly, J., Tirone, M., 1999. Diffusion closure temperature and age of a mineral with arbitrary extent of diffusion: theoretical formulation and applications. Earth Planet. Sci. Lett. 170 (1-2), 131–140.

Gautheron, C., Djimbi, D.M., Roques, J., Balout, H., Ketcham, R.A., Simoni, E., Tassan-Got, L., 2020. A multi-method, multi-scale theoretical study of He and Ne diffusion in zircon. Geochim. Cosmochim. Acta 268, 348–367.

Gerin, C., Gautheron, C., Oliviero, E., Bachelet, C., Djimbi, M.D., Seydoux-Guillaume, A. M., Tassan-Got, L., Sarda, P., Roques, J., Garrido, F., 2017. Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales. Geochim. Cosmochim. Acta 197, 87–103.

Gianzotti, P., Baroni, S., Bonini, N., Calandra, M., Car, R., Cavazzoni, C., Wentzcovitch, R. M., 2009. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys. Condens. Matter 21 (39), 395502.

Glorie, S., Alexandrov, I., Nixon, A., Jepson, G., Gillespie, J., Jahn, B.M., 2017. Thermal and exhumation history of Sakhalin Island (Russia) constrained by apatite U-Pb and fission track thermochronology. J. Asian Earth Sci. 143, 326–342.

Guzmics, T., Mitchell, R.H., Szabó, C., Berkesi, M., Milke, R., Abart, R., 2011. Carbonatite melt inclusions in coexisting magnetite, apatite and monticellite in Kerimasi calciocarbonate, Tanzania: melt evolution and petrogenesis. Contrib. Mineral. Petrol. 161, 177–196.

Guzmics, T., Zajacz, Z., Mitchell, R.H., Szabó, C., Wälle, M., 2015. The role of liquid-liquid immiscibility and crystal fractionation in the genesis of carbonatite magmas: insights from Kerimasi melt inclusions. Contrib. Mineral. Petrol. 169, 1–18.

Henrichs, I.A., O'Sullivan, G., Chew, D.M., Mark, C., Babechuk, M.G., McKenna, C., Emo, R., 2018. The trace element and U-Pb systematics of metamorphic apatite. Chem. Geol. 483, 218–238.

Hoffmann, M.J., Matera, S., Reuter, K., 2014. kmos: a lattice kinetic Monte Carlo framework. Comput. Phys. Commun. 185 (7), 2138–2150.

Hohenberg, P., Kohn, W., 1964. Inhomogeneous electron gas. Phys. Rev. 136 (3B), B864-B871.

Hughes, J.M., Cameron, M., Crowley, K.D., 1989. Structural variations in natural F, OH, and Cl apatites. Am. Mineral. 74 (7-8), 870–876.

Jo'nnsson, H., Mills, G., Jacobsen, K.W., Coker, D.F., 1998. Nudged elastic band method for finding minimum energy paths of transitions. In: Berne, B.J., Ciccotti, G. (Eds.), Classical and Quantum Dynamics in Condensed Phase Simulations. World Scientific, pp. 385–404.

Kirkland, C.L., Hollis, J., Danišk, M., Petersen, J., Evans, N.J., McDonald, B.J., 2017. Apatite and titanite from the Karrat Group, Greenland; implications for charting the thermal evolution of crust from the U-Pb geochronology of common Pb bearing phases. Precambrian Res. 300, 107–120.

Kirkland, C.L., Yakymchuk, C., Szilas, K., Evans, N., Hollis, J., McDonald, B., Gardiner, N. J., 2018. Apatite: a U-Pb thermochronometer or geochronometer? Lithos 318, 143–157.

Kohn, W., Sham, L.J., 1965. Self-consistent equations including exchange and correlation effects. Phys. Rev. 140, A1133–A1138.

Krogstad, E.J., Walker, R.J., 1994. High closure temperatures of the U-Pb system in large apatites from the Tin Mountain pegmatite, Black Hills, South Dakota, USA. Geochim. Cosmochim. Acta 58 (18), 3845–3853.

LaFlamme, C., McFarlane, C.R., Fisher, C.M., Kirkland, C.L., 2017. Multi-mineral geochronology: insights into crustal behaviour during exhumation of an orogenic root. Contrib. Mineral. Petrol. 172, 1–18.

Liferovich, R.P., Mitchell, R.H., 2006. Apatite-group minerals from napheline syenite, Pilansberg alkaline complex, South Africa. Mineral. Mag. 70 (5), 463–484.

Lu, J., Chen, W., Ying, Y., Jiang, S., Zhao, K., 2021. Apatite texture and trace element chemistry of carbonate-related REE deposits in China: implications for petrogenesis. Lithos 398, 106276.

Mark, C., Cogné, N., Chew, D., 2016. Tracking exhumation and drainage divide migration of the Western Alps: a test of the apatite U-Pb thermochronometer as a detrital provenance tool. Bulletin 128 (9-10), 1439–1460.

McGregor, M., Dence, M.R., McFarlane, C.R., Spray, J.G., 2020. U-Pb geochronology of apatite and zircon from the Brent impact structure, Canada: a Late Ordovician Sandbian-Katian boundary event associated with L-Chondrite parent body disruption. Contrib. Mineral. Petrol. 175, 1–21.

Mills, G., Jonsson, H., Schenter, G.K., 1995. Reversible work transition-state theory:

application to dissociative adsorption of hydrogen. Surf. Sci. 324, 305–337.

Miro, S., Studer, F., Costantini, J.M., Haussy, J., Troulard, P., Grob, J.J., 2006. Effect of composition on helium diffusion in fluoroapatites investigated with nuclear reaction analysis. J. Nucl. Mater. 355 (1-3), 1–9.

Pozun, Z.D., Henkelman, G., 2011. Hybrid density functional theory band structure engineering in hematite. J. Chem. Phys. 134 (22), 224706.

Pratt, L.R., 1986. A statistical method for identifying transition states in high dimensional problems. J. Chem. Phys. 85, 5045–5048.

Ptáček, P., 2016. Substituents and dopants in the structure of apatite. In: Apatites and their Synthetic Analogues—Synthesis, Structure, Properties and Applications. Open Acces Monograph. InTech, pp. 289–334.

Reich, M., Ewing, R.C., Ehlers, T., Becker, U., 2007. Lowtemperature anisotropic diffusion of helium in zircon: implications for zircon (U-Th)/He thermochronometry. Geochim. Cosmochim. Acta 71, 3119–3130.

Reiners, P.W., Brandon, M.T., 2006. Using thermochronology to understand orogenic erosion. Annu. Rev. Earth Planet. Sci. 34, 419–466.

Reiners, P.W., Spell, T.L., Nicolescu, S., Zanetti, K.A., 2004. Zircon (U-Th)/He thermochronometry: he diffusion and comparison with 40Ar/39Ar dating. Geochim. Cosmochim. Acta 68, 1857–1887.

Rønsbo, J.G., 1989. Coupled substitutions involving REEs and Na and Si in apatites in alkaline rocks from the Ilimaussaq intrusion, South Greenland, and the petrological implications. Am. Mineral. 74 (7-8), 896–901.

Saadoune, I., Purton, J.A., de Leeuw, N.H., 2009. He incorporation and diffusion pathways in pure and defective zircon ZrSiO4: a density functional theory study. Chem. Geol. 258, 182–196.

Schoene, B., 2014. In: Turekian, Karl K., Holland, H. (Eds.), U-Th-Pb Geochronology. Treatise on Geochemistry, Vol. 4. The Crust. 16 vols.

Schoene, B., Bowring, S.A., 2007. Determining accurate temperature-time paths from U-Pb thermochronology: an example from the Kaapvaal craton, southern Africa. Geochim. Cosmochim. Acta 71, 165–185.

Smye, A.J., Marsh, J.H., Vermeesch, P., Garber, J.M., Stockli, D.F., 2018. Applications and limitations of U-Pb thermochronology to middle and lower crustal thermal histories. Chem. Geol. 494, 1–18.

Spikings, R.A., Crowhurst, P.V., Winkler, W., Villagomez, D., 2010. Syn- and post-accretionary cooling history of the Ecuadorian Andes constrained by their in-situ and detrital thermochronometric record. J. S. Am. Earth Sci. 30, 121–133.

Tamm, T., Peld, M., 2006. Computational study of cation substitutions in apatites. J. Solid State Chem. 179 (5), 1581–1587.
Balout, H., Roques, J., Gautheron, C., Tassan-Got, L., Mbongo-Djimbi, D., 2017a. Helium diffusion in pure hematite (a-Fe2O3) for thermochronometric applications: a theoretical multi-scale study. Comput. Theor. Sci. 1099, 21–28.

Balout, H., Roques, J., Gautheron, C., Tassan-Got, L., 2017b. Computational investigation of the interstitial neon diffusion in pure hematite, a-Fe2O3. Comput. Mater. Sci. 128, 67–74.

Barnes, C.J., Majka, J., Jeanneret, P., Ziemiernik, G., Kooijman, E., Kośmińska, K., Schneider, D.A., 2021. Using Th-U-Pb geochronology to extract crystallization ages of Paleozoic metamorphic monazite contaminated by initial Pb. Chem. Geol. 582, 120450.

Bassal, F., Heller, B., Roques, J., Balout, H., Tassan-Got, L., Allard, T., Gautheron, C., 2022. Revealing the radiation damage and Al-content impacts on the diffusion in goethite. Chem. Geol. 611, 121118.

Bengtson, A., Ewing, R.C., Becker, U., 2012. He diffusion and closure temperatures in apatite and zircon: a density functional theory investigation. Geochim. Cosmochim. Acta 86, 228–238.

Blackburn, T., Bowring, S., Schoene, B., Mahan, K., Dudas, F., 2011. U-Pb thermochronology: creating a temporal record of lithosphere thermal evolution. Contrib. Mineral. Petrol. 162, 479–500.

Blackburn, T.J., Bowring, S.A., Perron, J.T., Mahan, K.H., Dudas, F.O., Barnhart, K.R., 2012. An exhumation history of continents over billion-year time scales. Science 335, 73–76. https://doi.org/10.1126/science.1213496.

Bortz, A.B., Kalos, M.H., Lebowitz, J.L., 1975. A new algorithm for Monte Carlo simulation of Ising spin systems. J. Comput. Phys. 17 (1), 10–18.

Chakhmouradian, A.R., Medici, L., 2006. Clinohydroxylapatite: a new apatite-group mineral from northwestern Ontario (Canada), and new data on the extent of Na-S substitution in natural apatites. Eur. J. Mineral. 18 (1), 105–112.

Chamberlain, K.R., Bowring, S.A., 2001. Apatite-feldspar U-Pb thermochronometer: a reliable, mid-range (~450°C), diffusion-controlled system. Chem. Geol. 172, 173–200.

Cherniak, D.J., 2005. Uranium and manganese diffusion in apatite. Chem. Geol. 219 (1-4), 297-308.

Cherniak, D.J., 2010. Diffusion in accessory minerals: zircon, titanite, apatite, monazite and xenotime. Rev. Mineral. Geochim. Acta 72 (1), 827–869.

Cherniak, D.J., Watson, E.B., 2001. Pb diffusion in zircon. Chem. Geol. 172 (1-2), 5–24.

Cherniak, D.J., Lanford, W.A., Ryerson, F.J., 1991. Lead diffusion in apatite and zircon using ion implantation and Rutherford backscattering techniques. Geochim. Cosmochim. Acta 55 (6), 1663–1673.

Cherniak, D.J., Watson, E.B., Grove, T.M., Harrison, T.M., 2004. Pb diffusion in monazite: a combined RBS/SIMS study. Geochim. Cosmochim. Acta 68 (4), 829–840.

Chew, D.M., Spikings, R.A., 2015. Geochronology and thermochronology using apatite: time and temperature, lower crust to surface. Elements 11 (3), 189–194.

Chew, D.M., Sylvester, P.J., Tubrett, M.N., 2011. U-Pb and Th-Pb dating of apatite by LA-ICPMS. Chem. Geol. 280, 200–216.

Chew, D.M., Babechuk, M.G., Cogne, N., Mark, C., O'Sullivan, G.J., Henrichs, I.A., Doepke, D., McKenna, C.A., 2016. (LA, Q)-ICPMS trace-element analyses of Durango and McClure Mountain apatite and implications for making natural LA-ICPMS mineral standards. Chem. Geol. 435, 35–48.

Cocrane, R., Spikings, R.A., Chew, D., Wotzlaw, J.F., Chiaradia, M., Tyrrell, S., Van der Lelij, R., 2014. High temperature (> 350 C) thermochronology and mechanisms of Pb loss in apatite. Geochim. Cosmochim. Acta 127, 39–56.

Danišk, M., McInnes, B.I., Kirkland, C.L., McDonald, B.J., Evans, N.J., Becker, T., 2017. Seeing is believing: Visualization of He distribution in zircon and implications for thermal history reconstruction on single crystals. Sci. Adv. 3 (2), e1601121.

**Digis Team, 2022. GEOROC Compilation: Melt Inclusions. GRO.data, Göttingen, Germany. https://doi.org/10.25625/7JW6XU.**

Djimbi, D.M., Gautheron, C., Roques, J., Tassan-Got, L., Gerin, C., Simoni, E., 2015. Impact of apatite chemical composition on (U-Th)/He thermochronometry: an atomistic point of view. Geochim. Cosmochim. Acta 167, 162–176.

Dodson, M.H., 1973. Closure temperature in cooling geochronological and petrological systems. Contrib. Mineral. Petrol. 40 (3), 259–274.

Farley, K.A., 2000. Helium diffusion from apatite: General behavior as illustrated by Durango fluorapatite. J. Geophys. Res. Solid Earth 105 (B2), 2903–2914.

Farley, K.A., Stockli, D.F., 2002. (U-Th)/He dating of phosphates: apatite, monazite, and xenotime. Rev. Mineral. Geochim. Acta 48 (1), 559–577.

Finger, L.W., Hazen, R.M., 1980. Crystal structure and isothermal compression of Fe2O3, Cr2O3, and V2O3 to 50 kbars. J. Appl. Phys. 51 (10), 5362–5367.

Flowers, R.M., Shuster, D.L., Wernicke, B.P., Farley, K.A., 2007. Radiation damage control on apatite (U-Th)/He dates from the Grand Canyon region, Colorado Plateau. Geology 35 (5), 447–450.

Ganguly, J., Tirone, M., 1999. Diffusion closure temperature and age of a mineral with arbitrary extent of diffusion: theoretical formulation and applications. Earth Planet. Sci. Lett. 170 (1-2), 131–140.

Gautheron, C., Djimbi, D.M., Roques, J., Balout, H., Ketcham, R.A., Simoni, E., Tassan-Got, L., 2020. A multi-method, multi-scale theoretical study of He and Ne diffusion in zircon. Geochim. Cosmochim. Acta 268, 348–367.

Gerin, C., Gautheron, C., Oliviero, E., Bachelet, C., Djimbi, M.D., Seydoux-Guillaume, A. M., Tassan-Got, L., Sarda, P., Roques, J., Garrido, F., 2017. Influence of vacancy damage on He diffusion in apatite investigated at atomic to mineralogical scales. Geochim. Cosmochim. Acta 197, 87–103.

Gianzotti, P., Baroni, S., Bonini, N., Calandra, M., Car, R., Cavazzoni, C., Wentzcovitch, R. M., 2009. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys. Condens. Matter 21 (39), 395502.

Glorie, S., Alexandrov, I., Nixon, A., Jepson, G., Gillespie, J., Jahn, B.M., 2017. Thermal and exhumation history of Sakhalin Island (Russia) constrained by apatite U-Pb and fission track thermochronology. J. Asian Earth Sci. 143, 326–342.

Guzmics, T., Mitchell, R.H., Szabó, C., Berkesi, M., Milke, R., Abart, R., 2011. Carbonatite melt inclusions in coexisting magnetite, apatite and monticellite in Kerimasi calciocarbonate, Tanzania: melt evolution and petrogenesis. Contrib. Mineral. Petrol. 161, 177–196.

Guzmics, T., Zajacz, Z., Mitchell, R.H., Szabó, C., Wälle, M., 2015. The role of liquid-liquid immiscibility and crystal fractionation in the genesis of carbonatite magmas: insights from Kerimasi melt inclusions. Contrib. Mineral. Petrol. 169, 1–18.

Henrichs, I.A., O'Sullivan, G., Chew, D.M., Mark, C., Babechuk, M.G., McKenna, C., Emo, R., 2018. The trace element and U-Pb systematics of metamorphic apatite. Chem. Geol. 483, 218–238.

Hoffmann, M.J., Matera, S., Reuter, K., 2014. kmos: a lattice kinetic Monte Carlo framework. Comput. Phys. Commun. 185 (7), 2138–2150.

Hohenberg, P., Kohn, W., 1964. Inhomogeneous electron gas. Phys. Rev. 136 (3B), B864-B871.

Hughes, J.M., Cameron, M., Crowley, K.D., 1989. Structural variations in natural F, OH, and Cl apatites. Am. Mineral. 74 (7-8), 870–876.

Jo'nnsson, H., Mills, G., Jacobsen, K.W., Coker, D.F., 1998. Nudged elastic band method for finding minimum energy paths of transitions. In: Berne, B.J., Ciccotti, G. (Eds.), Classical and Quantum Dynamics in Condensed Phase Simulations. World Scientific, pp. 385–404.

Kirkland, C.L., Hollis, J., Danišk, M., Petersen, J., Evans, N.J., McDonald, B.J., 2017. Apatite and titanite from the Karrat Group, Greenland; implications for charting the thermal evolution of crust from the U-Pb geochronology of common Pb bearing phases. Precambrian Res. 300, 107–120.

Kirkland, C.L., Yakymchuk, C., Szilas, K., Evans, N., Hollis, J., McDonald, B., Gardiner, N. J., 2018. Apatite: a U-Pb thermochronometer or geochronometer? Lithos 318, 143–157.

Kohn, W., Sham, L.J., 1965. Self-consistent equations including exchange and correlation effects. Phys. Rev. 140, A1133–A1138.

Krogstad, E.J., Walker, R.J., 1994. High closure temperatures of the U-Pb system in large apatites from the Tin Mountain pegmatite, Black Hills, South Dakota, USA. Geochim. Cosmochim. Acta 58 (18), 3845–3853.

LaFlamme, C., McFarlane, C.R., Fisher, C.M., Kirkland, C.L., 2017. Multi-mineral geochronology: insights into crustal behaviour during exhumation of an orogenic root. Contrib. Mineral. Petrol. 172, 1–18.

Liferovich, R.P., Mitchell, R.H., 2006. Apatite-group minerals from napheline syenite, Pilansberg alkaline complex, South Africa. Mineral. Mag. 70 (5), 463–484.

Lu, J., Chen, W., Ying, Y., Jiang, S., Zhao, K., 2021. Apatite texture and trace element chemistry of carbonate-related REE deposits in China: implications for petrogenesis. Lithos 398, 106276.

Mark, C., Cogné, N., Chew, D., 2016. Tracking exhumation and drainage divide migration of the Western Alps: a test of the apatite U-Pb thermochronometer as a detrital provenance tool. Bulletin 128 (9-10), 1439–1460.

McGregor, M., Dence, M.R., McFarlane, C.R., Spray, J.G., 2020. U-Pb geochronology of apatite and zircon from the Brent impact structure, Canada: a Late Ordovician Sandbian-Katian boundary event associated with L-Chondrite parent body disruption. Contrib. Mineral. Petrol. 175, 1–21.

Mills, G., Jonsson, H., Schenter, G.K., 1995. Reversible work transition-state theory:

application to dissociative adsorption of hydrogen. Surf. Sci. 324, 305–337.

Miro, S., Studer, F., Costantini, J.M., Haussy, J., Troulard, P., Grob, J.J., 2006. Effect of composition on helium diffusion in fluoroapatites investigated with nuclear reaction analysis. J. Nucl. Mater. 355 (1-3), 1–9.

Pozun, Z.D., Henkelman, G., 2011. Hybrid density functional theory band structure engineering in hematite. J. Chem. Phys. 134 (22), 224706.

Pratt, L.R., 1986. A statistical method for identifying transition states in high dimensional problems. J. Chem. Phys. 85, 5045–5048.

Pátek, P., 2016. Substituents and dopants in the structure of apatite. In: Apatites and their Synthetic Analogues—Synthesis, Structure, Properties and Applications. Open Acces Monograph. InTech, pp. 289–334.

Reich, M., Ewing, R.C., Ehlers, T., Becker, U., 2007. Lowtemperature anisotropic diffusion of helium in zircon: implications for zircon (U-Th)/He thermochronometry. Geochim. Cosmochim. Acta 71, 3119–3130.

Reiners, P.W., Brandon, M.T., 2006. Using thermochronology to understand orogenic erosion. Annu. Rev. Earth Planet. Sci. 34, 419–466.

Reiners, P.W., Spell, T.L., Nicolescu, S., Zanetti, K.A., 2004. Zircon (U-Th)/He thermochronometry: he diffusion and comparison with 40Ar/39Ar dating. Geochim. Cosmochim. Acta 68, 1857–1887.

Rønsbo, J.G., 1989. Coupled substitutions involving REEs and Na and Si in apatites in alkaline rocks from the Ilimaussaq intrusion, South Greenland, and the petrological implications. Am. Mineral. 74 (7-8), 896–901.

Saadoune, I., Purton, J.A., de Leeuw, N.H., 2009. He incorporation and diffusion pathways in pure and defective zircon ZrSiO4: a density functional theory study. Chem. Geol. 258, 182–196.

Schoene, B., 2014. In: Turekian, Karl K., Holland, H. (Eds.), U-Th-Pb Geochronology. Treatise on Geochemistry, Vol. 4. The Crust. 16 vols.

Schoene, B., Bowring, S.A., 2007. Determining accurate temperature-time paths from U-Pb thermochronology: an example from the Kaapvaal craton, southern Africa. Geochim. Cosmochim. Acta 71, 165–185.

Smye, A.J., Marsh, J.H., Vermeesch, P., Garber, J.M., Stockli, D.F., 2018. Applications and limitations of U-Pb thermochronology to middle and lower crustal thermal histories. Chem. Geol. 494, 1–18.

Spikings, R.A., Crowhurst, P.V., Winkler, W., Villagomez, D., 2010. Syn- and post-accretionary cooling history of the Ecuadorian Andes constrained by their in-situ and detrital thermochronometric record. J. S. Am. Earth Sci. 30, 121–133.

Tamm, T., Peld, M., 2006. Computational study of cation substitutions in apatites. J. Solid State Chem. 179 (5), 1581–1587.
Tang, M., Rudnick, R.L., McDonough, W.F., Bose, M., Goreva, Y., 2017. Multi-mode Li diffusion in natural zircons: evidence for diffusion in the presence of step-function concentration boundaries. Earth Planet. Sci. Lett. 474, 110–119.

Ullah, M., Klötzli, U., Wadood, B., Khubab, M., Islam, F.U., Shehzad, K., Ahmad, R., 2023. Pb diffusion in perfect and defective zircon for thermo/petro-chronological investigations: an atomistic approach. Chem. Geol. 121750.

Vineyard, G.H., 1957. Frequency factors and isotope effects in solid state rate processes. J. Phys. Chem. Solids 3, 121–127.

Voter, A.F., 1986. Classically exact overlayer dynamics: diffusion of rhodium clusters on Rh (100). Phys. Rev. B 34 (10), 6819.

Watson, E.B., Baxter, E.F., 2007. Diffusion in solid-Earth systems. Earth and Planetary Science Letters 253 (3–4), 307–327.

Watson, E.B., Harrison, T.M., Ryerson, F.J., 1985. Diffusion of Sm, Sr, and Pb in fluorapatite. Geochim. Cosmochim. Acta 49 (8), 1813–1823.

Webster, J.D., Piccoli, P.M., 2015. Magmatic apatite: a powerful, yet deceptive, mineral. Elements 11, 177–182.

Wert, C., Zener, C., 1949. Interstitial atomic diffusion coefficients. Phys. Rev. 76, 1169–1175.

Wohlgemuth-Ueberwasser, C.C., Tegner, C., Pease, V., 2017. LA-Q-ICP-MS apatite U/Pb geochronology using common Pb in plagioclase: examples from layered mafic intrusions. Am. Mineral. 102 (3), 571–579.

Yang, Y.H., Wu, F.Y., Yang, J.H., Chew, D.M., Xie, L.W., Chu, Z.Y., Huang, C., 2014. Sr and Nd isotopic compositions of apatite reference materials used in U–Th–Pb geochronology. Chem. Geol. 385, 35–55.

Zeitler, P.K., Herczeg, A.L., McDougall, I., Honda, M., 1987. U-Th-He dating of apatite: a potential thermochronometer. Geochim. Cosmochim. Acta 51 (10), 2865–2868.

Zeitler, P.K., Enkelmann, E., Thomas, J.B., Watson, E.B., Ancuta, L.D., Idleman, B.D., 2017. Solubility and trapping of helium in apatite. Geochim. Cosmochim. Acta 209, 1–8.

Zhang, J., Huang, W., Wu, J., Liang, H., Lin, S., 2021. Geological implications of apatite within the granite-related Jiepai W-(Cu) deposit, Guangxi, South China. Ore Geol. Rev. 139, 104548.
Tang, M., Rudnick, R.L., McDonough, W.F., Bose, M., Goreva, Y., 2017. Multi-mode Li diffusion in natural zircons: evidence for diffusion in the presence of step-function concentration boundaries. Earth Planet. Sci. Lett. 474, 110–119.

Ullah, M., Klötzli, U., Wadood, B., Khubab, M., Islam, F.U., Shehzad, K., Ahmad, R., 2023. Pb diffusion in perfect and defective zircon for thermo/petro-chronological investigations: an atomistic approach. Chem. Geol. 121750.

Vineyard, G.H., 1957. Frequency factors and isotope effects in solid state rate processes. J. Phys. Chem. Solids 3, 121–127.

Voter, A.F., 1986. Classically exact overlayer dynamics: diffusion of rhodium clusters on Rh (100). Phys. Rev. B 34 (10), 6819.

Watson, E.B., Baxter, E.F., 2007. Diffusion in solid-Earth systems. Earth and Planetary Science Letters 253 (3–4), 307–327.

Watson, E.B., Harrison, T.M., Ryerson, F.J., 1985. Diffusion of Sm, Sr, and Pb in fluorapatite. Geochim. Cosmochim. Acta 49 (8), 1813–1823.

Webster, J.D., Piccoli, P.M., 2015. Magmatic apatite: a powerful, yet deceptive, mineral. Elements 11, 177–182.

Wert, C., Zener, C., 1949. Interstitial atomic diffusion coefficients. Phys. Rev. 76, 1169–1175.

Wohlgemuth-Ueberwasser, C.C., Tegner, C., Pease, V., 2017. LA-Q-ICP-MS apatite U/Pb geochronology using common Pb in plagioclase: examples from layered mafic intrusions. Am. Mineral. 102 (3), 571–579.

Yang, Y.H., Wu, F.Y., Yang, J.H., Chew, D.M., Xie, L.W., Chu, Z.Y., Huang, C., 2014. Sr and Nd isotopic compositions of apatite reference materials used in U–Th–Pb geochronology. Chem. Geol. 385, 35–55.

Zeitler, P.K., Herczeg, A.L., McDougall, I., Honda, M., 1987. U-Th-He dating of apatite: a potential thermochronometer. Geochim. Cosmochim. Acta 51 (10), 2865–2868.

Zeitler, P.K., Enkelmann, E., Thomas, J.B., Watson, E.B., Ancuta, L.D., Idleman, B.D., 2017. Solubility and trapping of helium in apatite. Geochim. Cosmochim. Acta 209, 1–8.

Zhang, J., Huang, W., Wu, J., Liang, H., Lin, S., 2021. Geological implications of apatite within the granite-related Jiepai W-(Cu) deposit, Guangxi, South China. Ore Geol. Rev. 139, 104548.