Defect Properties and Lithium Incorporation in Li 2 ZrO 3

: Lithium zirconate is a candidate material in the design of electrochemical devices and tritium breeding blankets. Here we employ an atomistic simulation based on the classical pair-wise potentials to examine the defect energetics, diffusion of Li-ions, and solution of dopants. The Li-Frenkel is the lowest defect energy process. The Li-Zr anti-site defect cluster energy is slightly higher than the Li-Frenkel. The Li-ion diffuses along the c axis with an activation energy of 0.55 eV agreeing with experimental values. The most favorable isovalent dopants on the Li and Zr sites were Na and Ti respectively. The formation of additional Li in this material can be processed by doping of Ga on the Zr site. Incorporation of Li was studied using density functional theory simulation. Li incorporation is exoergic with respect to isolated gas phase Li. Furthermore, the semiconducting nature of LZO turns metallic upon Li incorporation.


Introduction
Lithium zirconate (Li 2 ZrO 3 or LZO) is a material of interest for many applications including electrode or electrolyte in Li-ion batteries [1][2][3], breeder blanket in nuclear reactors [4][5][6] and sorbent capture of CO 2 [7][8][9] due to its chemical and thermal stability. For battery applications, Li-ion conductivity should be high as battery performance is partly dependent on it during the charge-discharge process. Experimental work based on nuclear magnetic resonance (NMR) spectroscopy shows that the Li-ion conductivity in this material is high and the activation energy is calculated to be~0.51-0.65 eV [10]. Xu et al. [11] reported that the rate performance of LiNi 0.5 Co 0.2 Mn 0.3 O 2 can be improved by the surface coating of LZO. Electrochemical performance of LZO was examined experimentally and an excellent cyclic performance and rate capability reported [1]. In a density functional theory (DFT) simulation by Ferreira et al. [12], the relevance of vacancy-assisted Li-ion migration has been discussed.
Defects and dopants can influence a material's electrochemical property as diffusion of ions is directly related to the concentration of point defects [13,14]. Particularly, introducing dopants is a common way to form hyperstoichiometric dopant concentration of oxygen vacancies or interstitials, which in turn can impact oxygen self-diffusion as they effectively are the migrating defects [15]. Additionally, dopants can impact the migration energy barriers. In recent work, Kordatos et al. [16] performed DFT simulations to show that doping of aliovalent and isovalent cations on the Zr site is an efficient way to form oxygen vacancies and improve the diffusion and electronic properties of LZO. In the present study, atomistic simulations based on the classical pair-wise potentials to examine the defects, dopants and diffusion of Li-ions in LZO. In particular, we studied isovalent dopants on the Li and Zr sites for stabilizing LZO and aliovalent dopants on the Zr site for creating additional Li in LZO. In order to examine the incorporation, we used DFT simulations.

Computational Methods
Classical atomistic simulation was performed to calculate the defect energies, Li-ion migration pathways and solution energies for a variety of dopants. The generalized utilized lattice program (GULP) [17] was used. This method is based on the short-range (repulsive Pauli and attractive van der Waals) and long range (Coulombic) ionic interactions between ions. Short range interactions were modelled using Buckingham potentials (refer to Table 1) [18][19][20]. This potential method allows two unbonded atoms to interact and models the short-range interaction as a function of the intermolecular distance between those atoms. The Mott-Littleton method [21] was used to model point defects and migrating ions. In this method, the crystal surrounding the defect is divided into two regions (region I and region II). In region I, forces exerted by the defects are strong. Therefore, the ions in this region are relaxed explicitly. As forces in region II are relatively weak, quasi-continuum methods are used to relax the ions. Polarization of ions was modelled with the core-shell approach.
In previous work, we reported the methodology of calculating migration pathways and activation energies of migration in detail [22][23][24]. Table 1. Buckingham potential parameters used [18][19][20]. Incorporation of Li was studied using DFT simulation as implemented in the Vienna ab initio simulation program (VASP) code [25]. This code is based on projected augmented wave (PAW) potentials [26] and plane wave basis sets. We used a plane wave basis set with a cut-off of 500 eV in all calculations. For bulk LZO, we used 4 × 8 × 4 Monkhorst-Pack [27] k-points. A regular grid of k-points are distributed homogeneously in the Brillouin zone and sampled in the Monkhorst-Pack method. Reciprocal unit cell is first converted from the real space and the k-points are sampled based on the reciprocal basis vectors. A 2 × 2 × 2 super cell containing 192 atoms was used to model Li incorporation. For this supercell, 2 × 4 × 2 Monkhorst-Pack k-points were used. The generalized gradient approximation (GGA) parameterized by Perdew, Burke, and Ernzerhof (PBE) [28] was used to model exchange correlation. The conjugate gradient algorithm [29] was used to relax both atomic positions and cell parameters. Short range temporary dispersive interactions were modelled using zero damping DFT + D3 as parameterized by Grimme et al [30].

Structure of LZO
Li 2 ZrO 3 exhibits a monoclinic crystallographic structure with space group of C2/c (no. 15). Heiba et al. [31] reported the crystal structure of LZO with experimental values are a = 5.4089 Å, b = 9.0309 Å, c = 5.4144 Å, α = γ = 90 • and β = 112.50 • . The crystal structure has ZrO 6 octahedral units, which are interconnected by sharing their edges and corners (see Figure 1). Buckingham potentials used in the classical simulation and PAW potentials utilized in the DFT simulation were validated by performing geometry optimization calculations and comparing the relaxed lattice parameters with corresponding experimental values. There is a reasonable agreement between simulations and the experiment (refer to Table 2) as explained by Oberkampf et al. [32] Energies 2021, 14, 3963 3 of 11

Defect Energetics
In this section, we calculate the formation energies of point defects (vacancies and interstitials) using a classical simulation. Frenkel and Schottky defect energies were then calculated by combining point defect energies and appropriate lattice energies. Anti-site defects were also calculated. This defect is important as it can influence the diffusion property of a material. The following reaction equations were constructed to describe the intrinsic defect processes using the Kröger-Vink notation [33].
Li/Zr antisite (cluster) : The defect energies are tabulated in Table 3. The lowest defect formation energy is calculated for the Li Frenkel (1.09 eV/defect) as reported in a previous DFT simulation [13]. The Li-Zr anti-site cluster defect energy is higher only by 0.37 eV compared with the Li-Frenkel. In this defect, a small amount of Li and Zr atoms would exchange their positions simultaneously. However, this defect process would not take place spontaneously and would be dependent on the temperature and synthetic conditions. The difference in energy between the isolated and cluster forms of the anti-site defect is the binding energy. Here, the binding energy is calculated to be -1.42 eV. The presence of the anti-site defect has been reported in many experimental and theoretical studies [34][35][36]. This means isolated defects form clusters without any energy penalty. The Li 2 O Schottky defect is the lowest energy process among other Schottky processes. The formation of Li 2 O can destabilize the lattice. Such a destabilization can degrade the battery performance. The Zr Frenkel is higher in energy (by 3.59 eV) than the O Frenkel. This is due to the high formation energy of introducing Zr vacancy with +4 charge. Table 3. Defect formation energies of intrinsic defect process calculated in monoclinic LZO.

Defect Process Equation Number
Defect Energy (eV)

Diffusion of Li-Ions
Here, we discuss the diffusion of Li-ions and their activation energies. Classical atomistic simulations [37] were used to examine the Li-ion diffusion pathways and migration energies. Diffusion property of a materials is partly important in constructing a promising Li-ion battery. The identification of the diffusion ion pathways through experimental investigation is often challenging. Computational modelling has been shown to be a powerful tool to determine the migration pathways and activation energies [38][39][40]. Five different local Li hops (A-E) were identified, as shown in Figure 2. In Table 4, we report the Li hop distances and activation energies. Figure 3 shows the energy profile diagrams plotted for Li hops.   Figure 2. Li-ion diffusion pathways. Green, yellow, blue, grey and pink color atoms correspond to different Li hops. The hop C exhibits the lowest activation energy. A long range diffusion pathway (C→C→C→C) consisting of hop C was constructed. This pathway has a pattern of zigzag and Li-ion moves in the bc plane (see Figure 2). The difference in the activation energies calculated for hops A (0.65 eV) and B (0.63 eV) is very small. While there is a long range diffusion possible for hop B (B→B→B→B) in the ac plane with an overall activation energy of 0.63 eV, the diffusion of hop A is limited along the b axis. Another possible long range diffusion pathway (D→D→D→D) can be constructed by connecting local hops D in the bc plane. However, its activation energy is higher by 0.41 eV than that calculated The hop C exhibits the lowest activation energy. A long range diffusion pathway (C→C→C→C) consisting of hop C was constructed. This pathway has a pattern of zig-zag and Li-ion moves in the bc plane (see Figure 2). The difference in the activation energies calculated for hops A (0.65 eV) and B (0.63 eV) is very small. While there is a long range diffusion possible for hop B (B→B→B→B) in the ac plane with an overall activation energy of 0.63 eV, the diffusion of hop A is limited along the b axis. Another possible long range diffusion pathway (D→D→D→D) can be constructed by connecting local hops D in the bc plane. However, its activation energy is higher by 0.41 eV than that calculated for the lowest migration pathway (C→C→C→C). Finally, hops E were connected. In this long range diffusion (E→E→E→E), activation energy is relatively high. In a Li-NMR study carried out by Baklanova et al. [10], the activation energy for Li-ion migration was estimated to be 0.50-0.65 eV. This is in good agreement with the values calculated in this study. In a DFT simulation by Ferreira [12], a range of activation energies were reported with respect to reactants and products. Activation energies calculated from the products and reactants in the vacancy-assisted Li-ion diffusion are 0.651 eV and 0.749 eV respectively [12]. These values agree with those reported in this study, although DFT simulations treated Li vacancies as neutral in contrast to the current simulation where all defects were modelled as fully charged ions.

Solution of Dopants
Modification of material properties can be achieved through cation doping. In this strategy, isovalent doping requires no charge-compensating defect, whereas aliovalent doping can introduce point defects (vacancies and interstitials) in the lattice. The degree of disorder of a material is influenced by the concentration of point defects. Furthermore, point defect concentration increases exponentially with an increasing temperature [41]. Here we examined a variety of dopant ions using classical simulation. Table S1 in the Supplementary reports the Buckingham potentials used for dopants in this study (see Table S1).

Isovalent Dopants
First, doping of isovalent dopants (M = Na, K and Rb) on the Li site was considered as explained by the following reaction equation: Figure 4a reports the solution energies with respect to the radii of the dopant ions. Exothermic solution energy of −0.10 eV is calculated for Na + . The Li + has an ionic radius of 0.76 Å. The favorability of Na + can be partly due to its ionic radius (1.02 Å), closer to that of Li + . Solution energies are positive for both K + and Rb + . The highest solution of energy of 4.98 eV is calculated for Rb + , meaning that this dopant is highly unlikely for the doping process under normal temperatures. and reactants in the vacancy-assisted Li-ion diffusion are 0.651 eV and 0.749 eV respectively [12]. These values agree with those reported in this study, although DFT simulations treated Li vacancies as neutral in contrast to the current simulation where all defects were modelled as fully charged ions.

Solution of Dopants
Modification of material properties can be achieved through cation doping. In this strategy, isovalent doping requires no charge-compensating defect, whereas aliovalent doping can introduce point defects (vacancies and interstitials) in the lattice. The degree of disorder of a material is influenced by the concentration of point defects. Furthermore, point defect concentration increases exponentially with an increasing temperature [41]. Here we examined a variety of dopant ions using classical simulation. Table S1 in the Supplementary reports the Buckingham potentials used for dopants in this study (see Table S1).

Isovalent Dopants
First, doping of isovalent dopants (M = Na, K and Rb) on the Li site was considered as explained by the following reaction equation: Figure 4a reports the solution energies with respect to the radii of the dopant ions. Exothermic solution energy of −0.10 eV is calculated for Na + . The Li + has an ionic radius of 0.76 Å. The favorability of Na + can be partly due to its ionic radius (1.02 Å), closer to that of Li + . Solution energies are positive for both K + and Rb + . The highest solution of energy of 4.98 eV is calculated for Rb + , meaning that this dopant is highly unlikely for the doping process under normal temperatures. Next, tetravalent dopants (M = Si, Ge, Ti, Sn and Ce) were considered for doping on the Zr site. The following reaction equation describes the doping process.

MO + Zr → M + ZrO
Exothermic solution energies are calculated for Ge 4+ , Ti 4+ and Sn 4+ . The most favorable dopant is the Ti 4+ with a solution energy of −1.85 eV. Solution energy for Si 4+ is 0.18 eV. Next, tetravalent dopants (M = Si, Ge, Ti, Sn and Ce) were considered for doping on the Zr site. The following reaction equation describes the doping process.
Exothermic solution energies are calculated for Ge 4+ , Ti 4+ and Sn 4+ . The most favorable dopant is the Ti 4+ with a solution energy of −1.85 eV. Solution energy for Si 4+ is 0.18 eV. The endoergic solution energy for this dopant is owing to the ionic radius of Si 4+ (0.40 Å) deviating much from that of Zr 4+ (0.72 Å). The Ce 4+ exhibits the largest solution energy of 1.31 eV.

Aliovalent Dopants
The formation of extra Li + ions in LZO can be of interest to enhance the capacity of LZO for both battery and nuclear applications. Doping of trivalent cations (Al, Ga, Sc, In, Y, Gd, and La) was considered on the Zr site. This doping process will create Li interstitials based on the following reaction equation: Ga is the promising dopant for this process (refer to Figure 5). Solution energies calculated for Sc and In are higher only by 0.09 eV and 0.19 eV respectively compared

Li Incorporation in LZO
Simulations based on DFT were employed to examine the Li incorporated geometries and their electronic structures. Calculated total density of states (DOS) plot of bulk LZO is shown in Figure 6. Previous DFT simulation by Duan [42] shows that LZO is a wide gap semiconductor (band gap of 3.90 eV) which is in agreement with the calculated band gap of 3.60 eV in this study (see Figure 6). The net magnetic moment of the relaxed structure of bulk LZO is zero, as evidenced by the total DOS plot where spin-up states are the mirror image of spin-down states about the x-axis.  We have subsequently incorporated four Li atoms. Relaxed structures are shown in Figure 7. Incorporation is exoergic with respect to gas-phase Li atoms and endoergic with respect to Li metal (see Table 5). This is because of the higher endoergic dissociation energy to form Li atoms from bulk Li than the exoergic incorporation energy of gaseous Li atoms. The first two Li atom incorporation is highly negative compared with that of the third and fourth. Bader charge analysis [43] shows that Li atoms are positively charged and form strong bonds with oxygen in the lattice. Bader charge analysis uses partitioning method to calculate electronic charges on individual atoms in the lattice. Using Zero flux surfaces, atoms are divided and charge density is partitioned. There is a gradual increase in volume upon each Li incorporation.  Calculated total DOS plots of Li-incorporated structures are shown in Figure 8. Incorporation shifts the Fermi level towards the conduction band. This is because of an additional electron created by each Li incorporation. In all cases, Fermi level is occupied by states associated with Li and the resultant composites became metallic.

Summary
Defect properties, diffusion of Li-ion, solution of dopants and Li incorporation were investigated using classical and DFT simulations. Calculations find that the dominant defect in this material is the Li-Frenkel, ensuring the formation of Li vacancies. The next most favorable defect is the Li-Zr anti-site defect cluster. Li-ion diffusion is noted along the c axis and its activation energy is 0.55 eV. The isovalent dopants on the Li and Zr are calculated to be the Na and Ti respectively. The promising dopant on the Zr site to create extra Li in the lattice is Ga. DFT simulations show that incorporation of Li is exoergic with respect to gas phase Li. The metallic nature of the resultant composite is confirmed by the formation of Li + ions and electrons in the lattice. The results reported in this study can be validated by experiments for the development of Li-ion batteries and tritium breeding blankets.