Molecular Dynamics Approach for Predicting Release Temperatures of Noble Gases in Presolar Nanodiamonds

Alireza Aghajamali; Andrey A. Shiryaev; Nigel A. Marks

doi:10.3847/1538-4357/ac06cf

1. Introduction

Concentration of noble gases (NGs) in primitive chondrites may greatly exceed terrestrial abundances. Several NG components, different in chemical and isotopic composition, can be distinguished using step combustion experiments (see Ott 2014 for a review). A significant fraction of the anomalous components, i.e., those with an isotopic composition markedly different from the solar, reside in nanodiamond (ND) grains, which are relatively abundant in some primitive chondrites with content as high as 1500 ppm.

Step combustion and pyrolysis experiments allow one to distinguish three NG components present in NDs: P3, HL, and P6. Each of these components contains all five NGs (He, Ne, Ar, Kr, and Xe); understanding of Xe is, arguably, less ambiguous than the other elements. The P3 and, possibly, P6 components are not dramatically different from the solar elemental and isotopic patterns, but still, their composition cannot be derived from the solar by fractionation.

The P6 component is volumetrically the least important and yet is insufficiently characterized. At least a fraction of the P3- and P6-related ¹²⁹Xe is derived from the decay of ¹²⁹I, implying the initial presence of currently extinct iodine in NDs (e.g., Gilmour et al. 2005). The HL component is characterized by a marked excess in the relative abundance of heavy (H) and light (L) isotopes. The heavy isotopes—¹³⁴Xe and ¹³⁶Xe—are produced by r-process only, or, at least, overwhelmingly. The light isotopes—¹²⁴Xe and ¹²⁶Xe—are formed via the p-process. Historically, the only feasible explanation for the origin of the HL component was a mixture of NGs produced in p- and r-processes in supernova explosions of type I (Jørgensen 1988) or II (Heymann & Dziczkaniec 1979). However, none of the current models of the r-process fit the experimentally observed isotopic pattern, and models with variable degrees of robustness were proposed (Ott et al. 2012; Ott 2014). In any case, isotopic patterns suggest that NDs carrying the HL and likely also the P3 and P6 components should be of presolar origin. This view is also supported by the association of isotopically anomalous Pt, Ba, Te, and some other elements with ND-rich fractions. Importantly, all attempts to separate individual carriers of the L and H components failed, thus indicating thorough mixing of several independent populations of NDs in presolar nebulae or elsewhere. Since only approximately one out of 10⁶ ND grains contains a Xe ion, and all isotopic analyses require billions of grains, the mixing is not that surprising. The existence of several populations of NDs, possibly of different origins and/or stellar sources, is also suggested by C and N isotopic composition (Russell et al. 1996).

Step pyrolysis/oxidation experiments revealed complex patterns of NG release from NDs. For heavy NGs—Ar, Kr, and Xe—the gas release pattern is bimodal. Although some variations between meteorites with different thermal histories occur, in general, the P3 component is released mostly at low temperatures (200–900) with a peak at 450°C–550°C, the HL component peaks between 1300°C and 1450°C (overall range 1100°C–1600°C), and the P6 component evolves at yet higher temperatures. For Ne, the bimodality is much less pronounced, and for He, it is probably absent, although its composition changes with the release temperature.

The mechanism of incorporation of NGs into ND grains is a subject of ongoing debate. In principle, three mechanisms may be responsible for the introduction of an impurity atom into an ND grain: (1) the growth process, (2) diffusion and/or adsorption, and (3) ion implantation. The largest atom ever introduced into the diamond lattice purely by growth processes relevant for the formation process of NDs in space—plasma-assisted chemical vapor deposition—is tin (Westerhausen et al. 2020). However, in contrast to NGs, Sn, Ge, and other heavy elements belong to group IV of the periodic table, as well as C. Subsequently, the introduction of these elements is crystallochemically feasible even despite large differences in atomic radii with C. Diffusion in diamond is extremely sluggish (Koga et al. 2003), which effectively rules out any diffusion-related processes at low temperatures, relevant for meteoritic NDs. The surface adsorption of some elements is possible, and several works suggest that the P3 component might be trapped in the near-surface, perhaps partly graphitized layers of the grains (e.g., Fisenko & Semjonova 2010). However, ion implantation, where NGs are driven into preexisting NDs (Lewis & Anders 1981; Lewis et al. 1987; Verchovsky et al. 1998; Gilmour et al. 2005), seems to be the only feasible option (Huss & Lewis 1994a).

The implantation hypothesis was strongly supported by experimental work by Koscheev et al. (2001), who implanted low-energy (700 eV) NG ions into synthetic NDs with diameters of 4–5 nm. Pyrolysis measurements observed a single broad temperature distribution for He and Ne and a bimodal distribution for Ar, Kr, and Xe, corresponding closely to those found for meteoritic NDs (Huss & Lewis 1994a, 1994b). Similar studies by Verchovsky et al. (2000) also found a bimodal distribution, albeit in the latter case, the high-dose implantation most likely amorphized the ND grains, greatly complicating interpretation. However, the physical mechanism of the bimodal gas release pattern remains unexplained and raises important questions about whether the observed multimodality is real or results from experimental artifacts. For example, a single layer of ND grains was used as a target in Koscheev et al. (2001). However, formation of a uniform truly single layer is notoriously difficult, and it is impossible to exclude the possibility that an ion penetrates the ND particle in a top layer and is implanted with smaller energy into the underlying grain. In addition, the presence of doubly charged ions during the implantation experiment is not fully excluded.

In this work, we use the molecular dynamics (MD) approach and apply it to the question of the single and bimodal temperature release of NGs. To the best of our knowledge, atomistic simulation methods have not previously been applied to this problem. We simulate the implantation and pyrolysis processes for He and Xe and find excellent agreement with experimental data collected from NDs. We reproduce the single release peak for He and bimodal release for Xe. The simulations explain how the mass, implantation depth, and crystallographic location of the NGs give rise to the experimental observations and for the first time provide an atomistic explanation of the release mechanism. Although we do not provide an explanation for the origin of isotopic anomalies in NDs, our results are useful in the estimation of the composition of end-members of the P3 and HL components and in studies of isotopic fractionation of released NGs.

Our paper is structured as follows. In the methodology, we introduce our MD approach, which includes an Arrhenius framework to map the simulation temperatures onto their experimental equivalent. Our first results are qualitative, using visualization and a movie of the implantation and thermal release processes to highlight the differences between He and Xe. The second set of results involves robust statistical analysis of a large number of simulations spanning many implantation energies, implantation directions, and annealing temperatures. The final section uses the Arrhenius approach to make a direct comparison between meteoritic data and the simulations.

2. Methodology

2.1. Simulation Methods

The MD simulations are performed using the environment-dependent interaction potential (EDIP; Marks 2000) for C–C interactions in combination with the standard Ziegler–Biersack–Littmark (ZBL) potential (Ziegler et al. 1985) to describe close approaches (see Aghajamali 2020 for further details). EDIP has proved itself to be highly transferable (de Tomas et al. 2016, 2019) and has previously been successful in simulations of many different forms of C, such as diamond (Fairchild et al. 2012; Marks et al. 2012; Buchan et al. 2015; Regan et al. 2020), carbon onions (Lau et al. 2007), amorphous carbon (Best et al. 2000; Marks et al. 2002), nanotubes (Suarez-Martinez & Marks 2012), peapods (Suarez-Martinez et al. 2010), and carbide-derived carbons (de Tomas et al. 2017, 2018). One aspect of EDIP that is still being developed is the ability to describe hydrogen; accordingly, none of the simulations presented here include hydrogen. To describe He–C and Xe–C interactions, we use a Lennard–Jones (LJ) potential coupled with the ZBL potential, with interpolation between the two interactions controlled using Fermi-type switching functions as described in Buchan et al. (2015) and Christie et al. (2015). This combination of the LJ and ZBL potentials is identical to our recent paper (Fogg et al. 2019), where we studied the modification of NDs by Xe implantation. Full details are provided in Appendix A.

We perform our simulations using an in-house MD package. Implantation simulations are carried out in an NVE ensemble (meaning that the number of particles (N), volume (V), and energy (E) are conserved quantities) using Verlet integration and a variable time step (Marks & Robinson 2015). Annealing simulations are performed in an NVT ensemble using a velocity-rescaling thermostat. Periodic boundary conditions are not employed. All simulations use a 4999-atom ND with a diameter of 3.9 nm; this size is similar to synthetic and meteoritic NDs (Huss & Lewis 1994a, 1994b; Koscheev et al. 2001). (Meteoritic NDs are characterized by a lognormal size distribution with a median of 2.6 nm and maximal sizes up to 10 nm or more; Lewis et al. 1987.)

As shown in Figure 1, the ND has a truncated octahedral form, which is the stable geometry for dehydrogenated NDs (Barnard & Zapol 2004). The {100} faces of the ND are reconstructed in a 2 × 1 manner to eliminate dangling bonds. The local environment, bond length, and bond angles in this reconstructed geometry are shown in Figure 1(c) and compared with density functional theory (DFT) data from De La Pierre et al. (2014). Most of the bond lengths are very similar (within 0.04 Å), the only exception being the sp²–sp² bond, which is 0.16 Å larger with EDIP. The {111} faces are not reconstructed, as this is the stable configuration with EDIP. We note that EDIP does not predict the 2 × 1 π-bonded chain reconstruction (Pandey 1982) of the {111} faces; EDIP favors the 1 × 1 geometry by 3.6 J m⁻², whereas DFT calculations (De La Pierre et al. 2014) favor the 2 × 1 geometry by 0.93 J m⁻². We do not expect the 1 × 1 geometry preferred by EDIP impacts on our results because both reconstructions graphitize in a similar manner and produce the same final structure and orientation relative to the rest of the ND.

The coordinates of the ND were generated using a methodology described in our recent paper studying Xe implantation (Fogg et al. 2019). The key idea is to cut the ND out of an infinite crystal using clipping planes in the 〈100〉 and 〈111〉 directions. Using the notation of Fogg et al. (2019), the ND used in this work was generated with d₁₀₀ = 20 Å and d₁₁₁ = 17 Å. Prior to implantation, the ND is equilibrated at 300 K.

Implantation simulations are 1 ps in length, sufficient to model the ballistic phase of ion implantation into the ND. To generate a wide variety of implantation conditions, the initial position and direction of the implanted species (either He or Xe) were systematically varied. Following our previous studies (Buchan et al. 2015; Christie et al. 2015; Shiryaev et al. 2018; Fogg et al. 2019), the initial position is taken from a 25-point solution to the Thomson problem (Thomson 1904), which distributes coordinates uniformly on a sphere (the Cartesian coordinates are provided in Aghajamali 2020). Some implantations are performed directly toward the center of the ND, while others are directed slightly (up to 10 Å) away from the center of mass. After the system equilibrates, the implantation depth is computed relative to the nearest crystallographic face, {100} or {111}. For the He implantations, the mass was 4 amu, while the Xe implantations use a mass of 133 amu, the same as in Shiryaev et al. (2018) and Fogg et al. (2019) and slightly higher than the average isotopic value of 131.3 amu. Additional Xe simulations use masses of 124 and 136 amu, corresponding to the lightest and heaviest stable isotopes.

Annealing simulations extending up to 1 ns are performed to study the thermal release process. Coordination analysis is performed by counting the number of nearest neighbors within a cutoff of 1.85 Å. For the purposes of analysis and visualization, C atoms are considered to be sp, sp², and sp³ hybridized if they have two, three, and four neighbors, respectively. Visualization is performed using the OVITO package (Stukowski 2010).

2.2. Mapping to Experimental Temperatures

Generally speaking, the timescale of MD simulations is on the order of nanoseconds, around 13 orders shorter than the experimental annealing time. This enormous difference complicates the comparison between simulations and experiments, as thermally activated events will be suppressed in the simulations due to the vastly shorter time. Many different solutions to the MD timescale problem have been proposed (Voter 1997, 1998; Sørensen & Voter 2000), but here we employ a simple temperature-acceleration approach that we have used successfully on other C systems (de Tomas et al. 2017). The first step is to assume Arrhenius behavior and a single activation energy, which corresponds to the relation

$\begin{eqnarray}&&f=A\,\exp (-{E}_{a}/{k}_{B}T),\end{eqnarray} \tag{ 1 }$

where f is the frequency of events, A is the attempt frequency, T is the temperature, k_B is Boltzmann's constant, and E_a is the activation energy. The correspondence between the experimental and simulation temperature is determined by equating the time-frequency product (i.e., f_expt × t_expt = f_sim × t_sim) to ensure that, for the same activation energy, the same number of events occurs in both simulation and experiment. This yields the following expression:

$\begin{eqnarray}&&{T}_{\mathrm{sim}}=-\displaystyle \frac{{E}_{a}}{{k}_{B}}\times {\left[\mathrm{log}\left(\displaystyle \frac{{t}_{\mathrm{expt}}}{{t}_{\mathrm{sim}}}\right)-\displaystyle \frac{{E}_{a}}{{k}_{B}{T}_{\mathrm{expt}}}\right]}^{-1}.\end{eqnarray} \tag{ 2 }$

This equation links the experimental temperature and time (T_expt and t_expt) with those of the simulation (T_sim and t_sim), with the only parameter being the activation energy. To determine E_a, we make use of experimental data on graphitization of NDs, where it is known that ∼1500°C is required to convert NDs into a carbon onion (Kuznetsov et al. 1994a, 1994b; Tomita et al. 2002; Qiao et al. 2006; Xiao et al. 2014). To identify the corresponding temperature on the MD timescale, we performed a set of 1 ns simulations at different annealing temperatures and found that onionization of the ND occurs at around 3500 K (see Figure 2). Using these two temperatures and suitable times (t_sim = 10⁻⁹ and t_expt = 3600 s; Kuznetsov et al. 1994a, 1994b; Tomita et al. 2002; Qiao et al. 2006; Xiao et al. 2014), we obtain via Equation (2) a value of E_a = 9 eV. Having determined the activation energy, we can employ Equation (2) to map any simulation temperature to its experimental equivalent using the relation shown in Figure 3. For example, the pink lines in the figure show that a simulation temperature of 4500 K is equivalent to an experimental temperature of ∼1730°C, and a simulation at 2000 K is equivalent to an experiment at ∼1000°C. The latter also serves as a convenient dividing line between the main peaks of the Xe-P3 and Xe-HL components (Huss & Lewis 1994a, 1994b). To provide a sense of scale, Figure 3 also shows calibration curves for two other activation energies. Note that due to the logarithm term in Equation (2), the value of E_a has minimal sensitivity to the choice of t_expt. This point is discussed in detail in de Tomas et al. (2017). Finally, we note that throughout this paper, we adopt the convention that experimental temperatures are given in Celsius and simulation temperatures are in Kelvin. This helps conceptually separate the two quantities, which otherwise might be confused with one another.

**Figure 2.** Cross-sectional snapshots (1 nm slices) of pristine ND after annealing for 1 ns at various temperatures. Panel (f) shows that at 3500 K, the experimentally observed result, namely a carbon onion, is obtained. Red, green, and blue circles denote sp–C, sp²–C, and sp³–C hybridization, respectively. A full discussion of the onionization process is provided in Appendix B.
Download figure:
Standard image High-resolution image

**Figure 3.** Calibration curve between simulation and experimental temperatures via the Arrhenius approach. The orange dashed line indicates data for the onionization of an ND used to determine the value of E_a = 9 eV. The pink lines indicate two examples of the mapping process as described in the text.
Download figure:
Standard image High-resolution image

3. Results and Discussion

3.1. Individual Implantation and Thermal Release Events

Representative examples of the implantation and thermal release processes for He and Xe are shown in Figure 4. Panels (a) and (c) show that the implantation processes for He and Xe differ substantially, with the higher mass and size of Xe having a large effect. For He, only a relatively modest energy (∼100 eV) is needed to implant the atom into the center of the ND, and the He undergoes multiple deviations along the implantation trajectory (pink line), since the C atoms are three times heavier. In contrast, the implanting Xe simply slows down, and these are the C atoms that move. Around ∼800 eV of kinetic energy is needed to implant into the central region, similar to the value of circa 700 eV used by Koscheev et al. (2001) in their experiments.

Figure 4. Typical implantation (left panels) and thermal release (right panels) processes involving He (panels (a) and (b)) and Xe (panels (c) and (d)). The pink lines indicate the He and Xe trajectories, and the He and Xe are shown as orange circles. The color codings and cross-sectional slice details are the same as in Figure 2. An animation of the time evolution of the He and Xe release processes is available in the online journal. The animation is sequenced first showing the He release process (annealing time from 0–162 ps) followed by the Xe release process (annealing time from 0–433 ps).

(An animation of this figure is available.)

Download figure:

Video Standard image High-resolution image

Once the implanted ND system has equilibrated, the entire cluster is heated to find the temperature at which the NG atom escapes. Examples of this process are shown in panels (b) and (d). For both species, thermal release occurs at relatively high simulation temperatures: 3000 K for He and 3500 K for Xe. As seen in Figure 2, these temperatures are sufficient to transform portions of the ND into a onion-like structure. Due to its small size, the He is able to escape before the ND has completely onionized. The complexity of the process can be appreciated in the first part of the online movie (Figure 4), which shows how the He travels along multiple 〈110〉 channels and traverses a substantial fraction of the ND before escaping. In the case of Xe, the ND is fully transformed into a carbon onion, and the release process is more difficult, as the Xe must pass through the graphitic shells. The pink trajectory in panel (d) shows that the release process involves two local minima in which the Xe jiggles back and forth many times before escaping. An animation sequence of this process is provided in the second part of the online movie (Figure 4).

To determine the temperature at which He and Xe are released, we perform 1 ns annealing simulations at many different temperatures and monitor the distance between the NG species and the center of the mass of the system. An example of our methodology is provided in Figure 5, which shows this quantity (r_c.m.) for Xe as a function of time for three different annealing temperatures using the implanted ND in Figure 4(c) as the starting structure. At 3000 K (blue line), the Xe remains the same distance from the center of mass. This occurs because the Xe becomes trapped at the interface between the ND core and onionized outer layers (see Figure 2(e)). At 3250 K (violet line), the Xe moves toward the surface of the ND during the first ∼170 ps of annealing, but afterward, it is trapped between the graphitic sheets. Note that at around 30 and 160 ps, there are two significant surges in r_c.m. where the Xe atom passes through the graphitic sheets. At 3500 K (pink line), the Xe is initially trapped between the sheets, but after moving between them, it finally exits through a gap in the outer shell at 432 ps, producing a sharp jump in r_c.m..

**Figure 5.** Distance between Xe and the center of mass of the ND as a function of annealing time for three different annealing temperatures. The structure in Figure 4(c) is the starting point for all three simulations.
Download figure:
Standard image High-resolution image

3.2. Statistical Analysis

To collect a statistically significant data set, a large number of implantation and thermal release simulations need to be performed. The first step is to perform simulations that implant He and Xe into a wide variety of configurations within the ND. For each species, a total of 1875 implantations were performed: 25 different energies, each with 75 different initial conditions (i.e., directions and/or impact parameters). A summary of the resultant implantation depth and incorporation probability is shown in Figure 6. Panels (a) and (b) show the implantation depth of He and Xe as a function of implantation energy, where each dot indicates an implantation event whereby the NG species remains with the ND. Impacts where the He or Xe leave the ND are not shown. The fewest dots occur for high-He (where He passes through the ND) and low-Xe (where Xe is reflected) implantation energies. The distribution of depths differs considerably between the two species. For He, the implantation depth spans the maximum possible range, extending from the surface to the center of the ND, and the implantation depth is uncorrelated with energy. In contrast, low energies result in only shallow implantation of Xe, while nearly 600 eV is required to span the full range of depths. High-Xe implantation energies still produce a large number of shallow implantation depths, which is perhaps due to the small size of the ND.

Figure 6(c) quantifies the incorporation probability of He and Xe as a function of implantation energy. The probabilities for the two species are strikingly different and driven by the mass difference relative to C. These data show that low energies are optimal for He implantation, while efficient Xe implantation requires many hundreds of eV. Noting that He is the most abundant NG in meteoritic NDs (Huss 2005), these data impose constraints on the astrophysical conditions for He incorporation via implantation.

The second step in performing the statistical analysis uses the coordinates associated with each dot in Figure 6 as the starting structure of thousands of thermal release simulations. All simulations run for 1 ns, except for when r_c.m. indicates that release has occurred; in such cases, the simulation is terminated. If thermal release does not occur, the simulation is rerun at a higher temperature. Typically, the temperature increment between successive simulations is 100–250 K. The raw simulation data showing the relationship between the release temperature of the NG and the implantation depth are shown in Figure 7. The dots and error bars indicate the precision, meaning that if an NG atom releases at a temperature T₂ but not at a lower temperature T₁, then the dot denotes (T₁ + T₂)/2, and the error bar indicates the range [T₁: T₂]. Panel (b) shows that the Xe release temperatures cluster into two groups and correlate with the depth of the implanted atom. The average temperatures for these two clusters are indicated by the green lines, with a dividing line of 2000 K used to separate the groups. Varying this number by several hundred Kelvin makes little difference to the averages. For He, no clustering occurs, and the release temperature gradually increases with implantation depth. In this case, the green line indicates the average temperature for the full data set.

3.3. Comparison with Experiment

Figure 7(a) shows that a minimum of 1300 K is required to release He, and by 4000 K, all He is released. In contrast, the Xe data in panel (b) span a broader range, with a minimum release temperature of only 400 K and a maximum of nearly 4500 K, at which point the ND is effectively destroyed (see Figure 2(h)). The observation of two temperature clusters for Xe and a single broad distribution for He is in good qualitative agreement with the meteoritic ND observations. To map the simulation temperatures onto their experimental equivalents, we employ the Arrhenius approach explained in the methodology, using an activation energy of E_a = 9 eV. By histogram binning the simulation data in Figure 7 and applying the transformation in Equation (2), we obtain a data set shown as thick blue lines in Figure 8. On the same scale, we show experimental data for the Orgueil meteoritic NDs extracted from Huss & Lewis (1994a, 1994b). The agreement is remarkable, with the simulations reproducing all of the main meteoritic ND characteristics, including (i) the unimodal versus bimodal character, (ii) the position of the peaks, (iii) the widths of the distributions, and (iv) the maximum and minimum release temperatures. This is the first time that MD simulation has predicted these important effects.

**Figure 8.** Measured and simulated ND thermal release patterns for (a) ⁴He and (b) ¹³³Xe. Experimental values (red lines) are taken from data from the Orgueil meteorite ND (Huss & Lewis 1994a, 1994b), while the MD values (blue lines) are from this work.
Download figure:
Standard image High-resolution image

The high level of experimental detail reproduced by the simulations provides post hoc justification for the Arrhenius approach, confirming that the presumption of a dominant activation energy is reasonable for this class of problem. All of the predicted temperatures are correctly positioned relative to their experimental equivalents; this includes the onset of He and Xe release, the position of the release peaks, and even the upper limit at ∼1800°C, which corresponds to destruction of the ND and the release of any remaining gases. Regarding Figure 8(b), it is important to note that the simulations cannot predict the relative height of the peaks in the bimodal distribution, since this is a function of the balance between shallow and deeply implanted Xe that is not known.

Having reproduced the essential characteristics of NG release, we can address some of the most fundamental questions related to meteoritic NDs. Specifically, we can examine the atomistic origin of the low- and high-temperature peaks of Xe release and understand why He exhibits unimodal behavior. Considering first the question of Xe, the raw simulation data in Figure 7(b) provide clues as to why release occurs at two distinct temperatures. The low-temperature peak is seen to be strongly correlated with proximity to the surface, with 4.5 Å being the critical depth beyond which low-temperature release is impossible. In contrast, high-temperature release is possible for all implantation depths, and even for depths around 1 Å, there are two distinct temperature release populations. The origin of this behavior is that the surface of the ND progressively graphitizes with temperature, with the {111} faces transforming prior to the {100} faces (see Figure 2). As a result, Xe located close to the {111} faces can become trapped by the developing graphitic layers, and once the layer has fully formed, the Xe is too large to easily diffuse through the hexagonal graphene-like network. This effect is quantified in Figure 9, which replots the Xe data in Figure 7(b) for shallow depths. Panel (a) color codes each configuration according to the closest crystallographic face, while panel (b) shows a histogram of the thermal release temperatures. It is apparent that the low-temperature peak contains roughly equal contributions from Xe near the {100} and {111} faces, while the high-temperature component is dominated by Xe close to the {111} faces. Placing these observations in a gas release context, we can assert that a high-temperature peak (such as the main peak of Xe-HL) is associated with either deeply buried Xe or a shallow burial near the {111} faces, while the ions released at low temperatures (in particular, e.g., Xe-P3) sit just a few angstroms from the surface and have no crystallographic preference.

**Figure 9.** (a) Classification of shallow implantation Xe data in Figure 7 according to proximity to {100} faces (orange circles) or {111} faces (green circles). (b) Histogram analysis of the data in panel (a).
Download figure:
Standard image High-resolution image

The other major result in Figure 7 is the prediction of the unimodal release distribution for He. The animation sequences (online Figure 4) reveal that this behavior can be linked to the nature of the He defect within the ND. In particular, the simulations show that He prefers a tetrahedral interstitial (T-site) in the diamond lattice. Manual inspection of a large number of implantation events shows that the He is always in an interstitial position; substitutional defects and other lattice damage are not observed. The geometry of the T-site defect is shown in Figure 10(a), with the four orange lines between He and C highlighting the tetrahedral symmetry. The interstitial He diffuses among the 〈110〉 channels in the lattice, passing through a hexagonal interstitial (H-site) transition state en route to another T-site. This behavior means that in the low-temperature range, the He effectively performs a random walk around the ND. Even if the {111} faces have graphitized, which typically occurs at around 1500 K in the simulations (∼800°C in the experiments), the He cannot escape through the graphitic layer. Only once the {100} faces begin to graphitize does the He atom escape. These observations elegantly explain why the onset of He release is so much higher than Xe. Additionally, the smaller size of He means that it is more mobile than Xe, which in turn explains why the peak release temperature for He is lower than that of the corresponding Xe high-temperature peak.

**Figure 10.** (a) T-site of He and (b) Xe_s–V defect of Xe in ND. Color-coding details are the same as in Figure 4; He and Xe are shown as orange circles.
Download figure:
Standard image High-resolution image

The analysis of He migration in our MD simulations is supported by DFT calculations of NG defects in bulk diamond (Goss et al. 2009), which similarly conclude that He forms an interstitial defect on the T-site and diffuses via the H-site. The distances predicted by EDIP compare well with those from DFT. For the T-site, the EDIP (DFT) values are 1.61 (1.58) Å for C–He, 1.60 (1.57) Å for the three extended C–C bonds, and 1.52 (1.50) Å for the single compressed C–C bond. For the H-site, the EDIP and DFT bond lengths are identical to two significant figures: 1.57 Å for the C–He distance and 1.63 Å for the C–C bonds. The prediction for the relative energetics is only qualitative, with EDIP favoring the T-site by 0.14 eV, whereas with DFT, the T-site is preferred by 2.3 eV. This large difference means that He is mobile at room temperature in the MD simulations. The observation that the simulated He release temperature agrees well with experimental observations shows that He release is driven by the thermally induced modification of the ND.

For the case of Xe defects in diamond, DFT calculations by Goss et al. (2009) and Drumm et al. (2010) show that the preferred structure is a substitutional site involving an adjacent vacancy (Xe_s–V), with the Xe placed at the midpoint (Figure 10(b)). In the relaxed structure, the Xe is located equidistant from six carbon atoms; with EDIP, this distance is 2.25 Å, as compared to 2.15 Å in the DFT calculation of Drumm et al. (2010). Returning our attention to the simulation data in Figure 8(b), we can assign the low-temperature Xe release peak to an Xe_s–V defect near the surface. Since this defect compromises the stability of the diamond surface, the Xe is able to escape at modest temperatures of a few hundred degrees Celsius.

The DFT study by Goss et al. (2009) also studies other NGs and divides them into two groups: (i) He and Ne, which occupy the interstitial T-site and diffuse via the H-site, and (ii) Ar, Kr, and Xe, which occupy substitutional sites with vacancies. This distinction between the interstitial T-site and the substitutional vacancy provides a plausible explanation for the meteoritic ND data, where He and Ne have unimodal release peaks, while Ar, Kr, and Xe have bimodal distributions. To the best of our knowledge, this connection between defect type and temperature release behavior has not previously been made.

3.4. Isotopic Effects

Investigation of NG release kinetics from synthetic NDs reveals marked mass-dependent fractionation for all elements (Koscheev et al. 2001). Namely, the gases released at high temperatures become progressively enriched in heavy isotopes. At 1400°C, the magnitude of this effect reaches 0.99% per amu, and it is 1.44% per amu for Xe and Kr, respectively; even larger values are observed for lighter elements (Huss et al. 2008). Our MD approach naturally explains this observation.

Figure 11 shows thermal release data for ¹²⁴Xe and ¹³⁶Xe. As seen earlier, the release profile is again bimodal. For the low-temperature peak, the release temperature is similar across both isotopes, spanning a narrow range between 950 and 980 K, but for the high-temperature peak, there is a large isotopic effect, with the simulation release temperature varying from 3180 K for ¹²⁴Xe to 3485 K for ¹³⁶Xe. Figure 12 plots the mean release temperature for the high-temperature peak for three Xe isotopes, with the left axis showing the simulation temperature and the right axis showing the equivalent experimental value; error bars denote the standard error in the mean. Between the lightest and heaviest isotopes, the predicted temperature difference is over 80°C.

The isotopic effect in Figure 12 can be plausibly attributed to the effect of mass on the vibrational frequency during thermal release. An alternative explanation focusing on implantation is less attractive, since a spectrum of implantation energies is employed, and there is no obvious reason why a mass difference would generate different types of defects. Regarding the thermal release process, Figure 4(d) (and the online movie) shows that the ¹³³Xe atom escapes after an extended period of constant "jiggling." For about 350 ps (between t = 50 and 400 ps), the Xe is trapped between two graphitic shells and can be seen to move back and forth many times before eventually escaping. This observation helps explain why the lighter Xe isotopes release at lower temperatures, as the smaller mass implies a higher vibrational frequency and hence faster reaction rate.

4. Discussion and Implications for Meteoritic NDs

All calculations performed in our work are relevant for an ND grain with a well-defined morphology. It is shown that the proximity of the implanted atom to a given crystallographic face—{100} or {111}—can be important for the temperature of the gas release. Here an important question arises: what is the minimal size of an ND grain that develops crystallographic faces? For meteoritic NDs, the situation is not obvious, since all imaged grains appear rounded, and grains found in situ in meteorites (Garvie & Buseck 2006) are always associated with disordered sp²–C; spectroscopic data EELS (Garvie 2006) and XANES (Shiryaev et al. 2011) also indicate a significant amount of sp² bonding. The origin of this shell is unclear; it may result from severe thermochemical treatment employed for separation of the NDs from the meteorite and thus represent damaged diamond surfaces, and/or it may be a secondary material precipitated from the acids. Alternatively, it may reflect a deviation of the morphology of the small grains of meteoritic NDs from the assumed cuboctahedron shape. However, high-resolution TEM studies complemented by DFT modeling of synthetic ND grains showed well-developed faceting even for grains less than 2 nm, although the sp^2+x shell is always present (Stehlik et al. 2016; Chang et al. 2018). We note also that very small (1.6 nm in diameter) meteoritic NDs may possess intense photon emission of the Si–V defect (Vlasov et al. 2014), a fact indicating fairly perfect ordering of the nanograin lattice and little contribution of surface states.

An important related issue is the presence and eventual influence of surface functional groups on ND grains on the results of our calculations. As noted above, the environment of formation of NDs in space remains debatable, and chemical treatment complicates the discussion of their surface chemistry. In the case of NDs in outer space, the presence of significant concentrations of surface-bound N and O is not very plausible. However, hydrogen is almost certainly present, as suggested by both direct astronomical observations of C–H vibrations at 3.43 and 3.53 μm assigned to NDs (Guillois et al. 1999; Van Kerckhoven & Waelkens 2002) and the unusual isotopic composition of hydrogen in meteoritic NDs, possibly inherited from their formation (Virag et al. 1989). However, the extent of the surface hydrogenation is yet uncertain; see the discussion in Shiryaev et al. (2015). The influence of surficial hydrogen on the results of our calculations requires additional investigation, but we recall again that even for 1.6 nm ND grains, the influence of surface states appears to be minor.

The current work provides an atomistic explanation of a remarkable phenomenon: the release pattern of heavy NGs implanted into ND grains is bimodal and unimodal for He. We show that the release patterns are satisfactorily explained by such peculiarities of the behavior of the implanted ion in a nanograin as implantation depth, in-grain diffusion, and thermal evolution of the grain structure heated by the ion impact. Besides general interest for understanding size effects in the interaction of ions with matter, this work is relevant to carbon dust evolution in space. Two main scientific tasks are important: (i) the origin of isotopically different NG components and (ii) the presence of different ND populations, differing in origin, size, and/or morphology.

Xenon belonging to two of the three main NG components in NDs—P3 and HL—releases in two main peaks. It is usually assumed that the low-temperature peak is dominated by the P3 gases and the high-temperature one by the HL gases. The P3 gases appear to reside in traps with a broad range of activation energies (Huss & Lewis 1994b) and are readily lost during thermal metamorphism of the parent bodies of meteorites. Investigation of size-separated NDs revealed that the ratio of P3/HL decreases in small ND grains, but the absolute concentrations of both components increase with the size (Verchovsky et al. 1998; Fisenko et al. 2004).

Although the present work does not reject the possibility of association of the P3 component with species adsorbed in surficial sp²-containing shells, it shows that monoenergetic implantation alone nicely explains both components. In this case, both low- and high-temperature peaks should be present in a given irradiated population of NDs, implying that if implanted, the P3 and HL gases should contribute to the whole release pattern. In fact, Huss & Lewis (1994b) and Huss et al. (2008) suggested that an important fraction of HL gases—"labile HL"—is released roughly simultaneously with the P3 component and thus is lost during thermal metamorphism of the parent meteorite and/or due to heating in the nebula. In contrast, Koscheev et al. (2001) and Fisenko & Semjonova (2020) proposed that it is the P3 component, which dominates both low- and high-temperature release peaks. Taken together with the hypotheses discussed above, our work indicates that calculation of the Xe-HL isotopic composition should include corrections for both the high-temperature P3 release and preferential enrichment in heavy isotopes at these steps due to isotopic fractionation. Derivation of the isotopic composition of the "pure" HL component nevertheless remains difficult, since it is likely a product of more than one nucleosynthetic process.

The claim that the P3 component should have a high-temperature component also requires the assumption that isotopes of other elements should be implanted. In particular, the contribution of extinct iodine isotopes—¹²⁹I and, perhaps, ¹²⁸I—is invoked for the explanation of the Xe-P3 composition (see Gilmour et al. 2005, 2016 for details). Contribution of surface-bound ¹²⁹I can be ruled out, since the energy of the β-decay is far too small to drive the Xe recoil into an ND grain. Subsequently, implantation of radioactive iodine isotopes should have occurred. In recent work, Gilmour et al. (2016) suggested that a fraction of Xe released at high temperatures was also produced in situ from iodine. Although we have not explicitly analyzed the behavior of implanted iodine, we believe that the similarity in masses between I and Xe should lead to a qualitatively similar behavior of these elements.

The decrease of the P3/HL ratio in smaller ND grains is, however, more difficult to explain unambiguously. In small grains, almost all implants will be in close proximity to the surface, and direct application of our modeling would suggest a higher fraction of the low-temperature release. However, distortion of the shape of the smallest grains (Chang et al. 2018) influences the curvature of the facets and thus may strongly influence the diffusion paths of an implant. Variations of the curvature, as well as interaction of the implant with the sp^2+x surface shell, may also be responsible for the broad range of bonding energies of the P3 traps. In the same time, there is an ample evidence that NDs from every meteorite represent a mixture of several independent populations of grains (Russell et al. 1996; Verchovsky et al. 1998), which may have very different sources.

Examination of Figures 9 and 11 might give some hints for the origin of the high-temperature Xe-P6. The highest release temperatures were observed for a few events of very deep Xe implantation. Examination of size-separated NDs suggests that the P6 component appears to be present mostly in coarse grains (Gilmour et al. 2005). Our calculations addressed only relatively small ND grains with a diameter of 3.9 nm. Examination of larger grains is computationally very expensive and was not performed, but it is reasonable that for a given implantation energy and larger grain size, a larger fraction of the ions will stop at greater depths. The release temperature of these ions will also tend to cluster at higher values.

Last but not least, a result of our calculations is a confirmation of the stability of a Xe–V point defect in an ND grain. In annealed ion-implanted macroscopic diamond, this defect possess rather strong infrared photoluminescence (lines at 793–794 and 811–813 nm) when excited by visible light (Deshko & Gorokhovsky 2010, and references therein). In nanoparticles, the annealing required for conversion of the implant to the favorable configuration occurs during the implantation process itself. Together with the luminescence of the silicon-vacancy (Si–V) defects (Shiryaev et al. 2015), the Xe–V complex opens a new possibility for eventual observation of NDs in stellar environments.

5. Conclusion

In this work, we developed an MD approach to address the important unresolved question of the origin of the unimodal and bimodal thermal release patterns of NGs from meteoritic NDs. Our technique employs a large number (circa 10⁴) of small MD simulations to calculate the thermal release pattern of implanted NGs in ND and provides detailed atomistic insight. We reproduce the known experimental profiles for He and Xe and propose that the position and structure of the NG defect are responsible for the unimodal and bimodal patterns, respectively. In the case of Xe, we show that the low-temperature component is associated with shallow implantation on the {100} faces, while the Xe-HL component is a mixture of deeply buried defects and shallow implantation on the {111} faces. Although our simulation methodology does not capture every single experimental detail, it appears that a relatively straightforward description of the system captures the key features. Areas of extension in future work include the addition of hydrogen, modeling of the sp^2+x shell, and improved potential to describe the energetics of defect migration and the {111} surface reconstruction.

The MD methodology is a natural fit for cloud-computing facilities and can be easily extended to other systems where NGs are found in presolar grains. In NDs, obvious directions for future work include unimodal release in ³He, ²⁰Ne, ²¹Ne, and ²²Ne and bimodal behavior of ³⁸Ar, ³⁶Ar, ⁸⁴Kr, and ⁸⁶Kr. Success in replication of the isotopic effects may serve as a tool to computationally investigate the implantation processes of short-living isotopes, for example, spanning the whole range of known Xe isotopes from ¹⁰⁸Xe to ¹⁴⁸Xe.

Aside from NDs, there are several other materials relevant for interstellar dust studies where our MD methodology can be applied. These include silicon carbide (SiC) and graphite (Amari et al. 1990, 1993; Ott 1993, 2007, 2014; Davis 2011). These presolar materials have high melting points and therefore are suitable for our Arrhenius-based approach. The ballistic stage of ion implantation in nano- and micron-sized C and SiC grains was studied using a Monte Carlo–based approach (SRIM; Verchovsky et al. 2003). However, these calculations cannot address grain heating, annihilation of defects, and many other effects. Therefore, whereas the range of ions is rather accurately predicted, the overall impact of the implantation on nanograins remains elusive. A number of successful atomistic simulation studies address high-temperature effects in SiC (Fujisawa et al. 2019), graphite, and carbide-derived carbons (de Tomas et al. 2017, 2018). These structures provide a natural starting point for further atomistic simulations to reveal the astrophysical secrets of NGs in presolar grains.

The main findings of our work are summarized as follows.

(i)
Upon heating, monoenergetic Xe ions implanted into an ND grain will be released in two distinct peaks at markedly different temperatures.
(ii)
Monoenergetically implanted He ions will be released as a single broad peak with an onset shifted to higher temperatures relative to Xe.
(iii)
Isotopic fractionation favoring heavy isotopes accompanies high-temperature steps.
(iv)
Differences in the behavior of implanted NGs are explained by diverse types of formed lattice defects and diffusion paths.
(v)
The exact release pattern of the implanted NGs depends on the morphology of the grain, the position of the implant in the grain, and the peculiarities of the impact.
(vi)
Our results suggest that the implantation scenario can explain both the P3 and HL components of Xe in NDs. In addition, both components should possess low- and high-temperature release peaks. The absence of a low-temperature peak for the HL and, possibly, the P6 component is due to postimplantation annealing of the grains.

We gratefully acknowledge many helpful discussions with Ulrich Ott. N.A.M. acknowledges fellowship FT120100924 from the Australian Research Council. Computational resources were provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia.

Appendix A: He–C and Xe–C Interactions

As described in the methodology section of the main text, for the He–C and Xe–C interactions, we use the standard ZBL potential coupled with an LJ potential. The LJ potential has the form

$\begin{eqnarray}&&{U}_{\mathrm{LJ}}(r)=4\varepsilon \left\{{\left(\displaystyle \frac{\sigma }{r}\right)}^{12}-{\left(\displaystyle \frac{\sigma }{r}\right)}^{6}\right\},\end{eqnarray} \tag{ A1 }$

where for He–C interactions, ε_He−C = 0.0013 eV and σ_He–C = 2.98 Å (Nguyen et al. 2008), and for Xe–C interactions, ε_Xe−C = 0.0114 eV and σ_Xe–C = 3.332 Å (Simonyan et al. 2001). Following Buchan et al. (2015) and Christie et al. (2015), the total interaction energy is expressed as

$\begin{eqnarray}\begin{array}{rcl}{U}_{\mathrm{Total}}(r) & = & \left[{U}_{\mathrm{ZBL}}(r)\times f(r+\delta )\right]\\ & & +\left[{U}_{\mathrm{LJ}}(r)\times (1-f(r-\delta ))\right],\end{array}\end{eqnarray} \tag{ A2 }$

where f(r) is a Fermi-type switching function, δ_He–C = 0.09 Å, and δ_Xe–C = 0.07 Å (Fogg et al. 2019). The Fermi function is given by

$\begin{eqnarray}&&f(r)={\left[1+\exp \left({b}_{F}(r-{r}_{F})\right)\right]}^{-1},\end{eqnarray} \tag{ A3 }$

where b_F controls the sharpness of the transition and r_F is the cutoff distance; these parameters are chosen manually to ensure smoothness. For Xe–C, we use b_F = 8 Å and r_F = 2.7 Å, the same as our recent implantation studies (Shiryaev et al. 2018; Fogg et al. 2019), while for He–C, we use b_F = 9.2 Å and r_F = 2.2 Å. Figure 13 illustrates the He–C and Xe–C interaction energy covering the ZBL and LJ regimes. Further discussion and full details of the interpolation process are provided in Buchan et al. (2015) and Christie et al. (2015).

**Figure 13.** Pairwise interaction energies (blue line) for (a) He–C and (b) Xe–C. At close approach, the interaction is pure ZBL (orange line), while at distances around equilibrium and greater, an LJ expression is used (green line).
Download figure:
Standard image High-resolution image

Appendix B: Annealing of Pristine ND

As described in the main text, the simulation temperatures are mapped onto their experimental equivalent using the Arrhenius relation. Calibration is performed using the graphitization of NDs into C onions as a reference point. Figure 2 in the main text illustrates cross-sectional snapshots of the annealed ND for eight different temperatures. Panel (a) shows that at 1000 K, there is insufficient temperature to graphitize the ND on this timescale, while at 1500 and 2000 K (panels (b) and (c)), some of the bonds on the {111} faces are broken, and graphitic shells form. The reason this occurs is that sp³–C hybridized atoms near the surface can easily rearrange into graphite with sp²–C hybridization (Kuznetsov et al. 1994b; Wang et al. 2000; Pantea et al. 2002; Barnard et al. 2003). At higher annealing temperatures of 2500 and 3000 K (panels (d) and (e)), graphitization proceeds further inward, creating structures with a diamond core surrounded by a graphitic shell. This core–shell structure is in good agreement with previous experimental and simulation studies (Kuznetsov et al. 1994a, 1994b; Tomita et al. 2002; Los et al. 2005; Bródka et al. 2006, 2008; Qiao et al. 2006; Ganesh et al. 2011). At 3500 K (panel (f)), a pure carbon onion (or concentric fullerene) is obtained, and only a small number of atoms (∼1%) have been lost via evaporation. At this point, we have identified the temperature at which the experimental result is reproduced on the timescale that is affordable in the simulation. At the slightly higher temperature of 4000 K (panel (g)), around 8% of the atoms evaporate, and the C atom develops a hollow core; similar structures are seen in Figure 8(e) of Lau et al. (2007). At the highest temperature of 4500 K (panel (h)), the ND is completely destroyed and more than 29% of the atoms are evaporated, creating a large, disordered fullerene.

The time evolution of the graphitization process is shown in Figure 14 for two different annealing temperatures. The top row shows data at 2000 K, where as time increases, the {111} faces graphitize first, followed later by the {100} faces. At the higher temperature of 3500 K (bottom), a similar sequence occurs, with the addition of an unraveling process, which is facilitated by the development of spirals.

Appendix C: Online Movie

Time evolution of He (first part) and Xe (second part) release processes from ∼3.9 nm ND during the high-temperature annealing (as discussed in Figure 4 in the main text). The thickness of the slice is 1 nm. The pink lines indicate the He and Xe trajectories, and the He and Xe are shown as an orange circle. The C atoms are shown as red, green, and blue circles with sp–C, sp²–C, and sp³–C hybridizations, respectively.

Molecular Dynamics Approach for Predicting Release Temperatures of Noble Gases in Presolar Nanodiamonds

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction