Using Metadynamics to Obtain the Free Energy Landscape for Cation Diffusion in Functional Ceramics: Dopant Distribution Control in Rare Earth‐Doped BaTiO3

Barium titanate is the dielectric material of choice in most multilayer ceramic capacitors (MLCCs) and thus in the production of ≈3 trillion devices every year, with an estimated global market of ≈$8330 million per year. Rare earth dopants are regularly used to reduce leakage currents and improve the MLCC lifetime. Simulations are used to investigate the ability of yttrium, dysprosium, and gadolinium to reduce leakage currents by trapping mobile oxygen defects. All the rare earths investigated trap oxygen vacancies, however, dopant pairs are more effective traps than isolated dopants. The number of trapping sites increases with the ion size of the dopant, suggesting that gadolinium should be more effective than dysprosium, which contradicts experimental data. Additional simulations on diffusion of rare earths through the lattice during sintering show that dysprosium diffuses significantly faster than the other rare earths considered. As a consequence, its greater ability to reduce oxygen migration is a combination of thermodynamics (a strong ability to trap oxygen vacancies) and kinetics (sufficient distribution of the rare earth in the lattice to intercept the migrating defects).


Introduction
Functional materials are ubiquitous in modern technology and close control of their defect chemistry is central to the performance of many electronic, magnetic and optical devices. Simulation is a powerful tool for providing detailed information on the equilibrium geometry and energetics of the defects. Their distribution often controls their effectiveness within the material. Traditional simulation, however struggles with modeling some kinetics, particularly diffusion in solid state systems where and polycrystals [21] giving activation energies ranging from 0.5 to 1.28 eV. [20][21][22] Static and dynamic simulation studies using both density functional theory [23] and classical methods [24][25][26] have been used to look at the migration pathway, thus providing insight into diffusion in the undoped material.
Previous simulations of RE 3+ doping in BaTiO 3 [27] have shown that large RE 3+ ions (La, Ce, Pr, Nd, (Pm), Sm) substitute mainly on the A-site but small RE 3+ ions (Tm, Yb, Lu, Sc) prefer the B-site. Mid-sized RE 3+ ions (Eu, Gd, Tb, Dy, Y, Ho, Er) dope on both sites giving self-compensation. These results are in good agreement with experiment. [28] The influence of both divalent and trivalent dopant ions on oxygen migration has been studied in SrTiO 3 showing that all the trivalent ions strongly trap oxygen vacancies when doped solely on the B-site (due to the effective negative charge of the defects) and that the influence of the RE 3+ dopants extends further than two lattice sites into the surrounding crystal. [29] Quantum mechanical simulations [30] suggest that B-site dopants are most effective at trapping oxygen vacancies in BaTiO 3 but only the closest lattice sites to individual dopants have been considered due to the large computational demands of these simulations.
The results of these studies contradict much of the accepted wisdom associated with MLCC production within industry where it is well known that only certain mid-range RE 3+ cations, Ho, Dy, and Y, extend capacitor lifetime significantly. [11,14,28] Dy is generally considered to be the best, while the larger RE 3+ cations from La to Gd have little or no effect. This suggests the thermodynamics of oxygen vacancy trapping alone are not sufficient to explain how REs extend the lifetime of capacitors. We need to understand the distribution of these RE 3+ cations within the lattice and the kinetics of how they reach this distribution. Cation diffusion in BaTiO 3 is too slow to be simulated using standard methods. We therefore employ a bias potential approach coupled to a statistical mechanics population analysis to understand the distribution of these RE 3+ cations in BaTiO 3 , examine how this will affect the diffusion of oxygen vacancies and ultimately explain the enhancement in lifetime of MLCCs. This study illustrates the strength of using a range of simulation techniques in combination with experimental data to solve a puzzle of considerable industrial interest. It also demonstrates the power of bias-potential methods such as metadynamics to investigate slow processes like cation diffusion in many functional oxides which are difficult to address by simple dynamical simulations or indeed by experiment.

Defect Chemistry and Defect Clustering
Five defect compensation schemes (four published by Buscaglia et al. [31] and a fifth by Freeman et al. [27] ) were simulated for Y 3+ , Dy 3+ , and Gd 3+ dopants. The schemes are (Schemes 1-5): The dopant solution energies in the dilute limit, E s 0 (i), corresponding to the dissociated defect energies for each compensation scheme, i, were calculated along with the binding energies, E b (i), (the energy gained by associating the defects in the scheme adjacent to each other) by simulating the defects together. The binding energies were then used to calculate the final solution energy, E s (i). These are all given in Tables A3.1-A3.3 and Appendix S3 in the Supporting Information. The probability that a given RE ion in the lattice will exhibit a given compensation mechanism, i, at equilibrium over a range of temperatures, P RE (i) was calculated using statistical mechanics where E s (i) is the solution energy discussed above; k B is the Boltzmann constant and T is the temperature in Kelvin. The results for the associated and dissociated defect compensation scheme energies across a temperature range covering both working and sintering temperatures are shown in Figure 1a-c. This shows that Y 3+ , Dy 3+ , and Gd 3+ dope mainly by selfcompensation in agreement with experimental work of Lu and Cui [32] for Dy 3+ . Studies of the solubility of Ho 3+ , Dy 3+ , Y 3+ in BaTiO 3 also show [33] that these ions can dope onto either A or B sites; with site preference depending on the solution composition. All the dopants investigated here show associated and dissociated self-compensation RE Ba • and RE Ti ′ defects as well as RE Ti ′ defects compensated for by V O •• at high temperatures. As the ionic radii increase from Y 3+ to Gd 3+ , the probability of finding the RE Ti ′ defect decreases. Given the data in Figure 1, we would therefore expect the majority of the RE cations to adopt the selfcompensation scheme with a small percentage of B-site doping.

Dynamics of Oxygen Vacancy Migration
Nominally "undoped" BaTiO 3 crystals are weakly acceptor doped by impurities which are compensated for by oxygen vacancies in the bulk. [6] The enthalpies of migration for these vacancies have been reported as 0.68 ± 0.06 eV (0.70 ± 0.04 eV) from experiment [20] while simulations (discussed in ref. [20]) have given defect migration energies of ≈0.77 eV (lattice statics with classical potentials [34] ) and ≈0.89 eV (lattice statics with density functional theory (DFT) [23]  simulations using absolute position collective variables (for technical details see Section 2 and Appendix S2 in the Supporting Information) produce an oxide vacancy migration energy of 0.90 eV. This is a free energy of migration at the simulation temperature of 298 K (metadynamics) and not directly comparable with the measured enthalpy but the difference is small at low temperatures. [35] The free energy landscape for the oxide vacancy diffusion can be seen in Figure 2. By using the metadynamics methodology we are able to extract more than the barrier alone and can observe the pathway the diffusing ion takes without significantly biasing the choice of the pathway as many other techniques would do.
The positive oxygen vacancy will be trapped by the negative RE Ti ′ defect in the defect schemes identified by our simulations. For the less favorable 2RE Ti ′ + V O •• (Scheme 3) the binding energy of an additional oxygen vacancy is ≈−2.3 eV for all the RE defects. This represents weaker binding than the ≈−3.4 eV observed for the self-compensation schemes (Scheme 4). This demonstrates why the mid-sized RE cations (which predominantly dope with the self-compensation scheme) are most successful at trapping oxygen vacancies and extending the lifetime of BaTiO 3 based devices.
The effect of A and B site self-compensating rare earth pairs on the migration of oxygen vacancies in the lattice was further investigated using lattice statics methods. The 60 nearest sites for oxygen vacancies near RE dopants were simulated for two situations: an oxygen vacancy around a RE-doped self-compensating pair and an oxygen vacancy around a lone RE doped on a B site. Full details of the results can be seen in Appendix S5 in the Supporting Information. To summarize, the results show that both defect systems can trap oxide vacancies approaching from any direction. This demonstrates that both associated and disassociated self-compensation schemes can trap the oxide vacancies.
We use the potential energy surface of Figure SI5.1 in the Supporting Information to map out the lowest energy pathway for the diffusion of the oxide vacancy away from the RE defect. Each lattice site is labeled with a number (1 is closest to the RE and 6 is the furthest away). The results are shown in Figure 3.

For the RE Ti
′ case (Figure 3b), sites 1 and 2 are equally favorable, with a small barrier to migration of 0.2 eV between the sites. There is a much greater barrier to moving the vacancy further away, from site 2 to 3, so the oxide vacancy and RE are no longer neighbors. This is expected, due to the large coulombic attraction between the negative RE Ti ′ and the positive V O •• making it difficult for the vacancy to escape from the dopant. Conversely going from site 3 to 2, i.e., the oxide vacancy moving toward the RE defect, the barrier to migration is small or nonexistent.  Contour plot of the free energy landscape for oxide vacancy migration in BaTiO 3 . View from side on (left) and on top of surface (right). Energy profile indicated by colors with blue-purple lower in energy and yellow-orange higher in energy. Energies recorded in eV. The migration pathway follows the oxygen ion moving from the right-hand side well to the left-hand side well. Despite the wells representing identical sites the wells are not identical due to the diffusion process and metadynamics method. Therefore as the metadynamics simulation proceeds the surface evolves and regions of the free energy surface can be adjusted even with little contact from the moving ion. This means that the two minima will not necessarily match as the minima beyond the barrier evolves prior to the atom fully exploring this part of the surface.
In the RE associated pair case (Figure 3a) a similar behavior is seen with a small preference for site 1 compared to site 2. The small barriers present in moving the oxide vacancy toward the RE defect (i.e., from 6 to 1) imply that an oxide vacancy moving within four lattice sites of the RE defect will be pulled toward it and trapped. Assuming perfect mixing and a uniform distribution, this suggests that 2-3 at% of RE dopants could ensure that no oxide defect could pass through a barium titanate ceramic without moving within four lattice sites of a RE defect (see Appendix S7 in the Supporting Information for details of the argument). It can be seen as we move to the furthest distance between the RE and oxide vacancy (i.e., site 6) then the barrier to diffusion has largely returned to bulk values of ≈0.8 eV (which is the value we calculated for isolated oxygen vacancy diffusion in the RFO mechanism) indicating the final reaches of the influence of the RE dopants.
All the REs investigated (Figure 3a,b) have similar energetics for trapping oxygen vacancies in both doping scenarios despite the differences in ionic radii. This implies the ability of the RE to trap oxide vacancies should be primarily linked to the probability of the RE adopting the self-compensation scheme as the more REs present in this scheme the greater the number of potential trapping sites. Following this model, we would expect to see better trapping with increasing size of the mid-range RE defects which increasingly favors self-compensation and therefore Gd 3+ should be the most efficient. This is at odds with the experiments which all suggest that Dy 3+ is the most efficient. This raises the question of what other factor explains this discrepancy? If the different RE ions have a different distribution in the lattice then this could limit the ability of the RE cations to trap oxide vacancies. To study this, requires that the diffusion of the RE cations in the BaTiO 3 lattice must be modeled.

The Effect of RE Ion Diffusion
The literature on cation diffusion in BaTiO 3 is sparse. It is commonly assumed that the diffusion rate for titanium is much less than that for barium. However, work by Lee and coworkers [36,37] shows the diffusion rates of the two cations must be similar. Some experimental values for the diffusion activation energy in BaTiO 3 and SrTiO 3 exist. For A site vacancy diffusion there is the work of Garcia-Verduch and Lindner [38] (3.9 eV-Ba in BaTiO 3 ), Koerfer et al. [39] (5.6 ± 1.2 eV-Sr in BaTiO 3 ), Mayer and Waser [40]   there is the work of Preis and Sitte [45] (3.9 ± 0.7 eV-Ti in BaTiO 3 ); Koerfer et al. [39] (5.1 ± 0.6 eV-Zr in BaTiO 3 ). Although there is variation in these values they do not support the notion that vacancy diffusion rates on the A and B site lattices are very different. It should, however, be remembered that these are activation energies and could contain an energy contribution to create a mobile vacancy (release a vacancy from being bound to an impurity for example) as well as the migration energy of the mobile vacancy.
Several attempts have been made to calculate the diffusion migration energies. These are summarized in Table 1. All authors use static lattice methods to calculate the migration energy although ref. [26] made use of temperature-accelerated molecular dynamics (TAD) to identify the transition. Most authors used the nudged elastic band (which assumes that both the start and end points of the transition are known). Ref. [43] uses a mode-following method (RFO) which requires only the initial state. All results are defect energies (see Section 2). The results are clearly sensitive to the models used for the interatomic interactions (see the column for strontium vacancies in SrTiO 3 ). The results (for both simulation and experiment) are too scattered to permit a meaningful comparison except to note that (where results exist), the simulations predict a large difference in diffusion rate between A-site and B-site vacancy diffusion whereas the experimental results show the values should be similar for both the A and B-sites.
Our metadynamics results are shown in Table 2 below. The migration energies for intrinsic diffusion on the A and B sublattices are very similar-consistent with experiment although the migration energies are higher than the two experimental results (ref. [39] for A site diffusion; ref. [45] for B site diffusion but note the wide error bars). For all the REs, the simulations show that the lowest energy migration path involves alternating between the A and B sublattices. In this pathway, the rate-determining step is the B site -A site barrier crossing. Here the migration energy increases in the order Dy 3+ < Y 3+ < Gd 3+ , implying that Dy 3+ diffusion will be the fastest and Gd 3+ diffusion the slowest. Figure 4 shows the landscape of the vacancy diffusion mechanism in the three cases for Dy 3+ diffusion. Although all these RE ions can fit into both A and B sites, they are somewhat too small to fit comfortably into the A site and too large for the B site. It should be noted, however, that the B site is preferred for the RE ions due to the shorter cation-oxygen separations.
To understand why Dy 3+ has the lowest barrier we need to examine the shape of the free energy surfaces (Figure 5). The smaller size of Y 3+ compared to Dy 3+ means that it can fit better into the B site and therefore the well is deeper for diffusion from the B site to the A site. For the A sites a key effect is that the RE cations are too small for this site (reducing the coordination from 12 to 9) and therefore "rattle" in the site, generating local vibrations that can lead to diffusion. As we go to the larger Gd ion, this rattle is reduced as the ion begins to fit the site better and therefore the diffusion barrier is increased because the minimum is better defined. Dy 3+ shows the lowest migration barrier because it does not fit well into either site (a similar effect is known in intercalation compounds [46] ).
We can examine the effect of these different activation energies on the dispersion of RE cations into the grains of a BaTiO 3 ceramic using a simple diffusion model based on a solution to Fick's second law (full details can be seen in Appendix S9 in the Supporting Information). The calculated fractional distributions of the RE cations with their depth into the grain are shown in Figure 6.
We can immediately see a substantial difference between these distributions. Dy 3+ shows an almost homogeneous distribution across the first micron while the Y 3+ concentration drops much more rapidly over the first micron. The Gd 3+ concentration decays extremely rapidly and is effectively 0 at 0.5 Adv. Funct. Mater. 2020, 30,1905077 [25] TAD; AKMC; (NEB) forcefield from ref. [43] 3.95 Immobile a) -Immobile a) [22] NEB  micrometers. The change in the activation energy leads to significant differences in the diffusion constant that produce these very different profiles. These suggest that Dy 3+ will be able to pass through BaTiO 3 giving a uniform distribution of ions throughout the microstructure. The very poor diffusion of Gd 3+ will leave it largely trapped at the surfaces. Thus, despite its ability to distribute itself on both A and B sites, it will fail to disperse throughout the material and trap oxygen vacancies. Y 3+ could produce the desired dispersal but would require higher concentrations, higher temperatures, longer sintering times, or other sintering aids to make it as effective as Dy 3+ . In addition, more of the cations would be in the less effective defect traps arising from the scheme of exclusive B-site doping.

Conclusions
We have used simulations to demonstrate why Dy 3+ is the most effective RE ion to use for improving the lifetime of BaTiO 3based capacitors. Its ability to trap oxide vacancies is a combination of defect chemistry (thermodynamics) and diffusion behavior (kinetics). The size of the ion means that the dominant way by which the Dy 3+ ion is incorporated into BaTiO 3 is by self-compensation. Both the isolated Dy Ti ′ defect and the Dy Dy Ba • Ti − ′ pair are effective traps for oxide vacancies and form in high concentration. However, this is not sufficient to explain why Dy 3+ is better at preventing leakage currents than other rare earths like Gd 3+ . This requires a consideration of the kinetics of diffusion of rare earth ions in BaTiO 3 and here ion size also plays a role. The size of Dy 3+ is such that it fits poorly onto both the A site and the B site and this enables it to diffuse much faster than other rare earths and so produce a more uniform distribution of defects through the ceramic. Smaller quantities of Dy 3+ are sufficient to trap oxide vacancies throughout the material. Other RE ions could also do this, but would require higher sintering temperatures, higher concentrations, longer times, or different processing protocols to work.
Our use of complementary simulation methods, together with high-quality experimental data has resolved a longstanding conundrum in the defect chemistry of a ubiquitous and important electronic material. These methods could be applied to tackle similar issues that are found in many functional materials where the complete picture must be examined to understand the largescale behavior.

Computational Details
We have used both classical lattice statics and metadynamics molecular dynamics to model the defect behavior of BaTiO 3 . The potential set used for barium titanate was that of Freeman et al. [27] with RE potentials for Gd and Y from Lewis and Catlow. [47] A new Dy-O potential was fitted in-house to the Dy 2 O 3 structure. [48] Full details of the parameterization are given in Appendix S1 of the Supporting Information.
Static lattice defect simulations were performed using the GULP code (version 4.2) [49] using the Mott-Littleton approximation [50] whereby the volume of crystal surrounding the defect or cluster of defects is subdivided into two nested spherical regions. In the inner, defect-containing, region (I), the interactions between the atoms are calculated explicitly. In the outermost region (II) the interaction of the ions with the defect is based on the defect charge alone using a continuum method. The interactions between region (I) and the outer region (to a radius given by the sum of the radius of the inner region and the cut-off distance for the potential) are calculated explicitly. The region diameters used in all Mott-Littleton calculations were 21 Å for region 1 and 33 Å for the inner part of region II. These simulations give the internal energy at constant volume for the defect process which is a good Adv. Funct. Mater. 2020, 30,1905077   approximation to the constant pressure enthalpy and will be referred to as the "defect energy". [35] Both trapping (binding) energies and the activation energies (i.e., estimated barriers) for escape of the oxide vacancy were calculated as a function of the separation of the RE and oxide vacancy. Saddle points for an oxygen ion moving between vacant lattice sites were calculated using the rational functional optimization method [51] (RFO) which uses an uphill search method to identify the point at which only one eigenvalue of the hessian matrix is negative and is thus a saddle point. The difference in energy between a vacancy at a specific site and the saddle point is the activation energy, E A .
The migration rate of defects, particularly cation defects, is too low to use a simple molecular dynamics simulation and therefore we used the metadynamics technique to obtain the free energy migration pathway and hence the migration rate. Briefly, the method forces the simulation to explore new configurations defined by collective variables (which play the role of reaction coordinates in chemistry). They can be any simulation coordinates/parameters (e.g., atomic coordinates, interatomic separations, the energy of the system). A term (the bias) is added to the potential energy (and therefore to the interatomic forces) to prevent the system from revisiting a configuration already encountered. The system is therefore forced to visit all the configurations defined by the set of collective variables. The total bias potential added at each point can be used to construct a free energy surface in the space defined by the collective variables (e.g., along a diffusion pathway or reaction coordinate). This biasing enables the simulation to model infrequent events beyond the ability of standard molecular dynamics simulations but assumes that we can identify the right collective variables (reaction coordinates) for the process under discussion before we perform the simulation. A full introduction to the technique can be found in the review of Laio and Gervasio. [52] Simulations were performed using the DL_POLY Classic (version 1.9) [53] code in association with the PLUMED (version 1.3) plug-in. [54] An 8 × 8 × 8 supercell (2560 ions) was used. Each simulation was force-minimized and then equilibrated for 5000 timesteps (for a total of 2.5 ps) using the NPT ensemble (the relaxation times for the Nosé-Hoover thermostat and barostat were set to 0.10 and 0.05 ps, respectively). Production simulations were run for up to 6 ns using the NVT ensemble at a temperature of 298 K. The corners of the cell were pinned by fixing the positions of four atoms to prevent unphysical translations of the unit cell. Position coordinates (the absolute position of the moving ion referred to one or more of the Cartesian axes, x, y, z) were chosen as collective variables as this imposed no bias in the direction or pathway of migration (unlike using simply a distance between two ions). The Gaussian widths were chosen to correspond to the vibrational amplitude associated with the dopant in its potential well. Gaussian heights of either 0.05 or 0.10 eV were used. Gaussians were deposited every 300 timesteps (a time-step of 0.5 fs was used). Full details of the procedure are given in Appendix S2 in the Supporting Information.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.