Rock-salt structure lithium deuteride formation in liquid lithium with high-concentrations of deuterium: a first-principles molecular dynamics study

All simulations were performed within the Vienna Ab Initio Simulation Package (VASP) [32–34]. The nuclei (Li and D) and core electrons (1s electrons of Li) were replaced by the potentials generated by the projector augmented-wave method [35]. We believe that the ionic potentials of Li and D used here are suitable for studying a wide range of concentrations of D in liquid Li because, as we shall see, these potentials yield trends in ionic densities and melting points with varying D concentration that are consistent with measurements. The Perdew–Burke–Ernzerhof XC functional (one of the most widely used generalized gradient approximations) was employed [24]. Periodic boundary conditions were chosen and the k-point meshes in the Brillouin zone were generated with the Monkhorst-Pack scheme [36]. The kinetic energy cutoff for the plane wave basis was set to 250 eV (with $18\times 18\times 18$ k-points) to converge the total energy of a two-atom unit cell of body-centered-cubic (bcc) Li to within 0.5 meV atom⁻¹ and 400 eV (with $9\times 9\times 9$ k-points) to converge the total energy of an eight-atom unit cell of rock-salt (RS) LiD to within 1.2 meV atom⁻¹. The bulk moduli were computed by expanding and compressing the equilibrium volume by up to 2% and then the resulting energy-volume curves were fit to Murnaghan's equation of state [37]. We chose to use 400 eV with one k-point (0,0,0) for all MD simulations; the lack of k-point sampling is appropriate for disordered systems such as a liquid. Born–Oppenheimer MD was used with the Nosé–Hoover thermostat [38] in the NVT ensemble (constant number of atoms N, volume V, and temperature T). Natural abundance masses of 6.941 and 2.014 amu were chosen for Li and D, respectively. The mass for the Nosé–Hoover thermostat was set to 2.0 a.u., which controlled the fluctuations in total energy to be no more than 1.08% of the average value. A time step of 0.2 fs, which is smaller than the 0.5 fs used for pure Li in our previous OFDFT-MD simulations [14], was used to ensure the accuracy of the Verlet integration of Newton's equations of motion in the MD simulations. All of our simulations were done at zero pressure, which approximates the pressure condition in the high-flux D experiment (10⁻⁷ torr) [39].

We chose a $4\times 4\times 4$ bcc cell (128 atoms with two atoms in each cubic unit cell) as an initial geometry for liquid Li and a $3\times 3\times 3$ RS cell (216 atoms with eight atoms in each RS cell) for LiD as an initial geometry for MD simulations. We then inserted four different concentrations of D atoms into liquid Li at temperatures ranging from 470 to 1143 K. Specifically, the ratios of D to Li atoms, defined as β, were 0.25, 0.50, 0.75, and 1.00. All six systems were simulated separately at six temperatures, namely, 470, 570, 686, 793, 943, and 1143 K. Thus, there are in total 36 systems that were studied by KSDFT-MD.

We performed the following five steps to generate simulation cells with the desired densities of liquid Li with inserted D atoms at different concentrations and temperatures. Each density was chosen upon the condition that the external pressure of the simulation cell fluctuated around zero. First, we obtained the density of Li at different temperatures (from 400 to 1143 K) using KSDFT-NVT-MD by following these procedures: (1) We took the initial density of Li from our previous OFDFT-MD simulations, in which the NPT ensemble (constant number of particles N, pressure P, and temperature T) was used to obtain the density of liquid Li at temperatures ranging from 200 to 725 K. For liquid Li at temperatures higher than 725 K, the initial density of liquid Li was obtained by linearly extrapolating the densities of liquid Li below 725 K. (2) Once the initial density of Li was set, we increased the density by 1% in each Li system to get a second density. For both densities, we ran a 10 ps NVT-MD trajectory and obtained the equilibrium pressures by sampling the last 5 ps of the trajectory. (3) We used the two new pressures and densities to linearly extrapolate to a new density at which the pressure should be zero. (4) We always used the two newest pressures to do the next extrapolation until the equilibrium pressure obtained was smaller than 0.2 GPa. We found that two to three iterations were enough to obtain the final density. (5) Properties such as pair distribution functions, static structure factors, and diffusion coefficients (see below) of liquid Li at different temperatures were then calculated from a 20 ps production run and compared to existing experimental data to ensure that we had correctly simulated liquid Li at zero pressure. The density of LiD at different temperatures was obtained by heating up a RS cell from 0 K to a target temperature. The initial density of LiD was obtained from ${{\rho}_{\text{LiD}}}(T)={{\rho}_{\text{Li}}}(T)\times \frac{{{\rho}_{\text{LiD}}}(0)}{{{\rho}_{\text{Li}}}(0)}$, where ${{\rho}_{\text{Li}}}(T)$ and ${{\rho}_{\text{LiD}}}(T)$ are the densities of Li and LiD at temperature T. Once the initial density of LiD was obtained, we followed exactly the same procedures used for Li to obtain the final density of LiD.

To obtain homogeneous solutions of D in liquid Li, the four concentrations of D atoms (LiD$_{\beta}$ with $\beta =0.25$, 0.50, 0.75, and 1.00) were inserted randomly by a random number generator into the final configuration of liquid Li from the 20 ps production run above, heated up to 1500 K, and then equilibrated for 2 ps. After this step, we found that all of the cells experienced negative pressures (the cell was too large), ranging from around  −1 ($\beta =0.25$) to  −3 GPa ($\beta =1.00$). The negative pressure is due to the strong bonding between Li and D atoms, which shortens the bond lengths in Li samples with D atoms compared to those in pure liquid Li. The system was then quenched to a target temperature with a trial density, which was chosen by ${{\rho}_{\text{Li}}}(T)+\beta \times \left[{{\rho}_{\text{LiD}}}(T)-{{\rho}_{\text{Li}}}(T)\right]$ at each temperature. A trajectory ranging between 5 and 50 ps was run on the cooled-down system to determine the exact density of this system at each temperature. The length of the trajectory was dependent upon the external pressure of the simulation cell; if a non-zero external pressure was observed at equilibrium, we used the linear extrapolation method described above to obtain a more accurate density, and then the systems with inserted D atoms were cooled down yet again from 1500 K to a new density. The density of each system was updated manually until an accurate density was found to ensure only small external pressures (smaller than 0.2 GPa) existed throughout the entire trajectory. Finally, the correct density at zero pressure at a given temperature of each liquid Li sample with inserted D atoms was obtained; this was then followed by a 20 ps production run to collect statistics.

The statistics collected from a 20 ps production run were used to calculate a list of properties, which provide detailed structural and dynamical information. First, we computed the pair distribution function, defined as

$$\begin{eqnarray}&&g(r)=\frac{1}{\rho N}\langle \underset{i=1}{\overset{N}{\mathop \sum}}\,\underset{j=1,j\ne i}{\overset{N}{\mathop \sum}}\,\delta \left(\mathbf{r}+{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}\right)\rangle,\end{eqnarray} \tag{ 1 }$$

where N is the number of atoms and ρ is the atomic density. ${{\mathbf{R}}_{i}}$ and ${{\mathbf{R}}_{j}}$ are the atomic coordinates of atoms i and j, and $r=|\mathbf{r}|$ is the distance between the two atoms. g(r) describes a probability of finding any two atoms at a given inter-atomic distance, and is ideal for studying systems that do not exhibit long-range order, such as liquids. The Fourier transform of g(r) produces the static structure factor, defined as

$$\begin{eqnarray}&&S(q)=\frac{1}{N}\langle \underset{i=1}{\overset{N}{\mathop \sum}}\,\underset{j=1}{\overset{N}{\mathop \sum}}\,{{\text{e}}^{\text{i}\mathbf{q}\left({{\mathbf{R}}_{\mathbf{i}}}-{{\mathbf{R}}_{\mathbf{j}}}\right)}}\rangle.\end{eqnarray} \tag{ 2 }$$

Both g(r) and S(q) can be determined from x-ray or neutron diffraction experiments, so calculations of both properties can be good benchmarks to test the accuracy of computational models, as was done in our previous work [14]. However, the pair distribution function is used for one element; to analyze systems with two elements, i.e. Li and D, we used the partial pair distribution function [40]:

$$\begin{eqnarray}&&{{g}_{\alpha \beta}}(r)=\frac{N}{\rho {{N}_{\alpha}}{{N}_{\beta}}}\underset{i=1}{\overset{{{N}_{\alpha}}}{\mathop \sum}}\,\underset{j=1}{\overset{{{N}_{\beta}}}{\mathop \sum}}\,\langle \delta \left(\mathbf{r}+{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}\right)\rangle,\end{eqnarray} \tag{ 3 }$$

where species $\alpha \ne \beta$.

While the (partial) pair distribution function provides the relative placements of two atoms, the bond angle distribution function $\sigma (\theta )$ describes the relative positions of three neighboring atoms (or more, depending the choice of radial cutoff). The $\sigma (\theta )$ employed in our study characterizes local structures by giving the distribution of statistically averaged bond angles between each atom and its two nearest neighbors [41]. We used an infinite radius to search for the two nearest neighbors of each atom, recognizing that the two neighbors we find do not necessarily form bonds with the atom. The $\sigma (\theta )$ used in this work is defined such that $\sum_{\theta ={{0}^{{}^\circ}}}^{{{180}^{{}^\circ}}}\sigma (\theta )\Delta \theta =1$ with $\Delta \theta ={{0.05}^{{}^\circ}}$. We further divided $\sigma (\theta )$ into partial bond angle distribution functions, which were calculated by dividing each $\sigma (\theta )$ into different components, with each part representing one specific type of bond angle.

In addition to the aforementioned structural properties, we also evaluated dynamical properties, in particular, the diffusion coefficient, which is derived from the time evolution of the mean square displacement:

$$\begin{eqnarray}&&{{u}^{2}}(t)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathop \sum}}\,|{{\mathbf{R}}_{i}}(t)-{{\mathbf{R}}_{i}}(0){{|}^{2}},\end{eqnarray} \tag{ 4 }$$

where t is time. The diffusion coefficient is obtained from:

$$\begin{eqnarray}&&D=\underset{t\to \infty}{\mathop{\lim}}\,\frac{{{u}^{2}}(t)-{{u}^{2}}(0)}{6t}.\end{eqnarray} \tag{ 5 }$$

It is critical to keep fluctuations in an MD simulation within a normal range so as to maintain equilibrium when evaluating properties. Large fluctuations may lead to spurious phenomena, such as unwanted phase transitions. We therefore consider here the case expected to have the largest fluctuations, namely the heterogeneous liquid LiD simulated with KSDFT-NVT-MD at the highest temperature considered (1143 K). Figure 1 shows the fluctuations of the potential energy, the conserved Nosé–Hoover energy, and the external pressure of the liquid LiD cell in the 20 ps production trajectory. The potential energy is defined as the system's total energy calculated by KSDFT, while the conserved Nosé–Hoover energy is the sum of the potential energy, the kinetic energy of the atoms, and the thermostat's energy. The conserved energy shown in figure 1(a) keeps constant without drifting during the 20 ps equilibrium trajectory, while the potential-energy fluctuations are less than 1.08% of the average value. The measured phase transition pressures of Li and LiD at room temperature are 8 GPa [42] and higher than 94 GPa [43], respectively. Both pressure values are much larger than the observed fluctuating pressure values (around 1 GPa) shown in figure 1(b). We thus conclude that the fluctuations in our MD simulations are not large enough to induce phase transitions in Li and LiD.

Zoom In Zoom Out Reset image size

Solid and liquid Li have been studied extensively with KSDFT [14, 42, 45–47]. For example, KSDFT yielded 455 K as the melting point of Li at zero pressure [42], in excellent agreement with the experimental value of 453 K [48]. Our model's ability to reproduce properties of Li crystals had already been validated [14]. Four different crystal structures of Li were considered: bcc, face-centered-cubic, simple cubic, and diamond. Among them, the bcc structure was found to be the most stable, with a lattice constant of 3.439 Å and a bulk modulus of 14 GPa [14], in good agreement with the respective experimental values of 3.48 Å [49] and 13 GPa [50]. In the present work, we test the same model for liquid Li. Figure 2 shows the calculated densities of liquid Li at a range of temperatures, together with experimental data [44]. We see the density of pure liquid Li (pink solid circles) monotonically decreases as the temperature increases, consistent with the experimental data and corresponding to a simple thermal expansion of a single phase. However, our calculated Li densities are too large by 3.8% when compared to the experimental data at 470 K, with increasing overestimation up to 9.6% at 1143 K. The reason is three-fold: First, the imposed periodic boundary conditions in our simulations produce more ordered structures in liquid Li and consequently higher liquid densities. Second, our single phase simulations do not account for liquid-vapor coexistence present in the experiments. Higher vapor pressures at higher temperatures in the experiments—leading to lower densities—are probably partly responsible for the increased deviation at higher temperatures. Third, impurities in the experimental Li samples that are not present in our pure Li simulations also may affect the measured densities, as noted in [44]. We consider these deviations from experiment to be acceptable because they are much smaller than the density differences between Li and LiD. Additionally, all pure Li samples simulated at temperatures higher than 470 K are in the liquid phase, consistent with the experimental melting point of 453 K [48]. We also simulated pure Li at 400 K (vide infra) in order to compare properties of solid and liquid phases of Li.

Zoom In Zoom Out Reset image size

Figure 3 displays our KSDFT-derived pair distribution functions g(r) and static structure factors S(q) of Li at different temperatures. We expect these data to be quite reliable, given that our previous OFDFT-MD simulations of liquid Li already matched experimental data very well [14]. Here we emphasize the most notable changes in g(r) and S(q), which occur between 400 and 470 K. Specifically, the g(r) of pure Li at 470 K is smooth and does not change considerably as the temperature increases to 1143 K, while the g(r) of pure Li at 400 K has a significantly different shape. The latter exhibits features of solid bcc Li, matching the g(r) of solid bcc Li in [47]. Moreover, we see that the S(q) of pure Li at 400 K in figure 3(c) (inset in figure 3(b)) differs starkly from those S(q) at temperatures higher than 470 K. Specifically, the S(q) of Li at 400 K consists of discretized points whereas the S(q) of Li at 470 K and higher are all relatively smooth curves. The S(q) with discretized points also has peak heights about 20 times larger than the peak heights of S(q) with the smooth curves. We conclude that the S(q) in figure 3(b) correspond to liquid features (no sharp peaks but smooth curves of lower amplitude) whereas the S(q) in figure 3(c) has solid features (isolated high peaks). The four peaks of S(q) in figure 3(c) are located at 2.54, 3.60, 4.43, and 5.08 Å⁻¹, respectively. The location of each peak can be calculated according to the distances of neighbors of an atom in bcc Li. By studying the g(r) and S(q) of Li, we conclude that a solid bcc Li to liquid Li phase transition occurs between 400 and 470 K, which is consistent with both the experimental melting temperature of Li (453 K) and our earlier OFDFT-MD simulations [14].

Zoom In Zoom Out Reset image size

The bond angle distribution functions $\sigma (\theta )$ of pure Li are shown in figure 4(a). All of the $\sigma (\theta )$ of liquid Li exhibit two distinct peaks, with one peak located around ${{64}^{{}^\circ}}$ and the other peak around ${{117}^{{}^\circ}}$. Both peaks are not too different from the bond angles in a bcc structure at 0 K, which are ${{70.5}^{{}^\circ}}$ and ${{109.5}^{{}^\circ}}$, respectively. In general, liquids do not exhibit long-range periodicity but can have organized local structures; the $\sigma (\theta )$ can thus be used as fingerprints to identify pockets of pure liquid Li in a liquid mixture.

Zoom In Zoom Out Reset image size

The calculated diffusion coefficients of liquid Li together with experimental data are displayed in figure 5(b). Although slightly smaller than the measured values (solid line), these calculated diffusion coefficients (black circles) closely follow the experimental trend with increasing temperature [52], with the similar slopes suggesting that the measured activation energy for Li self-diffusion is captured well in our model. Specific predictions of diffusivity are also very accurate. For instance, the calculated diffusion coefficient of liquid Li at 470 K is 0.60 ${{\overset{\circ}{\mathop{\text{A}}}\,}^{2}}$ ps⁻¹ while the measured value is 0.63 ${{\overset{\circ}{\mathop{\text{A}}}\,}^{2}}$ ps⁻¹ [53].

Zoom In Zoom Out Reset image size

In a recent experiment [9], the sputtering yield of Li from thin Li films under high-flux D plasma bombardment was found to be about 5–10 times smaller than expected. To explain this discrepancy, the authors postulated that LiD precipitates might form in liquid Li with high concentrations of D atoms. Accurate KSDFT simulations of the solid and liquid phases of LiD can help to confirm or refute this hypothesis.

We therefore first simulated solid crystal structures of LiD. Three different phases of LiD were examined: the RS, cesium chloride, and zinc-blende crystal structures. LiD in the RS structure is observed at pressures up to 94 GPa at room temperature [43]. KSDFT properly predicts the RS structure to be the most stable (at 0 K and zero pressure), with a lattice constant of 3.998 Å and a bulk modulus of 35 GPa, in very good agreement with the measured values of 4.045 Å and 33.5 GPa [54].

We next simulated LiD at a range of temperatures—initialized in the RS structure—to obtain the theoretical melting point. The temperatures ranged from 470 to 1143 K; we chose 470 K as a lower bound to keep consistent with the liquid Li state we ultimately wish to study. The density of LiD obtained from KSDFT-MD simulations is shown in figure 2 (pink squares). We observe a discontinuity in the density versus temperature data, indicating that a solid–liquid phase transition occurs somewhere between 943 and 1143 K, consistent with the experimental melting point of $965\pm 2$ K [25]. This result suggests that our KSDFT-MD simulations to follow will properly capture the correct solidification of LiD at temperatures 943 K and below. We next discuss partial pair distribution functions, bond angle distribution functions, and diffusivities of LiD at various temperatures.

The three partial pair distribution functions of LiD, namely, ${{g}_{\text{LiLi}}}(r)$, ${{g}_{\text{LiD}}}(r)$, and ${{g}_{\text{DD}}}(r)$ are shown in figure 6. At temperatures lower than 1143 K, all three functions clearly exhibit the features of a solid RS structure. The evidence is twofold. First, ${{g}_{\text{LiLi}}}(r)$ and ${{g}_{\text{DD}}}(r)$ are very similar because the Li and D atoms occupy equivalent coordinates (related by simple translations) in a RS structure. Second, each Li atom is surrounded by D atoms, and vice versa. This is evidenced by the location of the first peak in ${{g}_{\text{LiD}}}(r)$, which represents the shortest distance between any two atoms in a RS structure. When the temperature is increased to 1143 K, dramatic changes simultaneously occur in ${{g}_{\text{LiD}}}(r)$, ${{g}_{\text{LiLi}}}(r)$, and ${{g}_{\text{DD}}}(r)$, suggesting that a solid–liquid phase transition occurs. Specifically, the location of the first peak in ${{g}_{\text{LiD}}}(r)$ shifts from 1.92 to 1.80 Å, while the second (3.52 Å) and third peaks (4.53 Å) disappear. The survival of the first peak means the bonding between Li and D atoms is still strong, while the disappearance of the second and third peaks suggests that the original long-range order in the solid RS structure is destroyed. Instead, a new peak at 4.18 Å indicates that a new local structure forms in the new phase of LiD at 1143 K. At the same time, in both ${{g}_{\text{LiLi}}}(r)$ and ${{g}_{\text{DD}}}(r)$, the magnitude of the first peak drops quickly, with a slight shift in its location from 2.84 to 2.77 Å. As with ${{g}_{\text{LiD}}}(r)$, the second (4.04 Å) and third (5.00 Å) peaks also disappear, and new peaks similarly emerge at 5.38 Å. Based on our observations that all the partial pair distribution functions at 1143 K are more smooth and lower amplitude (liquid-like) and are significantly different than those at 943 K and lower, we again conclude that a solid–liquid phase transition occurs above 943 K, with the melting point of RS LiD lying between 943 and 1143 K, consistent with the experimental value ($965\pm 2$ K [25]).

Zoom In Zoom Out Reset image size

The LiD bond angle distribution functions $\sigma (\theta )$ plotted in figure 4(b) provide further evidence for the melting of LiD at temperatures above 943 K. First, note that each atom in a perfect RS structure has six identical nearest neighbors (as can be seen in the inset of figure 4(b)). By choosing any two atoms from the six neighbors, 90° and 180° are the only two bond angles that can be found, in a ratio of $4:1.$ Similarly, at temperatures lower than 1143 K, the bond angle distribution functions in solid LiD are mainly located at 98° and 172°, with both angles deviating about 8° from those at 0 K due to thermal vibrations and their attendant anharmonicity that increases with temperature. We also observe that the magnitude of the 98° peak is much larger than the 172° peak, qualitatively consistent with the aforementioned $4:1$ ratio. Second, at temperatures of 943 K and below, essentially no bond angles are observed between 120° and 150°. However, when the temperature is increased to 1143 K, we see a sharp decrease of the 98° peak together with disappearance of the 172° peak. At the same time, a broad distribution of bond angles emerges from 120° to 150° decreasing towards 180°, indicating that more disordered structures (liquid-like features) form at this elevated temperature. These observations again support the conclusion that a solid–liquid phase transition occurs above 943 K. Moreover, when comparing the bond angle distribution functions $\sigma (\theta )$ between LiD and pure liquid Li, significant differences are apparent: the $\sigma (\theta )$ of LiD starts around 60° and has peaks at 98° and 172°, while the $\sigma (\theta )$ of pure Li starts around 45° and has peaks at ${{64}^{{}^\circ}}$ and ${{117}^{{}^\circ}}$. We thus conclude that these features in $\sigma (\theta )$ can act as fingerprints of the pure liquid and solid phases of Li and LiD.

The diffusivities of Li and D atoms in both solid and liquid phases of LiD are shown in figure 5. For LiD samples heated up from the RS structure, the diffusion coefficients in the solid phase (at temperatures 943 K and lower) are very small, as expected. Taking Li atoms as an example, the diffusion coefficients are less than 0.01 ${{\overset{\circ}{\mathop{\text{A}}}\,}^{2}}$ ps⁻¹ at both 470 and 943 K, significantly smaller than the diffusion coefficients of Li atoms in liquid Li, e.g. 0.60 ${{\overset{\circ}{\mathop{\text{A}}}\,}^{2}}$ ps⁻¹ at 470 K and 2.38 ${{\overset{\circ}{\mathop{\text{A}}}\,}^{2}}$ ps⁻¹ at 943 K, respectively. The difference is simply that the LiD structures have not melted and Li atoms cannot move freely. A sharp increase in the diffusion coefficient is observed in LiD when the temperature increases from 943 to 1143 K: from less than 0.01 to 1.41 ${{\overset{\circ}{\mathop{\text{A}}}\,}^{2}}$ ps⁻¹ for Li and from less than 0.01 to 1.44 ${{\overset{\circ}{\mathop{\text{A}}}\,}^{2}}$ ps⁻¹ for D. The two orders of magnitudes of increase in the diffusion coefficients provides yet another indicator that a solid-to-liquid phase transition has occurred. Importantly, we predict that the diffusion coefficients of Li and D atoms are strongly correlated, i.e. both types of atoms have almost identical diffusion coefficients at temperatures ranging from 470 to 1143 K. Combining this high correlation in the Li and D diffusivities in LiD with the fact that the first peak in the partial pair distribution function ${{g}_{\text{LiD}}}(r)$ survives the solid–liquid phase transition, we conclude that very strong bonding exists between Li atoms and their first neighbor D atoms in both solid and liquid phases of LiD. Indeed, it is likely that LiD diatomics are diffusing upon LiD melting.

The Li–LiD phase diagram [26] indicates that liquid Li and liquid LiD can coexist in an equilibrium state at 1143 K. However, much remains unknown about the liquid Li samples with inserted D atoms: does liquid LiD form? How do the D atoms affect the local structure of liquid Li? Do D₂ molecules form?

In this section, we first prove the existence of liquid LiD in our simulations and then further analyze its properties. The analysis is mainly based on the bond angle distribution functions $\sigma (\theta )$. Figure 4(c) shows the calculated $\sigma (\theta )$ for pure liquid Li, pure liquid LiD, and liquid Li samples with four different concentrations of D atoms. Several findings are described in what follows.

The $\beta =1.00$ liquid enters the equilibrium state by being impregnated with randomly distributed D atoms and then being cooled down from 1500 K to a target temperature. By contrast, the pure LiD sample is heated up from a solid RS structure at 0 K to a target temperature. The physical processes used to prepare the two samples serve different purposes: the former is used to study what will happen after inserting D atoms into liquid Li without any memory of solid LiD while the latter is used to gauge the melting point of solid RS LiD. Though the two systems were generated differently, we find that the $\sigma (\theta )$ of the $\beta =1.00$ liquid (cyan triangles in figure 4(c)) is almost identical to the $\sigma (\theta )$ of the LiD liquid (pink squares in figure 4(c)). In particular, the major peaks in both $\sigma (\theta )$ are located at ${{98}^{{}^\circ}}$. As previously discussed, we know the ${{98}^{{}^\circ}}$ peak in $\sigma (\theta )$ is a general feature of liquid LiD at zero pressure. We thus suspect that the $\beta =1.00$ system is mostly if not all liquid LiD, based on the observation that the $\sigma (\theta )$ of the $\beta =1.00$ sample is nearly identical to the $\sigma (\theta )$ of the LiD sample. To further confirm this hypothesis, we calculated the partial bond angle distribution functions: the $\sigma (\theta )$ was decomposed into six partial functions, with each part representing one specific type of bond angle. We selected and summed up the two bond angle types existing in a RS LiD structure, i.e. bond angles formed by connecting a Li atom with two D atoms, and vice versa. The results are shown in figure 7(b), from which we see that the ${{98}^{{}^\circ}}$ peak (a signature of liquid LiD) in $\sigma (\theta )$ increases with β, and reaches the highest magnitude in the $\beta =1.00$ sample. Based on these findings, we again conclude that the amount of liquid LiD formed at 1143 K increases with increasing concentrations of D atoms in liquid Li, and that primarily LiD liquid forms in the $\beta =1.00$ liquid. This finding suggests that full retention of D atoms in liquid Li is possible in mixtures of liquid Li and liquid LiD, at least up to a $1:1$ ratio of D to Li.

Zoom In Zoom Out Reset image size

The partial bond angle distribution functions in figure 7(a) shows the bond angle type formed by three Li atoms, obviously the only one possible in pure liquid Li. However, this type of bond angle is strongly suppressed with increasing β and hardly exists at all in the $\beta =0.50$ and higher samples. The rapidly decreasing magnitude of this bond angle distribution also implies that the pure Li local structures change a great deal with increasing β.

Figure 7(d) shows a type of bond angle formed by one D atom connecting to one Li atom and another D atom. This bond angle could occur either via formation of D₂ molecules solvated by Li atoms or via LiD molecules interacting with a solvated D atom. Note that the range of vertical axis in figure 7(d) is two orders of magnitude smaller than figures 7(a)–(c), suggesting the type of bond angle shown in figure 7(d) is much less probable than the bond angle types in figures 7(a)–(c). Thus we expect D₂ formation to be a rare event, though not impossible. Not a single bond angle ever was formed by three D atoms in any trajectory at any temperature considered. Thus, D clusters of three or more D atoms will not exist at these concentrations. This observation also implies that D atoms are distributed quite uniformly among Li atoms regardless of the D concentration in liquid Li.

Nuclei of solid LiD are expected to form in the liquid Li host (due to the liquid–solid phase transition) when liquid Li samples with inserted D atoms are cooled down below the melting point of RS LiD ($965\pm 2$ K). We will study these samples by analyzing their densities, pair distribution functions, bond angle distribution functions, diffusivities, geometries, and electron distributions. By analyzing all of these data, we are able to confirm that the liquid–solid transition occurs and RS LiD forms. We also observe some D₂ molecules in the $\beta =1.00$ sample. Details are provided in the following subsections.

3.5.1. Liquid to solid phase transition.

In general, samples at lower temperatures and higher D concentrations have the following four noteworthy features that furnish evidence of liquid–solid phase transformations.

First, the discontinuities observed in the density versus temperature plots for LiD and liquid Li samples with high concentrations of inserted D atoms are similar (figure 2). Because the discontinuity found in LiD's density versus temperature behavior was found to correspond to a solid–liquid phase transition, we suspect the latter samples also undergo this transition. Phase transition temperatures can be very roughly estimated from the midpoint of the temperature range where the discontinuities occur in the density curves during cool-down. We thus estimate them to be around 500, 750, and 900 K in the $\beta =0.50$, 0.75, and 1.00 samples, respectively. The $\beta =0.25$ sample may not have undergone a phase transition because the discontinuity in its density versus temperature curve is not significant compared to the others.

Second, we select temperatures near the discontinuities in the density curves and plot the partial pair distribution functions (${{g}_{\text{LiLi}}}$, ${{g}_{\text{LiD}}}$, and ${{g}_{\text{DD}}}$) in figure 8. Two temperatures are selected for each sample with a different value of β, i.e. one above and one below the discontinuity. We find that the ${{g}_{\text{LiLi}}}$ from $\beta =0.50$, 0.75, and 1.00 have solid-like features (more than two peaks) at selected lower temperatures and liquid-like features (smooth curves with only two peaks with lower amplitudes) at selected higher temperatures. This analysis provides more evidence of the liquid–solid phase transition. Note also that the ${{g}_{\text{LiLi}}}$ from $\beta =0.50$ and 0.75 in figure 8 closely resemble the ${{g}_{\text{LiLi}}}$ in both pure liquid Li (see figure 3(a)) and in solid LiD (see figure 6(a)), suggesting that at these intermediate concentrations we have evidence of nucleation of LiD in liquid Li.

Zoom In Zoom Out Reset image size

Third, the bond angle distribution functions $\sigma (\theta )$ for all samples with different β are illustrated in figure 9 for temperatures ranging from 470 to 943 K. Large changes in the $\sigma (\theta )$ are observed in all samples upon cooling. Specifically, the changes are observed at three locations in the $\sigma (\theta )$: (i) the height of the peak at 98° increases significantly; (ii) the peak at around 172° emerges; and (iii) contributions to $\sigma (\theta )$ between 120° and 150° disappear. In general, more localized and higher peaks in $\sigma (\theta )$ suggest more ordered local structures compared to those $\sigma (\theta )$ with relatively lower and broader peaks, suggesting that solid structures indeed formed in these samples as we cooled them.

Zoom In Zoom Out Reset image size

Finally, discontinuities in the diffusion coefficients of both Li and D atoms are observed (figure 5). The temperatures where these discontinuities occur are exactly the same as those temperatures found for the density curve discontinuities (figure 2). Notably, the diffusion coefficients have a much clearer boundary between solid and liquid phases than our earlier metrics, with 1–2 orders of magnitude smaller diffusion coefficients for the lower temperature phase compared to those of the higher temperature phase.

Overall, when liquid Li samples with D atoms are quenched down from 1500 K, liquid–solid phase transitions occurred at higher temperatures for higher concentrations of D. This trend is to be expected if more nuclei of LiD are forming, since LiD has a higher melting point than Li.

3.5.2. Formation of rock-salt LiD.

To reveal details of the LiD structures formed as a result of the liquid–solid phase transitions, geometries from the end of the 20 ps production runs are shown in figure 10. We plot two three-dimensional snapshots at temperatures above and below the liquid–solid phase transitions for samples at each D concentration. Figures 10(a)–(c) reveal that a subset of the atoms form nuclei of solid LiD because the number of Li atoms is larger than the number of D atoms in all three of these cases. Meanwhile, the disordered structure formed by the rest of the Li atoms is likely to be liquid Li, as evident in figure 10(a), though harder to infer from figures 10(b) and (c). Figure 10(d) shows the $\beta =1.00$ structure at 793 K, which mostly exhibits a typical RS structure where most Li and D atoms are bonded to the other type of atom, and the bond angles comprise one Li and two D atoms or vice versa. By contrast, figures 10(e)–(h) show the hallmarks of liquids, namely fully disordered structures, as we now describe in more detail.

Zoom In Zoom Out Reset image size

The partial pair distribution functions ${{g}_{\text{LiD}}}$ (in figure 8) in samples with different β keep the same shapes for all solid samples: three peaks in ${{g}_{\text{LiD}}}$ located around 1.91, 3.49, and 4.50 Å, respectively. Moreover, these peak locations are very close to those of RS LiD (1.92, 3.51, and 4.53 Å), as shown in figure 6. We conclude that both ${{g}_{\text{LiD}}}$ in RS LiD and the $\beta =1.00$ sample have the same local structures. Moreover, the bond angle distribution functions $\sigma (\theta )$ in figure 9 indicate that when the liquid–solid phase transition occurs, most atoms with bond angles between 120° and 150° are changing to ${{98}^{{}^\circ}}$ or 172° (typical bond angles of RS LiD at nonzero temperatures) and the resulting $\sigma (\theta )$ is very close in shape to the $\sigma (\theta )$ of RS LiD (figure 4(b)). All these findings suggest that RS LiD nuclei can form after inserting high concentrations of D atoms into liquid Li samples, as long as the temperature of the liquid Li is below the melting point of LiD.

3.5.3. Diffusivities.

3.5.4. Existence of D₂ molecules.

Evidence for the formation of D₂ molecules is discussed here. First, the partial pair distribution function ${{g}_{\text{DD}}}$ of the $\beta =1.00$ sample (figure 8(l)) has a small peak between 0.69 and 0.90 Å which is close to the bond length of a D₂ molecule (0.74 Å). We suspect that this peak represents D₂ molecules vibrating in the $\beta =1.00$ sample at 793 K. Because the charge state of a D atom is very different in a D₂ molecule and LiD, we evaluated charge states of all D atoms in all samples to verify the existence of D₂ molecules. This analysis provides us a second piece of evidence for D₂ formation. Specifically, we used Bader charge analysis [55] of the electron distributions, which decomposes the total electron density of a system into atomic contributions. Figure 11 displays the valence electron distributions of Li and D atoms in the $\beta =1.00$ sample at two temperatures. It is immediately evident that the majority of Li and D are present as ions in this sample, as roughly Li⁺ and D⁻, judging from large peaks at small (Li) and large (D) numbers of valence electrons. As a reference, the electron numbers for Li and D atoms in RS LiD at 0 K were calculated to be 0.17 for Li and 1.83 for D; those values appear at the center of the large peaks mentioned above. Thus the bonding in the $\beta =1.00$ sample is ionic, very similar to that of RS LiD. However, we do see that some very small peaks exist near 1.0 electron, corresponding to neutral D, which would be present only in D₂. Note that we find D₂ molecules only in the $\beta =1.00$ samples. Bader charge analyses of samples containing lower concentrations of D reveal peaks with broader distributions. For example, Li atoms with Bader charges between 0.15 and 1.00 can be found, suggesting both ionic and metallic Li are present in mixtures of LiD and Li metal. By contrast, the D atoms in these lower concentration samples all have Bader charges larger than 1.75, corresponding to D anions in the LiD regions, as in the majority of the $\beta =1.00$ material, but with no evidence of formation of D₂ molecules.

Zoom In Zoom Out Reset image size

Interestingly, we find that these D₂ molecules are not stable at 1143 K; they survive in the first few ps simulations but then dissociate when the system approaches equilibrium at this higher temperature. Similar conclusions were reached in previous work [56] examining liquid LiD at temperatures of 1015 K and higher.

As discussed earlier, the diffusion coefficients in the $\beta =1.00$ sample are consistently larger than those found in the pure LiD sample, despite the concentrations of D atoms being the same. The existence of D₂ molecules in the $\beta =1.00$ sample is responsible for the different diffusivities. The compositions in the two samples are different: the $\beta =1.00$ sample is a mixture of mostly LiD with a few D₂ molecules while the pure LiD sample only has LiD. The D₂ molecules contribute some extra diffusivity to the Li and D atoms, resulting in larger diffusion coefficients for both type of atoms. Lastly, the D₂ molecules are responsible for the larger volume of the $\beta =1.00$ sample when compared to the volume of the LiD sample (as is clear from the density curves in figure 2).

3.5.5. A new local structure.

3.5.6. Metastable states.