Cu Diffusion in Amorphous Ta₂O₅ Studied with a Simplified Neural Network Potential

To construct the simplified NN potential, 2000 training data were obtained via DFT. The resulting data was randomly split into a training set (90%) and an independent testing set (10%). The training set was used to fit the parameters of the NN, which consisted of two hidden layers with 9 nodes each, one input layer with 36 nodes, and one output node. The NN parameters (corresponding to the number of hidden layers and number of nodes in each layer) were determined by performing a special testing procedure.⁵⁸⁾

The optimization of NN potential parameters was done iteratively with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. In each iteration during the optimization, the root-mean-square errors (RMSEs) in both the training set and testing set were monitored. The training process was truncated when testing RMSE stops decreasing. The learning curve of such NN potential was shown in Supplemental Material (Fig. S2).⁵⁹⁾ The RMSEs for energies after training process were 23 and 39 meV per structure for the reference training and independent testing set, respectively. The mean absolute errors (MAEs) obtained for the training and test structures were 16 and 21 meV/structure, respectively. According to the results presented in Fig. 1, the energies of the training and test structures calculated using the NN potential matched the values obtained via DFT simulations (the corresponding correlation coefficients were equal to 0.998 and 0.997, respectively). In contrast with the high-dimensional NN potentials, whose remaining total energy errors are typically under 5 meV/atom, the accuracy of simplified NN potential appears to be worse. However, it should be pointed out that the RMSEs in the present study are given as values per whole supercell, while they are given as the values per atom in the high-dimensional NN potentials. So the comparison of the potential accuracy between the simplified and high-dimensional NN potentials is not straightforward. On the other hand, the fitting accuracy of NN can be constantly improved by increasing training data at the expense of additional DFT calculations. Considering the fact that the satisfactory accuracy was obtained for diffusion barriers (which will be shown in Sect. 3.2), we decided to perform the following investigation based on the present NN potential.

»View larger version

Figure 1. (Color online) A comparison between the total energies of the training set and testing set structures obtained using the DFT and NN approaches. The inset contains the MAE distributions (\(|\mathrm{E}_{\text{NNP}}-\mathrm{E}_{\text{DFT}}|\)) calculated for these two data sets.

It should be noted that the calculation of the total energy of a single geometric configuration, which was conducted using the current implementation of the NN potential model on a single CPU core (Intel^® Xeon^® W3565 processor), took only 0.001 s (for comparison, performing the same task using the VASP code and 12 CPU cores of the same type took about 1 h). Hence, the computational speed obtained using the NN potential was \(10^{7}\) times greater than that observed during DFT simulations.

To test the reliability of the NN potential for describing the atomic diffusion process, the former was used to predict a long path of the Cu diffusion in the amorphous Ta₂O₅ matrix and the corresponding energy profile via NEB simulations. During this procedure, the overall motion of the Cu atom was described by several hops between the equilibrium sites. The locations of images after optimization (representing the diffusion path) are shown in Fig. 2(a). The obtained NN potential was in good agreement with the DFT results. The lengths of the diffusion paths defined as the sums of distances between the adjacent Cu transient positions were equal to 28.6 and 28.9 Å for the NN potential and DFT simulations, respectively. The corresponding energy variations along the long-range diffusion path, which are plotted against the diffusion length in Fig. 2(b), also demonstrate a good agreement between the two utilized methods. The largest mismatch between the data obtained via the NN potential and DFT calculations was about 0.15 eV, which corresponded to the saddle point of the highest barrier. Apart from the inherent inaccuracy of the NN potential method, a large discrepancy observed for this point might also result from the NEB images that missed the saddle point of the minimum energy path, which represented a common problem for NEB calculations and could be fixed by using a climbing image NEB method.⁶⁰⁾ It is worth noting that the range of Cu–O coordination number is wider in the training structures than the transient images and the range of Cu–Ta coordination number is the same between them (for details, see Sect. B of Supplemental Material).⁶¹⁾ This suggests that the training dataset includes a variety of local environments sufficient to describe the transient image structures, which may be one of the reasons why the simplified NN potential method can achieve a good accuracy in the NEB calculation.

»View larger version»View larger version

Figure 2. (Color online) A comparison between the NEB calculation results obtained for the selected diffusion pathway of a Cu atom in amorphous bulk Ta₂O₅ using the DFT and NN potential approaches. (a) Transient images of the Cu atom obtained during the NEB calculations, from which the locations of diffusion pathways can be acquired. The convex polyhedrons represent Ta–O species in the amorphous matrix, while O and Ta atoms are not shown for clarity. (b) Energy profiles computed along the corresponding diffusion paths.

The advantages of the simplified NN potential are as follows: 1) it is based on the multi-layer perceptron NN, which can be easily constructed using various open-source codes; and 2) the number of data sets required for training the NN is about 1–2 orders of magnitude smaller than that utilized for a typical ternary high-dimensional NN potential.

However, the proposed simplified method also has limitations. First, it cannot be extended to MD techniques because of the absence of atomic descriptions (especially forces acting on atoms) in the host matrix. Second, although the NN potential can be trained in a particular amorphous network, it is often desirable to transfer it to another amorphous matrix. However, it was found that the transferability of the simplified potential was not very high (except for the cases with a very similar matrix structure) concerning the energy of the system including a Cu atom inserted at a randomly chosen point.⁶²⁾ For the calculation of metastable Cu positions, migration paths and migration energy barriers, the transferability seems much better, but in some cases the errors are considerably large.⁶²⁾ We expect that the transferability of the proposed method can be significantly improved by training the NN potential with a wide variety of the local environments of Cu atoms (which would require the use of various amorphous host matrices) and/or larger cutoff radii of the symmetry functions (this hypothesis will be tested in future studies).

On the basis of the above, we think that the NN potential method proposed in this paper can be complementary to the high-dimensional NN potential in studying atomic diffusion. Highly reliable description of atomic diffusion behavior is possible by a well-constructed high-dimensional NN potential, while at least rough overview can be obtained by the present NN potential with much lower computational cost.

Determining potential equilibrium sites in amorphous bulk materials is a challenging task because of structural disorder. To avoid bias, it is important to locate all interstitial sites without making assumptions about their structure. Owing to the high computation speed of the NN potential, all possible Cu positions in the amorphous Ta₂O₅ matrix structure were investigated in this work. First, Cu atoms were inserted at each of the points defining the \(50{\times} 50{\times} 50\) uniform grid in the matrix supercell, and their positions were optimized using the BFGS algorithm. After optimization, many of the final Cu positions were close both spatially and energetically. This phenomenon represented an artifact originated from the lax force convergence criterion, which terminated Cu optimization slightly before reaching energy minimum. Therefore, hierarchical clustering analysis was performed for the optimized interstitial positions using the distances between two Cu atoms as the metric⁶³⁾ (it should be noted that the described technique was used in previous studies to search for equilibrium sites in amorphous bulk materials^64,65)). As a result, the calculated interstitial positions were successively merged until all the distances between clusters became longer than the cutoff distance of 0.5 Å (during this procedure, 29 different clusters were obtained). For each cluster, the energy variance of all configurations was below 0.04 eV, and the lowest energy configuration was used to represent the metastable Cu interstitial position.

According to the results obtained in the previous section, 29 metastable interstitial positions of Cu atoms were determined using the NN potential. Assuming that a Cu atom can hop between the two adjacent interstitial positions with a distance between them less than 5 Å, 67 diffusion paths have been obtained (during the search, the paths that can be considered combinations of two or more individual pathways were eliminated). After that, the activation energies of the resulting pathways were determined using the NEB method (the spatial locations of all the obtained diffusion paths are depicted in Fig. 3).

»View larger version

Figure 3. (Color online) Equilibrium sites (blue spheres) and diffusion pathways (blue lines) determined for the Cu atoms in the bulk amorphous Ta₂O₅ structure. The polyhedrons represent Ta–O polyhedrons, while O and Ta atoms have been removed for clarity.

For each individual diffusion path, two different activation energies corresponding to the two opposite hopping directions were computed. In general, for the path connecting sites i and j with the saddle point energy \(E_{ij}\) (which represents the maximum energy of the NEB energy profile), the activation energy of Cu diffusion from i to j is equal to \(E_{ij}\)–\(E_{i}\) (here \(E_{i}\) is the energy of site i), while the activation energy of the diffusion in the opposite direction equals \(E_{ij}\)–\(E_{j}\). The calculated barrier energies of the diffusion paths span across a wide range between 0.01 and 1.79 eV with an average value equal to 0.55 eV, while their lengths vary from 0.91 to 4.98 Å. Figure 4 displays the distribution of barrier energies determined for single hops. Although most of the pathways exhibit relatively small barrier energies of less than 0.2 eV, it does not imply that the corresponding Cu diffusion activation energies are very low because a particular Cu atom may not be able to escape from the narrow region between the sites connected via the paths with such small energies. For this reason, the periodical diffusion paths (which are defined as the pathways connecting equilibrium sites and their equivalents in the adjacent supercells) have been sampled. As a result, the shortest pathways (determined from the obtained number of nearest-neighbor hops) and the second shortest ones were selected. The barrier energy of a periodical path was equivalent to the largest barrier energy calculated for individual hops. The results of the statistical analysis summarized in Fig. 4 reveal that the obtained barrier energies of the periodical paths vary from 0.68 to 1.79 eV. Interestingly, the lowest value (0.68 eV) corresponds to the activation energy of Cu diffusion in the diffusion network obtained via kinetic Monte Carlo (KMC) simulations, which will be described in detail in the next subsection.

»View larger version

Figure 4. (Color online) Distributions of the barrier energies calculated for all diffusion pathways of atomic Cu in bulk amorphous Ta₂O₅. The black columns correspond to individual hops (single NEB calculations), and the red columns represent the periodical paths, which connect equilibrium sites with their equivalent sites in the adjacent supercells. The frequencies were normalized with respect to the total number of diffusion pathways.

The diffusivity and diffusion activation energy of Cu atoms were examined by performing KMC simulations. During the utilized procedure, the probability for the Cu atom hopping from interstitial site i to site j was determined using the harmonic transition state theory: \begin{equation} \mathrm{p}=k\exp\left(-\frac{\Delta E_{ij}}{k_{\text{B}}T}\right), \end{equation} (6) where the prefactor k represents the jump frequency (equal to about \(10^{13}\) s⁻¹), and \(\Delta E_{ij}\) is the hopping barrier between sites i and j. At the beginning of the simulations, Cu atoms were distributed randomly among the interstitial sites according to the Boltzmann-weighted energy relationship. Afterwards, the initial state was equilibrated by running 5000 MC events per atom. After reaching equilibrium, additional 10 million MC events were run for 1000 atoms, and the self-diffusivity of Cu atoms was evaluated from their respective trajectories using the Einstein expression: \begin{equation} \mathrm{D}=\lim_{t\to\infty}\frac{1}{6t}\langle|\boldsymbol{{r}}(t)-\boldsymbol{{r}}(0)|^{2}\rangle. \end{equation} (7) Figure 5(a) describes the self-diffusivity of Cu atoms in a-Ta₂O₅ estimated at temperatures ranging from 500 to 1400 K. After fitting the obtained temperature dependence of the Cu self-diffusivity with the Arrhenius expression \(\mathrm{D}=\mathrm{D}_{0}\exp(-E_{\text{eff}}/k_{\text{B}}T)\), where D₀ and \(E_{\text{eff}}\) were the diffusion prefactor and effective activation energy, respectively, slight differences between the activation energies calculated for the low-temperature (500–800 K) and high-temperature (900–1500 K) regions were observed. The activation energy at low temperature is 0.67 eV, which is 0.2 eV lower than the high temperature activation energy 0.87 eV. In addition, the obtained low-temperature activation energy values matched the lowest barrier energies of the periodical paths obtained in the previous subsection (4.2), suggesting that at low temperatures, one or several periodical paths with the lowest barrier energies determined the hopping behavior of Cu atoms, while the contribution of the pathways with higher barriers [which were calculated via Eq. (6)] increased rapidly with increasing temperature.

»View larger version»View larger version

Figure 5. (Color online) (a) Cu self-diffusivity in the amorphous bulk Ta₂O₅ structure calculated from the results of KMC simulations. The two black lines represent the linear fits of the simulation data obtained at high temperatures (900–1400 K) and low temperatures (500–800 K). (b) The comparison of Cu diffusivity between experimental observations and KMC simulation. The data of experiments 1 and 2 are taken from Refs. 6 and 7, respectively.

Lastly, we compared the diffusion behaviors extracted from KMC simulation with the experimental observations. Two independent experiments about the Cu diffusion in amorphous Ta₂O₅ thin film have been reported.^6,7) The diffusivities evaluated from the experiments are shown in Fig. 5(b), as well as the KMC simulation results in the same temperature range. In the first experiment, the Cu diffusivity at 200–500 °C was measured with secondary ion mass spectroscopy.⁶⁾ By fitting the measured data to the Arrhenius relation, the activation energy for Cu diffusion is estimated to be 0.64 eV, while the room temperature diffusivity is estimated to be \(4.9\times 10^{-21}\) cm²/s. In the other experiment, the Cu diffusivity in amorphous Ta₂O₅ film deposited by electron-beam was evaluated to be \(10^{-13}\) cm²/s from cyclic voltammogram measurements.⁷⁾ The discrepancy of about eight orders of magnitude in the room temperature diffusivity may be ascribed to the difference in the sample thickness and the structure of the deposited Ta₂O₅ film.⁷⁾ Though the large discrepancy between the two experiments makes the quantitative comparison with our simulation results difficult, we may be able to say that our KMC results are reasonable. Firstly, the low temperature activation energy obtained from KMC simulation, 0.67 eV, is close to the one estimated from the first experiment, 0.64 eV. Secondly, the room temperature diffusivity extrapolated from the KMC results, \(5.7\times 10^{-17}\) cm²/s, is in between the two experimental estimates.