Influence of elastic strain on the thermodynamics and kinetics of lithium vacancy in bulk LiCoO2

TL;DR

This study uses first-principles calculations to show how elastic strain significantly impacts lithium vacancy behavior and ionic conductivity in LiCoO2, affecting battery performance.

Contribution

It provides a detailed analysis of strain effects on vacancy formation and migration in LiCoO2 using DFT and elasticity theory, highlighting their impact on battery efficiency.

Findings

Elastic strain can change ionic conductivity by an order of magnitude.

Strain influences vacancy formation and migration energies.

Elastic dipole tensor effectively predicts strain effects.

Abstract

The influence of elastic strain on the lithium vacancy formation and migration in bulk LiCoO2 is evaluated by means of first-principles calculations within density functional theory (DFT). Strain dependent energies are determined directly from defective cells and also within linear elasticity theory from the elastic dipole tensor (Gij) for ground state and saddle point configurations. We analyze finite size-effects in the calculation of Gij, compare the predictions of the linear elastic model with those obtained from direct calculations of defective cells under strain and discuss the differences. Based on our data, we calculate the variations in vacancy concentration and mobility due to the presence of external strain in bulk LiCoO2 cathodes. Our results reveal that elastic in-plane and out-of-plane strains can significantly change the ionic conductivity of bulk LiCoO2 by an order of magnitude and thus strongly affect the performance of Li-secondary batteries.

Equations

$$E(n_d,ϵ)$$

$$+\frac{\partial E}{\partial n_d}{n}_d+\sum _{i,j}\frac{{\partial }^{2}E}{\partial n_d\partial {ϵ}_{ij}}{ϵ}_{ij}{n}_d+\dots$$

$$=E_o+\frac{1}{2}\sum _{i,j}{C}_{ijkl}{ϵ}_{ij}{ϵ}_{kl}+n_d\left(E_d+\sum _{i,j}\frac{{\partial }^{2}E}{\partial {ϵ}_{ij}\partial n_d}{ϵ}_{ij}\right)+\dots .$$

$$=E_o+\frac{1}{2}\sum _{i,j}{C}_{ijkl}{ϵ}_{ij}{ϵ}_{kl}+n_d\left(E_d+\sum _{i,j}\frac{\partial {\sigma }_{ij}}{\partial n_d}{ϵ}_{ij}\right)+\dots .$$

$$=E_o+\frac{1}{2}\sum _{i,j}{C}_{ijkl}{ϵ}_{ij}{ϵ}_{kl}+n_d\left(E_d-\sum _{i,j}{G}_{ij}{ϵ}_{ij}\right)+\dots .$$

$$\Delta E=-\sum _{ij}{G}_{ij}{ϵ}_{ij}.$$

$${\sigma }_{ij}^d\equiv \frac{\partial E(n_d,{ϵ}_{ij})}{\partial {ϵ}_{ij}}=\sum _{kl}{C}_{ijkl}{ϵ}_{kl}-n_d{G}_{ij}={\sigma }_{ij}^0-n_d{G}_{ij},$$

$$G_{ij}=-\frac{\partial E_d}{\partial {ϵ}_{ij}}=-{\frac{\partial {\sigma }_{ij}}{\partial n_d}|}_{{ϵ}_{ij}}=-\frac{1}{n_d}({\sigma }_{ij}^d-{\sigma }_{ij}^0)=-V_0\Delta {\sigma }_{ij},$$

$$E_d({ϵ}_{ij})=E_d(0)-\sum _{ij}{G}_{ij}{ϵ}_{ij}.$$

$$E_d({ϵ}_{ij})$$

$$E_b(ϵ)=E^S({ϵ}_{ij})-E_d({ϵ}_{ij}),$$

$$E_b({ϵ}_{ij})$$

$$=E_b(0)-\sum _{ij}\Delta G_{ij}{ϵ}_{ij},$$

$$\displaystyle D\approx{\rm exp}\bigg{(}{-\dfrac{\Delta G_{\rm A}}{k_{B}T}}\bigg{)},$$

Full text

I Abstract

The influence of elastic strain on the lithium vacancy formation and migration in bulk LiCoO2 is evaluated by means of first-principles calculations within density functional theory (DFT). Strain dependent energies are determined directly from defective cells and also within linear elasticity theory from the elastic dipole tensor ($G_{ij}$) for ground state and saddle point configurations. We analyze finite size-effects in the calculation of $G_{ij}$, compare the predictions of the linear elastic model with those obtained from direct calculations of defective cells under strain and discuss the differences. Based on our data, we calculate the variations in vacancy concentration and mobility due to the presence of external strain in bulk $\rm LiCoO_{2}$ cathodes. Our results reveal that elastic in-plane and out-of-plane strains can significantly change the ionic conductivity of bulk LiCoO2 by an order of magnitude and thus strongly affect the performance of Li-secondary batteries.

Stress/strain, defect thermodynamics and kinetics, elastic dipole tensor, finite-size-effect, ionic mobility

II I. Introduction

The presence of strain fields can significantly influence the efficiency and lifetime of functional materials such as semiconductors, solar cells and Li-ion batteries Sun et al. (2007); Penev et al. (2001); Hubbard et al. (2008); Sheu et al. (2005); Zhang et al. (2011); Deshpande et al. (2010); Tealdi et al. (2016); Malik et al. (2010); Wang et al. (2007); Diehm et al. (2012). In principle, strain can be internally generated by structural defects within the bulk and/or interfaces or can be induced by external loads. These strains or corresponding stresses can cause lattice deformations and distortions and therefore also affect the formation and migration of point defects. As a result, the conductivity of an ionic conductor can be significantly changed Walukiewicz (1994); Azpiroz et al. (2015); Fisher et al. (2008); Freedman et al. (2009); Goyal et al. (2015). This coupling is most relevant in Li-ion batteries, where due to charging/discharging processes, bulk and interfacial as well as thermal strains can occur. The fact that induced stresses, which raise during intercalation, could weaken the interface between electrode/electrolyte and finally degrade the battery performance has been extensively discussed in the past Moradabadi and Kaghazchi (2015, 2016); Hoang and Johannes (2016, 2012, 2014); Guo et al. (2016); Chung et al. (2002); Ning et al. (2014); Bower et al. (2011); Garcia et al. (2005); Islam and Fisher (2014); Okubo et al. (2007). Much less is known, however, about the coupling of strain fields to the formation and migration of Li vacancies in the bulk part of the material.

As one of the first, Choi $et$ $al.$ Choi and Pyun (1997) have examined the effect of intercalation-induced stress on Li migration in LiCoO2 using experimental techniques such as electrochemical impedance spectroscopy together with cyclic voltammetry. They showed that the mismatch strain between intercalated and deintercalated states in LiCoO2 is due to a phase transition between the $\alpha$ and $\beta$ phases. This results in intercalation-induced stress in the phase boundary region. Since an inserted Li ion causes structural disorder in the host material, it introduces a strain field which affects the next intercalated Li and leads to a strain-induced elastic interaction between all intercalated Li ions. This elastic interaction has both short and long range effects and it was also found that as the particle size decreases, the stress gradient across LiCoO2 particles increases.

Transmission electron microscopy (TEM), X-ray diffraction and electrochemical strain microscopy (ESM) have revealed that upon delithiation the c-axis in LixCoO2 (x=0.5) expands by about 2% Wang (1999); Takahashi et al. (2007) and rather large stresses occur during charging/discharging process at different charge/discharge rates (c-rates) Diercks et al. (2014); Kim et al. (2004); Mukhopadhyay and Sheldon (2014); Balke et al. (2010). Garcia $et.$ $al.$, for example, have shown based on a continuum model that large stress values (up to $\pm$ 200 MPa) are generated at particle contacts and these values increase with increasing c-rate of discharge Garcia et al. (2005). Xiong $et$ $al.$ have shown in an ab-initio based study that during deintercalation of LixCoO2, the $c$-axis increases up to 3.25% Xiong et al. (2012). Li $et~{}al.$ also reported that compressive stress raises while the Li concentration is decreasing and at x=0.5, the measured strain at the O-Co-O octahedral slabs is 4.8% Li et al. (2009). A mathematical model was developed by Renganathan $et$ $al.$ to reveal the mechanical stresses generated during the discharge process in carbon and LiCoO2 Renganathan et al. (2010). Their findings show that at high discharge c-rates, stresses also increase and the stress caused by phase transformations is related to the amount of each phase present in the electrode. They also concluded that particle size and distribution can affect the generated stress Renganathan et al. (2010). Critical rates of charging and particle size below which fracture of LiCoO2 occurs were predicted by Zhao $et$ $al.$ using a kinetic and fracture mechanics model. They have shown that as the discharge rate increases the particle size should decrease so that fracture is prevented (for c-rates more than 5 C, particle sizes less than 200 nm) Zhao et al. (2010).

While all these studies point to the fact that significant stress and strain levels can occur in cathode materials of Li-secondary batteries, their influence on Li-ion diffusion has hardly been studied. In a recent theoretical work Ning $et$ $al.$ Ning et al. (2014) have shown that by applying uniaxial tensile strain along the c-axis of bulk LiCoO2, the Li diffusion barrier decreases Ning et al. (2014), but did not derive consider the case of a more complex tensorial strain field.

In this study, we calculate the elastic dipole tensor Leibfried and Breuer (1978) in order to characterize the coupling of strain fields to the formation and migration energies of Li vacancy in bulk-LiCoO2. We obtain the components of the elastic dipole tensor for defect formation and also migration from total energy calculations within density functional theory. Similar calculations of the defect dipole tensor have been recently reported for defects in metals and silicon by Varvenne $et\,al.$ Varvenne et al. (2013) and in UO2 by Goyal $et\,al.$ Goyal et al. (2015). In order to demonstrate the degree of coupling between lateral and longitudinal components of the stress tensor in LiCoO2, various supercell geometries are studied. Moreover, we compare the predictions from linear elasticity theory with directly calculated formation and migration energies of vacancy under strain and show that the elastic dipole-tensor allows us to quantify the influence of strain in a computationally efficient manner. Finally, we estimate the effect of lateral and longitudinal strains on ionic conductivity in bulk LiCoO2.

III II. Theory and Computational methods

III.1 Elastic dipole tensor

The insertion of a point defect into a crystal produces local elastic distortions. Moreover, there will be an interaction between this defect and a stress or strain field present in the crystal. This is similar to the interaction of an electric dipole with an applied electric field. Therefore, a defect inducing local distortions is called an elastic dipole, which is -contrary to the electric dipole- characterized by a second-rank tensor. This elastic dipole tensor, which is also called double force tensor Leibfried and Breuer (1978), is the negative derivative of the defect formation energy $E_{d}$ with respect to an imposed bulk strain

$G_{ij}=-\dfrac{\partial E_{d}}{\partial\epsilon_{ij}},$

if, as usual, we ignore entropy contributions. Thus, $G_{ij}$ is relating the atomic structure of a point defect and its elastic field. In case of a purely dilatational strain the defect relaxation volume $\Delta V=\dfrac{1}{3B}{\rm Tr}\{G_{ij}\}$, where $B$ is the bulk modulus, can be directly obtained from $G_{ij}$. In principle, the concept of the elastic dipole tensor can be conveniently understood by expanding the free energy per volume in terms of the density of defects $n_{d}=N_{d}/V$ and the strain tensor $\epsilon_{ij}$ (both being intensive quantities) Leslie and Gillan (1985). If entropy contributions are neglected the expansion of the energy density (T=0 K) reads as:

[TABLE]

Here $E_{o}$ is the total energy of the non-strained defect-free system, $n_{d}$ the number of defects, respectively, $E_{d}$ is the formation energy of a defect, $\sigma_{ij}$ is the stress, $C_{ijkl}$ are the components of the stiffness tensor and $G_{ij}$ is the elastic dipole tensor. As such $G_{ij}$ describes the interaction of the defect with a strain field. The change in energy under strain that is exclusively due to the presence of a defect is given by

[TABLE]

Thus, for example, a positive lattice strain would lower the formation energy of a defect having a positive relaxation volume. The stress in a material under strain $\epsilon_{ij}$ with a defect density $n_{d}$ is eventually given by

[TABLE]

if we only consider first-order components of the elastic dipole tensor. Note that -depending on the stress definition- different sign conventions have been proposed in literature. Here we stick to the original one used by Leibfried and Breuer Leibfried and Breuer (1978).

With relation (3) it is straightforward to compute the components of the elastic dipole tensor numerically using atomistic methods. In a given periodic supercell, an individual defect is introduced. While the cell parameters are fixed, the atomic positions are relaxed and the induced stress is calculated. Then, the elastic dipole tensor can be obtained from the relation

[TABLE]

where $V_{0}$ is the volume of the supercell containing one defect, $\sigma_{ij}^{d}$ is the stress of the defective cell and $\sigma_{ij}^{0}$ is the stress of the defect-free cell (which in principle should be or very close to zero).

In supercell calculations, there will be a certain unwanted contribution to the energy described in Eq. 1 from the interaction between the defects in their periodic images. It is therefore necessary, either to correct for this interaction, or to increase the size of the repeating unit to make it negligible. In the case of charged defects, another contribution comes from the Coulomb interaction between the defects, which must be corrected as shown by Leslie and Gillan Leslie and Gillan (1985). Elastic contributions can usually be made small enough by increasing the cell size. If ab-initio methods are used, however, the accessible cell sizes are rather limited and elastic interactions can induce higher order effects, that are not covered by linear elasticity. As can be seen from Eq. 4, the variation of stress by a defect follows the relation $\Delta\sigma_{ij}=-G_{ij}/V_{o}$ and thus goes to zero in the dilute limit ($\lim_{V_{o}\rightarrow\infty}\Delta\sigma_{ij}=0$). Therefore, by plotting the components of $\Delta\sigma_{ij}$ as function of the inverse volume, one can test the convergence behavior of the calculated stresses. For small supercell sizes higher order effects will affect the linearity of $\Delta\sigma_{ij}$ with respect to the inverse volume. In practice, a polynomial fit including higher order terms are required. The coefficient of the linear term dominates at small values of inverse volume and thus still corresponds to the elastic dipole tensor.

III.2 Defect formation and migration energies under strain

The defect formation energy for a given homogeneous strain $\epsilon_{ij}$, considering Eq. 2, is then given by

[TABLE]

Note, that $E_{d}(\epsilon_{ij})$ can also be directly calculated as shown later.

In the present work, we only consider a neutral Li vacancy. The reason is that we want to disentangle the electrostatic and elastic interactions. The conventional way to calculate the formation energy of defects in the neutral state and in the presence of strain is to use the following equation Zhang and Northrup (1991)

[TABLE]

In this equation, the first term is the total energy of the strained system containing a single neutral Li vacancy, the second term is the total energy of the strained pristine system and the last term is the strain-free chemical potential of the Li reservoir. This can be compared to the defect formation energy as a function of strain calculated from Eq. 5 for a specific concentration. It should be noted here that the vacancy formation energy depends on the chemical potential of lithium in the reservoir, which is taking up the removed Li atom. Thus, in principle one has to consider the fact that strain might also affect the reservoir.

We now move to the coupling of defect migration and strain fields. The energy barrier that is required for an ion to jump between two sites is obtained by the energy difference between the saddle point and the initial configuration

[TABLE]

where $E^{S}$ is the energy of the defective system in the saddle point (transition state) and $E_{d}$ the energy of the defective initial state.

For a system under strain, the energy barrier can be computed by applying a particular strain to the simulation cell and direct calculation of this energy difference. Alternatively, one can also calculate the defect dipole tensor $G^{S}_{ij}$ at the saddle point configuration and obtain the strain dependent barrier from

[TABLE]

where $E_{b}(\epsilon_{ij})$ is the migration energy barrier in the strained material, $E_{b}(0)$ is the migration energy barrier in the unstrained material, while $\Delta$G${}_{ij}^{b}$ is the change of the elastic dipole tensor by going from the initial to the transition state. For calculating the $\Delta$Gij, we perform two single point calculations, one for the initial and one for the transition state in order to obtain $\Delta\sigma^{d}_{ij}$. Afterwards, using Eq. 4, the elastic dipole tensor for the initial and transition state can be calculated. The atomic coordinations used in the single point calculations are obtained from the NEB calculation for the unstrained case. Therefore, by performing only one NEB calculation at $\epsilon_{ij}=0$, the strain-dependent activation barriers for Li-diffusion can be investigated.

III.3 Model structures

Fig. 1 (a) shows the unit cell of bulk LiCoO2. Li atoms (violet color) are coordinated within the octahedrals of CoO*-2*. During the delithiation, Li vacancies are formed and start to migrate. Considering the layered structure of LiCoO2, we expect a significant influence of strain on the formation and migration of these vacancies. The strain can either originate from the expansion/contraction during the deintercalation/intercalation, respectively, or can be applied externally (e.g. interfacial strains from SEI, solid electrolyte or binder).

In Fig. 1(b), the direct mechanism of Li diffusion towards the single vacancy is indicated (or migration of vacancy towards Li). During this process, the Li ion must overcome an energy barrier $E_{b}$ by pushing the two nearby Li ions and passing through them. This mechanism is shown in Fig. 1(b) and (c) by gray arrows. Therefore, Li on the saddle point is under compressive stress (tensile on its neighbors). While there are also other Li migration mechanisms reported der Ven and Ceder (1999); Moradabadi and Kaghazchi (2015), in the following, we will focus on the strain dependence of the formation and migration of a neutral single Li vacancy moving on a direct pathway. Since we are only interested in the strain dependency, we deliberately study neutral cells in order to disentangle electrostatic image interactions of the charged defects from the elastic image interactions being also present.

III.4 Computational details

The calculations were performed using the local atomic-orbital DFT-code SeqQuest http://www.cs.sandia.gov/ paschul/Quest/ (2015) with norm-conserving pseudopotentials and the generalized-gradient approximation of Perdew, Wang and Ernzerhof (PBE) Perdew et al. (1996) for exchange and correlation. Spin optimization is performed during all geometry relaxations. The energy convergence criterion for all calculations is $1\times 10^{-5}$ eV. For smearing, we used the Gaussian method, together with a narrow width of smearing (0.005 eV) to make sure that the correct spin polarization can be achieved. Diffusion pathways are investigated using the nudge elastic band (NEB) method Jónsson et al. (1998) (as implemented in the SeqQuest code). A 24$\times$24$\times$5 Monkhorst-Pack $k$-point mesh for the $1\times 1\times 1$ unit cell is considered and for larger supercells, it is adjusted accordingly. For all calculations, grid spacing for the charge density integration is set to 0.16 $\AA$. Convergence tests with respect to k-point sampling and grid spacing show that calculated stresses are well converged (less than 0.0002 GPa). For the charge density difference calculation, the relaxed structures from SeqQuest calculations were used in the DFT-code VASP Kresse and Furthmüller (1996). The unit cells for all VASP calculations were optimized again while the energy-forces criteria and spin polarization were chosen similar to ones used in SeqQuest calculations.

IV Results and discussion

IV.1 I. Size dependence of elastic dipole tensor

($G_{ij}$ in dilute limit)

In order to capture finite-size effects, we used supercells with different volumes and aspect ratios. We also investigated the degree of coupling between lateral and longitudinal components by studying various cell geometries. Figure 2 (left) shows the diagonal elements of the calculated stresses for isotropically repeated $n\times n\times n$ ($1\leq n\leq 3$) supercells containing a single (neutral) Li vacancy and those of non-isotropic replicas $n\times n\times 1$ ($1\leq n\leq 5$, middle) and $n\times n\times 2$ ($1\leq n\leq 4$, right) together with the second-order fit of the stresses and the diagonal components of the elastic dipole tensor. For the isotropic case (Fig. 2 - left), due to computational limitations, the $4\times 4\times 4$ supercell is not calculated. In all cases, the non-diagonal elements exhibit small values. This is due to the fact that a Jahn–Teller distortion of MO2 octahedrals occurs Xiong et al. (2012); Zhang et al. (2013), which can affect these non-diagonal elements of the elastic dipole tensor. The separation of defects along the $z$-axis between periodic images for $1\times 1\times 1$ is equal to 13.96 Å, while across the $xy$-plane, it is equal to 2.85 Å. Therefore, we expect a stronger defect-defect interaction between the periodic images in the $xy$-direction rather than along the $z$-axis.

In all geometries, where the $z$-extension was varied (see Fig. 3), the $xx$ and $yy$ components of the elastic dipole tensor are comparable (about 1.2 eV), while the $zz$ component is strongly varying and is even changing its sign in case of the $1\times 1\times n$ cells indicating a strong coupling between $xx$, $yy$ and $zz$ components.

We determined the components of the elastic dipole tensor by fitting the relation $\Delta\sigma_{ij}=-G_{ij}/V+\beta/V^{2}$ to the calculated data. The second term is accounting for the fact that higher order contributions to the elastic dipole tensor might be significant. Using the mechanical definition of stress, we describe outward stress on the cell boundaries with positive and inward stress with negative signs. This means that the absence of positive (outward) stress would lead to cell contraction and vice versa ($\Delta V=\dfrac{1}{3B}{\rm Tr}\{G_{ij}\}$).

All data provide evidence for positive components of the dipole tensor in $x$- and $y$-direction and negative components in the $z$-direction. Therefore, the presence of a vacancy (during the deintercalation process) leads to a contraction in the $xy$-plane and expansion in the $z$-direction. These results are in agreement with the data reported by Xiong $et$ $al.$ Xiong et al. (2012). The results for the $1\times 1\times 1$ cells deviate from the expected scaling behavior, since due to the high defect concentration, non-linear contributions prevail. Thus, these data are not shown in the plots. It is also evident that the data for the non-isotropically replicated supercells exhibit significant non-linear contributions for the $xx$- and $yy$-components of the defect elastic dipole tensor.

In Figure 4 we compare the values of the diagonal components of $G_{ij}$ using histogram plots. These plots show the scaling behavior from lateral (or longitudinal) extension alone towards isotropic one. It can be seen that the $G_{xx}$/$G_{yy}$ and $G_{zz}$ components do not show a linear scaling behavior. This is due to the fact that in these relatively small cells image-image interactions are still strong (especially across the $xy$-plane).

Figure 5 reveals the origin of the strong coupling between Gxx-Gyy and Gzz components. For this, we have plotted the charge density differences ($\rho_{\rm defective}-\rho_{\rm pristine}$) for two supercell dimensions, namely 3$\times$3$\times$1 and 3$\times$3$\times$2. For both supercells, similar isosurface value of 0.001 $|e|/{\AA}^{3}$ is considered. Figure 5 indicates the charge distribution difference between defective and pristine structures as the lattice is increasing along $z$-axis. It can be seen that by lattice extension along the $z$-axis, the overlap of orbitals is decreasing which leads to the minimization of defect-defect interactions. However, due to the especial geometry of LiCoO2 unit cell, despite the 3 times repeated cell along the $xy$-plane, the two defects between nearby images still show significant electrostatic interactions (although the overlap of orbitals along the $xy$-plane is also decreasing but with a less progress) which can be translated to the rather strong coupling of $xy$-plane stress components.

To sum up this part, our results reveal that size and geometry effects have a massive influence on the calculation of the elastic dipole tensor because of the coupling between stress components. We see that the fits to the data scale more linearly as the cell sizes are increasing homogeneously ($n\times n\times n$). Therefore, non-linear terms affect the result the least for the isotropically scaled cell.

IV.2 II. Defect formation energy under strain

In order to investigate the effect of strain on the defect formation energy in LiCoO2, we used Eqs. 5 (dipole-tensor method) and 6 (direct method). With SeqQuest we calculated a value of 2.86  eV for the formation energy of a neutral Li vacancy in the 3$\times$3$\times$1 supercell as our reference for the unstrained case, which is in agreement with data previously reported by Hoang $et$ $al.$ Hoang and Johannes (2014).

Fig. 6 shows the change in defect formation energy as a function of external strain (laterally, in the $xy-$plane) with two methods and two DFT codes. For all cases, the red points are obtained from the direct calculation in which the total energies of pristine and defective are strained. In this case, we can directly calculate the formation energies under strain. Since we need to calculate the total energies in each strain regime, this is computationally time-consuming. The second method is shown with the blue points which are obtained from Eq. 5. In this case, we only need the total energies of unstrained pristine and defective supercells, together with the already-calculated $G_{ij}$ for that specific supercell.

The directly calculated strain-dependent defect formation energies show a slightly non-linear behavior and do not follow the prediction of linear elastic theory. Moreover, this non-linear behavior is not symmetric in case of 3$\times$3$\times$1 supercell using either DFT methods, as can be seen in Fig. 6.

A similar trend for other systems has also been previously reported. For example Zhu $et\,al.$ investigated the effect of strain on the formation energy of Cu vacancy in Cu2ZnSn(S,Se)4 system Zhu et al. (2014) in which they show how symmetry breaking will lead to this deviation from linear elasticity. Similarly Aschauer $et\,al.$ found a non-linear behavior for the strain dependence of oxygen formation in MnO Aschauer et al. (2015). The reason why linear elastic theory fails to predict this trend is first due to the fact that there are higher order terms relevant in the Taylor expansion given in Eq. 1. Another reason is that the elastic constants of the defective system are not the same as in pristine system. Regardless of the methods and supercell sizes, it can also be concluded that under the influence of lateral elastic strain, the formation energy of a single neutral vacancy in bulk LiCoO2 varies by about 0.02 eV (with $\epsilon=1\%$).

IV.3 III. Migration barrier of Li vacancy under strain

We now focus on the effect of strain on the Li diffusion for the pathway indicated in Figs. 1 (b) and (c). In the first section of results and discussion, we calculated the elastic dipole tensor in the dilute limit and revealed that size and geometry of the supercell have a significant influence. This is why we first investigated the influence of cell size on the migration barrier for Li vacancy jumps using the NEB method. The results are shown in Fig. 7. The calculated values for the 3$\times$3$\times$1 and 4$\times$4$\times$1 deviate by about 10%, while the smaller cell shows a significantly smaller migration barrier. In order to reduced computational efforts, we have chosen the 3$\times$3$\times$1 supercell in order to investigate the strain effect on migration energy barrier.

Figure 8 shows migration barriers for vacancy hopping as function of strain using the directly calculated data (red in the $xy$-plane) and those obtained from the elastic dipole tensor method (blue in the $xy$-plane and green along the $z$-axis) for 3$\times$3$\times$1 supercell using Eq. 8. For the direct calculations, four lateral strain states are introduced homogeneously on the lattice parameters of $a$ and $b$. The calculation of the elastic dipole tensor in Eq. 8 was performed at the initial and saddle points for the 3$\times$3$\times$1 supercell. Since we already calculated $G_{ij}$ for the initial state (Fig. 2), we used the geometry of the saddle point from the NEB calculation for the unstrained case to obtain the $G_{ij}$ for the transition state. We obtain the following values for the elastic dipole tensor in the initial state and at the saddle point:

$G_{\rm initial}=\begin{bmatrix}\ 1.66&0.01&0.01\\ 0.01&1.65&-0.02\\ 0.02&-0.02&-1.72\\ \end{bmatrix}$

$G_{\rm saddle}=\begin{bmatrix}\ -0.08&-0.81&1.48\\ -0.81&0.77&-2.30\\ 1.60&-2.49&8.24\\ \end{bmatrix}$

By considering the signs of diagonal components of elastic dipole tensor for the initial state, the direction of the gray arrows in Figs. 1-b and 1-c are justified.

Fig. 8 shows the very good agreement between the directly calculated data and those obtained from elastic-dipole tensor in the $xy$-plane. Due to computational limitations, we do not compare the two methods for the longitudinal stain, however, we postulate that the good agreement observed for the $xy$-plane strain is also established for the $z$-axis strain. It can also be seen that all curves behave linearly but with different slopes between lateral and longitudinal strains (between red-blue and green plots). By applying a positive strain, the energy barrier decreases and vice versa, which is also in agreement with previous results from Ning et. al. Ning et al. (2014). It can be seen that, as a result of tensile strain, there is more Li-Li separation and less Li-Li repulsion and therefore, moving Li atom can intercalate/deintercalate more easily. From another point of view, applying strain could disturb the octahedrals orientation and therefore affects the potential energy surface which directly influences the energy barrier.

Since we see a good agreement between two methods in Fig. 8 for the lateral strain, we also plotted the effect of longitudinal strain on the migration energy barrier which is indicated with green color. It is clear from Fig. 8 that the longitudinal strain has a much stronger effect on the energy barrier, compared to lateral strain. Therefore, this shows that instead of significant computational efforts due to performing NEB calculations at each strain regime, it is possible to use Eq. 8 and obtain the same results with much less computational resources.

A comparison between Figs. 6 and 8 also indicates that the effect of lateral ($xy$-plane) strain on the Li vacancy formation and migration are opposite of each other. This means that, while a lateral compressive strain decreases the energy barrier (due to above mentioned reasons), the same strain results in a larger Li vacancy formation energy. However, as can be seen from the plots, this effect on formation and migration energy is not equal. Therefore in the activation energy, which is the sum of formation and migration, at a lateral strain regime equals to $+1\%$, for example, the migration term is decreased by 0.04 eV while formation term is increased by 0.02 eV. Thus, at this strain value, the overall trend is a decrease in the activation energy by 0.02 eV. Considering the relation of diffusivity,

[TABLE]

decreasing of 0.02 eV in activation energy leads to almost five times increase in the diffusivity and hence the conductivity. Therefore, the massive strain effect on ionic conductivity in bulk LiCoO2 is evident. We again note that the effect of longitudinal strain is even more dominant than the lateral one. Particularly, only 1% strain along the $z$-axis can change the conductivity up to ten times compared to unstrained case. Moreover, according to theoretical findings Li et al. (2009), 4-5% and experimental evidences Wang (1999); Takahashi et al. (2007), 2-3% strain is expected during the lithiation/delithiation of LiCoO2 to Li0.5CoO2. Therefore, strain fields can significantly influence the ionic conductivity in bulk LiCoO2 which shows that it can be employed to tune the particle mobility in battery materials.

V Summary and CONCLUSIONS

In summary, we have performed a detailed comparison for the variations of defect formation and migration energies with respect to strain in bulk LiCO2 using (i) direct evaluation of strained supercells (at each strain regime separately) and (ii) classical elasticity theory by computing the elastic dipole tensor ($G_{ij}$). The latter method requires only three DFT calculations (pristine and defected cells for formation energies plus defected cell in transition state for energy barriers) for evaluation of both formation and migration energies to obtain the full variations as a function of any strain state.

We found that the calculated formation energies using the elastic dipole tensor method deviate slightly when they are compared with the direct method (by less then 1% for $\epsilon<0.015$). We note tha, however, deviation to some degree depends on particular basis sets and/or pseudopotentials. Moreover, the mentioned deviation largely cancels out for migration barriers and it makes this method even more error-free in case of strain-induced diffusion analysis. Therefore, using the elastic dipole tensor method for the analysis of strained-induced formation and migration energies is computationally very efficient.

We also highlight that estimating the elastic dipole tensor in the dilute limit can be affected by finite-size effects and coupling between the stress components. Moreover, we found that the contribution of migration energy to the total activation energy when a lateral strain regime is applied, is more dominant than the contribution of formation energy. Finally, we can conclude that the effect of even small strains on ionic transport properties in bulk LiCoO2 is very significant. We showed that only 1% in-plane strain can change the conductivity by a factor of 5, while the presence of out-of-plane strains can change the conductivity by an order of magnitude.

$\bf{Acknowledgments}$

AM and PK gratefully acknowledge support from the ”Bundesministerium für Bildung und Forschung” (BMBF), the computing time granted on the Hessian high performance computer ”LICHTENBERG” and Zentraleinrichtung für Datenverarbeitung (ZEDAT) at the Freie Universität Berlin. JR and KA acknowledge support through SPP 1473 of the German Research Foundation and project DFG-578/19-1.