20998
J. Phys. Chem. C 2009, 113, 20998–21007
Dynamics of Coupled Lithium/Electron Diffusion in TiO2 Polymorphs Sebastien Kerisit,* Kevin M. Rosso, Zhenguo Yang, and Jun Liu Pacific Northwest National Laboratory, Richland, Washington 99352 ReceiVed: July 8, 2009; ReVised Manuscript ReceiVed: October 22, 2009
Molecular dynamics simulations were performed with a potential shell model to investigate the diffusion of lithium ions and electron polarons in rutile and anatase. Simulations of an isolated lithium ion in rutile predict fast diffusion in the c channels with an activation energy of 0.05 eV, which corresponds to a jump rate of 4 × 1011 s-1 and a diffusion coefficient of 9 × 10-5 cm2 · s-1 at room temperature. In anatase, the activation energies for intra- and interoctahedron lithium hopping are 0.02 and 0.39 eV, respectively, and the lithium diffusion coefficient is 4-5 orders of magnitude slower than in rutile. When in the presence of an electron polaron, lithium hopping is predicted to be affected up to four hops away. The effects are more pronounced in rutile; whereby the first energy minimum along the c direction is absent due to the strong lithium-electron electrostatic interactions along the open c channels. Combining the lithium and electron polaron hopping rates, a coupled diffusion mechanism emerges whereby the electron polarons hop rapidly back and forth around the lithium ions. This process can lead to the occurrence of an instantaneous driving force for lithium hopping. The lithium ion-electron polaron binding energies are found to be large, with a stronger binding in rutile than in anatase, 0.45 and 0.28 eV, respectively, suggesting that, at low lithium mole fractions, lithium ions and electron polarons will form strongly correlated pairs. Introduction Titania (TiO2) is an attractive alternative for carbon-based anode materials in lithium-ion batteries due to its high surface area, chemical stability, and high theoretical capacity. Moreover, several polymorphs of titania have been shown to significantly benefit from nanostructuring or from the incorporation of nanoporosity. Increased rate capability,1 capacity,2 and tolerance for strain due to lithium insertion/extraction3 have all been reported for TiO2 nanomaterials. Although it has been recognized that nanostructuring and nanoporosity offer improved performance due to higher contact areas between electrolyte and electrode and shorter diffusion distances for lithium ions and electrons, there is a lack of fundamental understanding of the effects on charge and ion transport efficiency of nanosizing electrode-electrolyte interfaces, of the interdependence of electrode materials and electrolytes, contact areas, and diffusion distances, thus preventing reliable predictions of performance. Therefore, a long-term goal of our research is the discovery of design principles for nanostructured metal oxide electrodes that maximize charge and ion diffusivity. The first step in this process is to improve our understanding of charge and ion transport in electrode materials such as the titania polymorphs. In particular, two of the most fundamental questions, which we wish to address in this paper, are how is charge carrier diffusivity correlated to ion diffusivity at the level of individual electrons and ions and how does this coupling affect the overall transport efficiency and conductivity? In this paper, we concentrate on the two most stable and most studied TiO2 polymorphs, namely, rutile and anatase. Lithium insertion in both polymorphs to form LixTiO2 compounds has been studied extensively over the past few decades. Li insertion in macro- and microsized rutile crystals has consistently been shown to be limited at room temperature with published lithium * To whom correspondence should be addressed. E-mail: sebastien.kerisit@ pnl.gov.
mole fractions of x ) 0.0008,4 0.02,5 0.03,3,6 0.07,7 and 0.15.8 Higher mole fractions can be reached at higher temperatures. For example, Zachau-Christiansen et al.8 reported a Li mole fraction of 0.5 at 120 °C and Macklin and Neat9 achieved the insertion of one lithium ion per TiO2 unit at the same temperature. However, in the latter study, strong evidence for a phase transformation was also reported. In contrast to lithium insertion in macro- and microsized rutile, insertion in rutile nanoparticles is more extensive at room temperature with Li mole fractions of x ) 0.35,5 0.8,2 0.85,7 and 1.03,10 or greater.1 Similar results have been obtained by Wang et al.11 and Qiao et al.12 with mesoporous rutile and porous nanorods, respectively. It should be noted that, in all cases where nanoparticles were used as starting materials, various phase transformations were observed upon lithium insertion. Lithium insertion is more extensive in bulk anatase relative to bulk rutile at room temperature with reported lithium mole fractions of x ) 0.6,13,14 0.8,15 and 1.0.16,17 Importantly, it has been shown that lithium insertion in anatase induces a phase transformation18-20 to an orthorhombic phase21 and that phase separation20 between a Li-poor anatase phase and a Li-rich orthorhombic Li-titanate phase can occur. Wagemaker et al.22 subsequently showed that the two phases were in equilibrium with a continuous Li flux operating between them. As observed for rutile, anatase nanoparticles have been shown to incorporate lithium to a larger extent than their bulk counterparts.23-25 Neutron diffraction data obtained by Wagermaker et al.23,24 indicate that both the Li-poor anatase and Li-rich orthorhombic phases have higher Li capacities when the particle size reaches the nanoregime. In addition, for small enough particle sizes, they observed the occurrence of a new phase with a Li/Ti ratio of 1.24 Lithium insertion is accompanied by the addition of chargecompensating electrons in the titania lattice. In many transition metal oxides, electrons tend to be highly localized via strong interaction with phonons, forming small polarons.26-30 Polaron
10.1021/jp9064517 2009 American Chemical Society Published on Web 11/17/2009
Coupled Lithium/Electron Diffusion in TiO2 Polymorphs motion occurs by thermally activated hopping from one site to the next. Hence, in the case of very strong interaction, added electrons in lithiated titania polymorphs should manifest in the form of Ti3+ cations, and this model has proven useful for predicting temperature-dependent electron mobilities in these materials.31-34 Consistent with this model, X-ray photoelectron spectroscopy shows the formation of distinct Ti3+ states upon Li insertion in anatase as evidenced by the appearance and increase of a new peak in the Ti 2p spectrum.35,36 Additionally, Richter et al.37 evaluated the amount of charge transfer between lithium ions and the Ti 3d states to be 0.85 e, although this result was sensitive to the cross section of the ionization of a Ti 3d electron, which was not well-constrained. However, the amount of charge transfer derived experimentally is in good agreement with electronic structure calculations at the HartreeFock (HF)38 (0.79 e), density functional theory (DFT) using hybrid functionals (e.g., B3LYP)39 (0.69 e), and DFT/GGA19 (0.86 e) levels of theory. Because these localized chargecompensating electrons create relatively electronegative Ti cations, binding of Li+ to these sites is possible and the motion of Li+ through the lattice during (de)lithiation could in principle be strongly coupled to motion of electron polarons, and vice versa, both of which are thermally activated. However, this potential ambipolar diffusion aspect of Li insertion has largely been deemphasized in previous experimental work on Li diffusivity given that, while it is conceptually straightforward to measure lithium hopping rates in these materials, the converse is true for electron polarons. Measured diffusion coefficients of lithium in rutile and anatase vary significantly between studies owing, for the most part, to the difference in the form of the materials and the spatial scale at which the diffusion is probed. Macroscopic diffusion coefficients in anatase films40-45 vary from 10-10 to 10-17 cm2 · s-1. Kavan et al.46 reported diffusion coefficients of ∼10-13 cm2 · s-1 in anatase single crystals. Wagemaker et al.20,23 probed Li diffusion at the microscopic scale and found diffusion coefficients on the order of 10-12 cm2 · s-1. The data on Li diffusion in rutile is scarcer but consistently shows a high anisotropy between diffusion in the c channels (∼10-6 cm2 · s-1) and the ab plane (∼10-15 cm2 · s-1).4,47-49 Dependence of the lithium diffusion on the orientation of single crystals and therefore on the crystal faces exposed to the electrolyte was also reported.50 In parallel to the extensive experimental effort briefly summarized above, many computational studies have been carried out to determine the energetics of the different insertion sites in both rutile47-49,51-56 and anatase.19,38,47,51,54-57 These studies show that two sites are available for lithium insertion in rutile with tetrahedral and octahedral symmetries. The tetrahedral site is, however, much higher in energy than the octahedral site.48 Koudriachova et al.48 showed that lithium can adopt two positions when in the octahedral site, an on-center position with two short and four long Li-O bonds and an offcenter position with four short and two very long Li-O bonds. In anatase, lithium can occupy the vacant octahedral site. At x ) 0.5, Koudriachova et al.19 found that the lithium ion is in 5-fold coordination with Li-O distances from 1.97 to 2.05 Å. There are therefore two symmetrically equivalent positions inside the oxygen octahedron and quasi-elastic neutron scattering data presented by Wagemaker et al.58 showed that the lithium ions hop between the two positions on a picosecond time scale. Koudriachova et al.47,48,51 also calculated the open circuit equilibrium voltage profiles of rutile and anatase as a function of Li content and found that they correctly reproduce the
J. Phys. Chem. C, Vol. 113, No. 49, 2009 20999 principal features of the experimental discharge curves of TiO2-Li metal cells.8,9 Computational studies of lithium diffusivity in rutile47-49,59-61 and anatase38,47,57,58 have also been reported. In a series of papers, Koudriachova et al.47-49 studied lithium diffusion in rutile using the plane-wave pseudopotential approach and the spin-polarized generalized gradient approximation. The approach was to perform a series of static calculations at different points along the c direction and the ab plane. They found energy barriers of 0.04 and 0.8 eV, which correspond to diffusion coefficients of 10-6 and 10-14 cm2 · s-1, in the c direction and ab plane, respectively. They also analyzed the electronic structure as a function of Li content48 and found that ∼50% of the extra charge was localized on neighboring titanium ions and up to 25% was on oxygen and lithium ions. Other predictions of the lithium charge in rutile include the HF calculations of Stashans et al.56 (0.7-0.8 e) and Mackrodt55 (0.95 e) and the variable-charge calculations of Gligor and de Leeuw59 (0.45 e). Kingsbury et al.60 and Ajayi et al.61 calculated the activation energy for lithium diffusion in the c channels with potential models but found largely different values: 0.11 and 0.75 eV, respectively. To the best of our knowledge, the only computational study that made use of dynamic simulations to model lithium diffusion in rutile is that of Gligor and de Leeuw.59 Gligor and de Leeuw carried out molecular dynamics simulations with a variable-charge model based on the model of Streitz and Mintmire.62 They found that the lithium ions remained within their initial c channels. In other words, their simulations predict a negligible amount of diffusion in the ab plane, an observation consistent with the static calculations of Koudriachova et al.47-49 and experimental data.4 At low Li concentrations, they obtained a diffusion coefficient of ∼10-6 cm2 · s-1 with an activation energy of 0.22 eV. Early computational work on Li diffusion in anatase includes a study by Lunell et al.,38 who used unrestricted periodic Hartree-Fock (UHF) with the linear combination of atomic orbital approach and semiempirical intermediate neglect of differential overlap (INDO) method to calculate the energy barrier for lithium hopping between two octahedral sites in anatase. They found energy barriers of 0.60 eV at x ) 0.5 with the UHF method and an increasing energy barrier with the INDO approach from 0.51 to 0.56 eV between x ) 0.0625 and 0.5. In addition, they found that the unpaired electron occupied a previously empty d orbital of a neighboring Ti atom. However, symmetry constraints were applied in these calculations whereby the Li ion was constrained to remain equidistant to four Ti atoms. Koudriachova et al.47 reported an activation energy of ∼0.5 eV at x ) 0.5 and 1 eV at x ) 0.75 using, as previously described for rutile, a plane-wave pseudopotential approach without symmetry constraints. Tielens et al.57 also used a planewave pseudopotential approach to calculate the energy barrier for lithium hopping between octahedral sites. Surprisingly, they found that the activation energy decreases with increasing Li content. It is not clear why the trend is different from the other two computational studies. Li diffusion in anatase was also investigated with a potential shell model by Olson et al.,63 who calculated activation energies between 0.45 and 0.65 eV for values of x e 0.1. The only computational study making use of dynamic simulations is that of Wagemaker et al.;58 however, they focused solely on the intraoctahedron hopping. In this paper, we employ molecular dynamics techniques using a potential shell model32,64 based on the Matsui and Akaogi (MA) model64 to investigate the coupled lithium/electron polaron diffusion in anatase and rutile lattices. Here we specifically focus
21000
J. Phys. Chem. C, Vol. 113, No. 49, 2009
Kerisit et al.
on assessing the prospective role of ambipolar ion/electron electrostatic coupling on the dynamics of net diffusion. Although not a fully established mechanism, coupled lithium ion/electron polaron diffusion was shown by Ellis et al.65,66 to take place in another material relevant to Li-ion batteries, namely LixFeO4, and was confirmed computationally by Maxisch et al.67 The molecular dynamics approach employed here for electron polaron transport dynamics is based on our previous successful implementation of an umbrella sampling strategy applied to a two-state model of polaron transport developed for a number of transition metal oxides.32,68,69 The potential model at the heart of this work is conceptually simpler than the electronic structure calculations reviewed above, and it has been shown to be able to reproduce the crystal structures of a range of TiO2 polymorphs.32,64,70 In particular, Swamy et al.70 showed that the MA potential successfully reproduced the structure of TiO2 polymorphs that were not included in the original derivation set and that had titanium coordination numbers ranging from 6 to 9, whereas the polymorphs used for the derivation were only examples of 6-fold coordinated titanium. In addition to structural aspects, this potential model provides an accurate description of the relative stability of TiO2 polymorphs and, as shown by Dubrovinskaia et al.,71 was able to predict the existence of an intermediate phase at high pressure, which was confirmed experimentally. Furthermore, the potential model has been recently augmented to allow for the simulation of electron and hole polaron transfers. The use of a potential-based molecular dynamics approach offers several advantages. First, it allows us to perform dynamical simulations, for long periods of time, of lithium diffusion in TiO2 polymorphs, which have been particularly scarce.58,59 Second, it allows us to treat much larger lattices than used previously for the materials of interest. This, in turn, allows us to consider isolated lithium ions, electron polarons, and lithium-electron pairs and hence to probe lithium-electron interactions at essentially the infinite dilution limit. Finally, it provides a straightforward way to position and localize electron polarons on the individual sites of interest. Computational Methods In this section, we first introduce the potential shell model used throughout this work and the origin of the potential parameters. The models and simulation techniques used for computing the rates of electron polaron hopping and lithium hopping are then described. Details of the molecular dynamics simulation parameters are given in the last subsection. Potential Model. The calculations reported in this work are based on the Born model of solids,72 in which the atoms of a system are represented as point-charge particles that interact via long-range Coulombic forces and short-range interactions. The latter are described by parametrized functions and represent the repulsion between electron-charge clouds and the van der Waals attraction forces. In this work, the short-range interactions are described using a Buckingham potential and therefore the pairwise interaction energy takes the following form:
Uij )
( )
rij Cij 1 qiqj + Aij exp - 6 4πε0 rij Fij rij
(1)
The parameters (q, A, F, and C) used in this study are those optimized by Matsui and Akaogi64 for the TiO2 polymorphs rutile, anatase, brookite, and TiO2-II. In a previous publication, we derived a shell model version of the MA potential that allows for modeling electron small polaron hopping in these materials
TABLE 1: Core and Shell Charges Used in This Worka species 4+
Ti Ti3+ Li+ O2O2-(e)
core (e)
shell (e)
2.196 1.647 0.549 0.500 0.500
-1.598000 -1.673167
a The subscript (e) indicates the oxygen ions that belong to an isolated electron polaron. All oxygen core-shell spring constants, k, are set to 44.3 eV · Å-2. Core-shell interaction potential form: V ) k × rc-s2.
TABLE 2: Buckingham Potential Parameters Used in This Worka ion pair (ij)
Aij (eV)
Fij (Å)
Cij (eV · Å6)
Li-Li Ti-Li Ti-Ti Ti-O Li-O O-O
38533.955 33089.570 31120.528 16957.710 15465.549 11782.884
0.100 0.127 0.154 0.194 0.167 0.234
0.00 0.00 5.25 12.59 0.00 30.22
a
Buckingham potential form: Vij ) Aij × exp(-rij/Fij) - Cijrij-6.
and gives excellent agreement with the reorganization energies obtained from DFT calculations of several polaron transfers in titania polymorphs. A full description of the approach taken to derive these potential parameters is given in that publication.32 In the shell model,73 a polarizable ion is composed of two particles, a core and a shell, which share the ion’s charge and are linked by a harmonic spring, k: 2 Uc-s ) k × rc-s
(2)
where rc-s is the core-shell separation distance. The potential parameters used to describe the interaction between lithium ions and the TiO2 lattices were determined by fitting the potential parameters for Li to the lattice parameters, lattice constants, and bulk modulus of Li2O while keeping the MA potential parameters constant. The lithium partial charge was determined from the Li2O stoichiometry while keeping the oxygen partial charge from the MA model constant (i.e., qLi ) -qO/2 ) +0.549 e). Details of the derivation of the Li parameters and their performance are provided in another publication.74 All the potential parameters used in this work are summarized in Tables 1 and 2. Polaron Hopping Model. The aim of this section is to provide a brief description of the model used to describe small polaron hopping in this work. The reader is referred to refs 75 and 76 and references therein for a complete discussion of the model. Activated thermal hopping of electrons localized as small polarons is modeled as a valence interchange electron transfer reaction. In the polaron model, electrons introduced in a crystal lattice interact strongly with phonons and lower their energy by localizing at a lattice site, thereby creating a local distortion of the lattice. Polaron transfer is described by a two-state model whereby the reactants (before electron transfer) and products (after electron transfer) potential-energy surfaces are assumed to be parabolic with respect to the reaction coordinate, which is a function of the nuclear coordinates. Electron transfer occurs when thermal fluctuations cause the reactants and products states to be energetically degenerate. This conceptual model provides a basis for several theories, such as those of Holstein,27-30 Austin and Mott,26 and Marcus,77,78 which describe the transfer of polarons in various systems. In what follows, for simplicity, we use the terminology established by Marcus. A set of three
Coupled Lithium/Electron Diffusion in TiO2 Polymorphs
J. Phys. Chem. C, Vol. 113, No. 49, 2009 21001
P(X) [ P(〈∆E〉) ]
∆G(X) ) -RT ln
Figure 1. Free energy diagram of polaron hopping process modeled as an electron transfer reaction as a function of the energy gap.80-82
parameters is required to describe the kinetics and thermodynamics of electron transfer: the reorganization energy, the free energy of electron transfer, and the electronic coupling matrix element. The reorganization energy, λ, is the energy to distort the configuration of the reactants into that of the products, or vice versa, without changing the electronic distribution; the free energy of reaction, ∆G0, is the change in free energy upon electron transfer; and the electronic coupling matrix element, VAB, is the amount of electronic interaction between the reactants and products states at the transition state. The electron transfer parameters are illustrated in Figure 1. Electron transfer can take place in two different regimes: adiabatic and nonadiabatic, depending on the value of the adiabaticity criterion, which, in turn, depends mainly on the magnitude of VAB.79 When the coupling is strong, the electron transfer is adiabatic and the rate expression is as follows
ket ) υn exp(-∆G*/kBT)
(3)
where
∆G* )
(λ + ∆G0)2 - VAB 4λ
(4)
kB is the Boltzmann constant, T is the temperature, and νn is the effective vibrational frequency. However, if the electronic coupling is weak, the electron transfer is nonadiabatic and the rate constant takes the following expression:78
ket )
[
1 2π (∆G0 + λ)2 |VAB | 2 exp p 4λkBT √4πλkBT
]
(5)
To determine λ and ∆G0 using computer simulations, we follow an approach introduced by Warshel and co-workers80-82 whereby the free energy surfaces of the reactants and products states are calculated as a function of the energy gap, ∆E, which is the energy difference between the reactants and products charge distributions for a particular configuration of the system. Following Warshel’s work,80 the free energy curves are calculated from the probability of finding the system in a state for which ∆E ) X versus the probability of being in a state for which ∆E )〈∆E〉:
(6)
where 〈∆E〉 is the ensemble average of the energy gap. For most electron-transfer reactions, the activation free energy is too high for the region of phase space in the vicinity ∆E ) 0 to be significantly sampled with a single molecular dynamics calculation. To overcome this problem, we employ the umbrella sampling technique as first introduced, for an electron transfer reaction, in seminal papers by Hwang and Warshel,81 Kuharski et al.,83 and King and Warshel.82 In this method, simulations are carried out for a series of discrete values of the mixing parameter, θ, which allows for the construction of a new system potential from a linear combination of the reactants and products state potentials. This biases the system to generate configurations that have values of ∆E significantly different from 〈∆E〉 thereby allowing for a thorough sampling of the phase space along the reaction coordinate. The bias thus introduced is removed, once the simulation has been run, by adjusting the free energy.83 An important point to note is that electron transfer is a Franck-Condon process in that the electronic transfer is instantaneous relative to the nuclei motion, and hence, the system polarizability needs to be calculated for both the reactants and products states. As our model accounts for the polarizability of anions by means of a shell model, the position of the shells must be relaxed for each reactants and products state configuration in the calculation of the energy gap. The energy minimization of the shells is carried out using the steepest descent method. Finally, the determination of the third parameter, VAB, requires the use of molecular orbital electronic structure calculations. VAB is calculated using the approach introduced by Farazdel et al.:84
VAB )
|HAB - SAB(HAA + HBB)/2| 2 1 - SAB
(7)
where HAB is the interaction energy (HAB ) 〈ψA|H|ψB〉), SAB is the overlap between the reactant and product states (SAB ) 〈ψA|ψB〉), and H is the total electronic Hamiltonian. Cluster models excised from representative ∆E ) 0 configurations in the MD trajectories are used for this purpose. In this work, we only carried out a limited number of VAB calculations to supplement previously published values for rutile and anatase.31 Constrained Molecular Dynamics Simulations. The free energy profiles of lithium hopping between lattice sites were determined using the thermodynamic integration technique. In this approach, the free energy profile for lithium hopping is constructed from a series of simulations whereby a lithium ion is constrained at different points along the reaction coordinate. The reaction coordinate was defined by the vector formed by the final and initial positions of the lithium ion. In each simulation, the force along the reaction coordinate was removed to keep the lithium ion fixed in that direction but free to move normal to the reaction coordinate. The force required to keep the lithium ion fixed along the reaction coordinate, f(r), was recorded and used to determine the free energy profile via the following integration:
∆F(r) ) F(r) - F(r0) )
∫rr 〈f(r)〉 dr 0
(8)
where r0 is the initial lithium position and F(r0) is the zero of energy. Atomistic Simulations. All atomistic simulations were carried out with the computer program DL_POLY85 at 300 K and zero-applied pressure. Trajectories were generated in the NPT ensemble (constant number of particles, constant pressure, and
21002
J. Phys. Chem. C, Vol. 113, No. 49, 2009
Kerisit et al.
Figure 2. Free energy profiles for lithium hopping in rutile and anatase. In rutile, the reaction coordinate is the c direction. In anatase, the reaction coordinate is the vector defined by the centers of two neighboring oxygen octahedra.
constant temperature). In these simulations, the temperature and pressure were kept constant by use of the Nose´-Hoover thermostat86 and the Hoover barostat,87 respectively. The electrostatic forces were calculated by means of the Ewald summation method.88 A 9 Å cutoff was used for the short-range interactions and the real part of the Ewald sum. The Verlet leapfrog integration algorithm was used to integrate the equations of motion with a time step of 0.2 fs. The shells were given a mass of 0.2 au and their motion treated as that of the cores following the adiabatic shell model first introduced by Mitchell and Fincham.89 Results and Discussion Lithium Diffusion in Rutile and Anatase. We first considered lithium diffusion along the c-direction channels (referred to as “c channels” hereafter) in rutile. The simulation cell consisted of one lithium ion in a 5 × 5 × 8 supercell (400 TiO2 units). A series of 13 NPT simulations was carried out whereby the lithium ion was constrained in the c direction. A distance interval of 0.25 Å along the c direction was used for a total diffusion distance of 3.0 Å, which spans one unit cell. Each simulation was run for 200 ps after a 5-ps equilibration period. One titanium atom some distance away from the lithium ion was kept frozen to avoid any potential overall translation of the lattice to compensate for the force due to the lithium ion being displaced from its free energy minimum. The free energy profile obtained from integration of the force along the c direction is shown in Figure 2 and displays two free energy minima at positions r ) 0.0 and 0.5. These sites correspond to the two symmetrically equivalent octahedral sites. As already mentioned in the Introduction, Koudriachova et al.48 obtained two different coordination environments in the octahedral site with similar energies, namely, an on-center position with Li-O distances of 1.87 (×2) and 2.08 Å (×4) and an offcenter position with four oxygen nearest neighbors between 1.83 and 1.87 Å and two oxygen ions at 2.57 and 2.99 Å. The lithium-oxygen radial distribution function reveals that only the on-center position is occupied with two oxygen atoms at 1.83 and four at 2.08 Å. The calculated free energy barriers are 0.050 and 0.045 eV. The difference of 0.005 eV between the two equivalent barriers gives an indication of the uncertainty on the activation free energy due to the approach. The small activation energy for Li
diffusion along the c channels indicates that Li hopping is fast enough to be followed by direct simulation. Therefore, a series of NPT simulations was carried out whereby the same simulation cell was used but the lithium ion was not constrained and let free to diffuse. The temperature was varied from 300 to 500 K with 50 K intervals. The simulations were run for 500 ps with a 5-ps equilibration period. For each simulation, the lithium mean square displacement was used to determine its diffusion coefficient and generate an Arrhenius plot. From this plot, an activation energy of 0.051 eV was obtained with a preexponential factor of 2.85 × 1012 s-1 (R2 ) 0.991). As expected, the activation energy is in excellent agreement with the constrained simulations. The pre-exponential factor will be used in the calculations of Li hopping rates and diffusion coefficients. The predicted activation energy is in very good agreement with the work of Koudriachova et al.,47-49 who computed an activation energy of 0.04 eV from static DFT calculations. Other published activation energies, which were all obtained with potential models, vary quite significantly. Early work by Kingsbury et al.60 using a rigid-ion model predicted an activation energy of 0.11 eV. Later, Ajayi et al.61 obtained an activation energy of 0.75 eV with a polarizable model. More recently, Glidor and de Leeuw59 calculated the activation energy to be 0.22 eV from the temperature dependence of the Li diffusion coefficient at low Li concentrations using a charge-variable model. The range of activation energies published thus far highlights the importance of validating potential parameters against experimental data and/or electronic structure calculations. In our case, the close agreement between the predicted activation energy and the DFT calculations of Koudriachova et al.47-49 gives us confidence in our model. However, there is a large difference with the activation energy derived experimentally by Johnson4 from experiments with rutile single crystals (0.33 eV). It should be stressed that these experiments were done with millimeter-sized single crystals and that they therefore provide an effective macroscopic activation energy, although the reproducibility between samples suggests that the reported activation energy could be due to intrinsic properties of the rutile lattice. We will investigate in a subsequent section whether considering the diffusion of lithium in the presence of a chargecompensating electron improves the agreement with the experimental data of Johnson.
Coupled Lithium/Electron Diffusion in TiO2 Polymorphs
J. Phys. Chem. C, Vol. 113, No. 49, 2009 21003
Throughout this work, we calculated the diffusion coefficients using the following equation: n
D)
∑ miRi2ki i)1
(9)
2N
where n is the number of nonequivalent hops, mi is the number of accessible sites for hop i, R is the hopping distance, k is the hopping rate, and N is the number of dimensions in which diffusion takes place. For lithium diffusion along the c channels, we used n ) 1, mi ) 2, R ) 1.5 Å, and N ) 1 and obtained a diffusion coefficient of 9.3 × 10-5 cm2 · s-1 at 300 K (see Table 3). As expected, the same value was obtained with the direct molecular dynamics simulations. As mentioned in the Introduction, lithium diffusion in the ab plane has been shown experimentally4 and computationally47 to be much slower and was therefore not considered in this work. For anatase, a single lithium ion in a 6 × 6 × 3 supercell (432 TiO2 units) was considered. A series of 13 constrained molecular dynamics calculations was carried out to simulate lithium hopping between two octahedral sites. As was done for rutile, each simulation was run for 200 ps with a 5-ps equilibration period. Figure 2 shows that the free energy barrier for interoctahedron hopping is 0.39 eV. Wagemaker et al.20 reported a microscopic activation energy of 0.2 eV in the Lipoor anatase phase derived from 7Li nuclear magnetic resonance (NMR) measurements and a macroscopic activation energy of ∼0.5 eV for lithium insertion. They later showed that the macroscopic activation energy was related to that for lithium diffusion across the boundary between the Li-poor anatase and Li-rich titanate phases.22 Other macroscopic activation energies for lithium insertion and extraction were reported by Lindstro¨m et al. to be 0.35 and 0.38 eV for nanoporous films and 0.54 and 0.78 eV for chemical vapor deposited films from chronoamperometry43 with similar activation energies obtained with voltammetry.44 As mentioned in the Introduction section, electronic structure calculations, at the HF38 and DFT47,57 levels of theory, of lithium diffusion in Li0.5TiO2 anatase gave activation energies between 0.50 and 0.67 eV. In addition, activation energies between 0.45 and 0.65 eV were obtained with a potential shell model by Olson et al.63 for x e 0.1. The INDO calculations of Lunell et al.38 and the DFT calculations of Koudriachova et al.47 showed an increase of the diffusion barrier height with increasing Li mole fractions, whereas Tielens et al.57 and Olson et al.63 predicted the opposite trend. It is encouraging to note that, among the theoretically derived activation energies, our predicted value is the closest to the experimental work of Wagemaker et al.,20 which constitutes the only available activation energy obtained with an atomistic-level probe. In addition, a 1-ns MD simulation of a single lithium ion in anatase shows that the intraoctahedron hopping described by Wage-
maker et al.58 is reproduced by the potential shell model. Indeed, two minima within the oxygen octahedron were found separated by 0.8 Å along the c direction. Although this distance is shorter than the value reported by Wagemaker et al.58 (1.61 Å), the difference in temperatures between the two studies (exptl, 10 K; MD, 300 K) could explain in part the discrepancy between the experimental and calculated distances. Although the experimental activation energy for intraoctahedron hopping in lithium anatase was not available, the calculated activation energy of 17 meV lies, as expected, between the intraoctahedron barrier in lithium titanate and the interoctahedron barrier in lithium anatase.20,58 The diffusion coefficient of lithium in anatase was calculated using eq 9. The activation energy for Li interoctahedron hopping in anatase is too high to be observed with direct simulation. As a result, the pre-exponential factor could not be obtained from direct MD simulations and that obtained for rutile was used instead. The calculated interoctahedron jump rate at room temperature is 8 × 105 s-1, which translates to a correlation time of ∼1 µs, in agreement with the experimental value of Wagemaker et al.20 (47 µs). The parameters n ) 1, mi ) 4, R ) 2.54 Å, k ) 8 × 105 s-1, and N ) 3 yield a diffusion coefficient of 3.4 × 10-10 cm2 · s-1 at 300 K (see Table 3). The calculated diffusion coefficient is at the faster end of the experimental range (10-10-10-17 cm2 · s-1, see Introduction). Even though the calculated activation energy is higher than that derived experimentally by Wagemaker et al.,20 the calculated lithium diffusion coefficient is ∼2 orders of magnitude higher due to the much lower pre-exponential factor derived from the NMR experiments. In conclusion, the potential shell model used in this work gives satisfactory agreement with experimental data and electronic structure calculations when comparable data is available. The discrepancies with activation energies obtained from macroscopic experiments could be the result of other processes playing a role in the overall lithium diffusion. Finally, our model predicts that the diffusion of an isolated lithium ion in an otherwise pure TiO2 lattice is approximately 4-5 orders of magnitude faster in rutile than it is in anatase. Polaron Hopping in Rutile and Anatase. The calculations of the reorganization energy and free energy of reaction of several polaron hopping processes in rutile and anatase have been published in a separate paper.32 Therefore, in this section, we limit ourselves to briefly summarizing these results. However, since estimates of the polaron hopping rates were not included in that paper, we present their calculations in this section. As stated in the Computational Methods section, small polaron hopping was modeled as valence interchange electron transfer reactions between nearly completely localized electrons on the Ti sublattice. For both rutile and anatase, two types of electron transfer were considered, namely, those between edgeand corner-sharing TiO6 octahedra. The reorganization energies
TABLE 3: Activation Energies, Jump Rates, and Diffusion Coefficients of Isolated Lithium Ions and Electron Polarons in Rutile and Anatase electron/ion
direction
ET regime
νn (s-1)
∆G* (eV)
k (s-1)
electron electron ion
edge corner c channels
adiabatic nonadiabatic -
Rutile 2.4 × 1013 1.6 × 1012 2.9 × 1012
0.033 0.258 0.050
6.9 × 1012 7.8 × 107 4.1 × 1011
electron electron ion
edge corner interoctahedron
nonadiabatic nonadiabatic -
Anatase 1.5 × 1013 6.4 × 1012 2.9 × 1012
0.250 0.270 0.390
9.5 × 108 1.9 × 108 8.0 × 105
D (cm2 · s-1) 2.1 × 10-3 9.3 × 10-5 7.7 × 10-4 3.4 × 10-10
21004
J. Phys. Chem. C, Vol. 113, No. 49, 2009
Kerisit et al.
Figure 3. Lithium diffusion in the presence of an electron polaron in rutile (a) and anatase (b). Titanium atoms are shown in white, oxygen atoms in red, lithium atoms in green, and the titanium ion on which the electron polaron is localized in gray.
obtained at 0 K were 1.12 and 1.21 eV for the edge- and cornersharing transfers in rutile, respectively, and 1.15 and 1.24 eV for the edge- and corner-sharing transfers in anatase, respectively. The results were in excellent agreement with electronic structure calculations of Deskins and Dupuis31 performed at the DFT+U level of theory for the same four electron transfers. Then, umbrella sampling calculations were carried out to generate free energy diagrams. For each electron transfer, four values of θ were used to construct the free energy surfaces: 0.00, 0.17, 0.33, and 0.50. In each case, a 200 ps simulation was run with a 5-ps equilibration period and the energy gap distributions were obtained by collecting a configuration every 0.2 ps. The umbrella sampling calculations at 300 K gave activation energies of 0.23 and 0.26 eV of the edge- and cornersharing transfers in rutile and 0.25 and 0.27 eV for the edgeand corner-sharing transfers in anatase. The electronic coupling matrix elements obtained with small clusters treated at the HF level by Deskins and Dupuis31 were used to calculate the rates of electron transfer. Deskins and Dupuis31 showed that the edgesharing transfer in rutile has a VAB value of 0.20 eV, whereas the electronic coupling matrix elements of the other three transfers are much lower, between 0.01 and 0.03 eV. Diffusion coefficients were calculated using eq 9 and are shown in Table 3. In summary, the electron transfer between edge-sharing TiO6 octahedra in rutile is the only electron transfer in the adiabatic regime and is therefore by far the fastest electron polaron hopping process. However, since the corner-sharing electron transfer is rather slow, the electron polaron diffusion coefficient in rutile is only about three times higher than that in anatase. Comparing now with the lithium diffusion coefficients obtained in the previous section, it can be seen that electron diffusion is much higher in both materials with the difference being greater in anatase. In the next section, we evaluate to what extent the diffusion of isolated lithium ions and electron polarons is affected when both processes are coupled. Coupled Lithium/Electron Polaron Dynamics in Rutile and Anatase. For lithium diffusion in rutile, a series of 27 NPT simulations was performed with the lithium ion constrained in the c direction with 0.25-Å intervals, which corresponds to two unit cells. An electron polaron was placed on one of the nearestneighbor titanium ions of the initial position of the lithium ion. Figure 3 illustrates the simulation series. Each simulation was run for 200 ps after a 5-ps equilibration period. The resulting
free energy profile is shown in Figure 4. As expected the free energy profile shows a continuously uphill process due to the incomplete screening of the electrostatic interaction between the Li ion and the electron polaron. On the basis of the free energy profile in Figure 2, free energy minima are expected every 1.5 Å starting from position 0.0. Interestingly, the free energy minimum at position 0.5 is not present. At this distance, the electrostatic interaction between the lithium ion and the electron polaron is strong enough to make the free energy minimum disappear. Consequently, the diffusion barrier height increases significantly to 0.29 eV. It should be noted that it is now closer to the experimental value of 0.33 eV derived by Johnson.4 The barrier heights for the next two hops (0.13 and 0.08 eV) are also greater than that in the absence of an electron polaron and tend to the isolated case with increasing lithiumelectron distance. For lithium diffusion in anatase, a series of 36 NPT simulations was carried out with interval of 0.158 Å along the b direction and alternating (0.125 Å along the c direction. An electron polaron was localized on a nearest-neighbor titanium ion. Figure 3 illustrates the lithium diffusion path in the presence of an electron polaron. Figure 4 shows that the energy barrier increases by 36%, which is less significant that the increase predicted for rutile. In addition, unlike in rutile, all expected minima are actually observed. The next two hops show an increase of 3% and 8% with respect to the energy barrier in the absence of an electron polaron. Therefore, our calculations indicate that the presence of an electron polaron impacts the lithium diffusion to a greater extent in rutile than in anatase. We now consider the inverse case whereby an electron polaron diffuses in the presence of a fixed lithium ion. For rutile, we considered four consecutive edge-sharing polaron hops in the presence of a fixed lithium ion. As previously mentioned, the edge-sharing transfer takes place along the c direction. The four transfers correspond to a transport distance of ∼12 Å or 4 unit cells. Figure 5 shows an increase in activation energy for the first hop away from the lithium ion from 0.23 to 0.29 eV. The effect of the lithium ion diminishes with distance and reduces to 0.12 eV after four hops. The activation energy obtained for an isolated polaron, with the charge distribution used when a lithium ion is present, is 0.11 eV. This indicates that, by the fourth hop, the presence of the lithium ion does not affect the size of the activation barrier. The difference in activation energy between this case and that described in the
Coupled Lithium/Electron Diffusion in TiO2 Polymorphs
J. Phys. Chem. C, Vol. 113, No. 49, 2009 21005
Figure 4. Free energy profiles for lithium diffusion in the presence of an electron polaron in rutile (left) and anatase (right). Energy barriers are shown in red and the free energy relative to the initial position in black.
Figure 5. Free energy profiles for electron polaron diffusion in the presence of a lithium ion in rutile (left) and anatase (right). Energy barriers are shown in red and the free energy relative to the initial position in black.
previous section is due to the fact that, when in the presence of a lithium ion, some of the charge remains on the lithium ion, as described in the Computational Methods section. The reorganization energy is little affected by the presence of the lithium ion. It is predicted to differ only for the first hop away from the lithium ion and to converge to its bulk value for subsequent hops: 0.49 and 0.46 eV, respectively. For electron transfer in anatase, we also considered four edgesharing polaron hops in the presence of a fixed lithium ion. The four transfers correspond to a transport distance of ∼7.5 Å or 2 unit cells in the b direction. Figure 5 shows that the increase in activation energy for the first hop is similar to that calculated for rutile together with a similar free energy of electron transfer. However, unlike what was predicted for rutile, the change in free energy is comparatively small beyond the first hop. The difference between rutile and anatase could be explained by the fact that, in rutile, the lithium-electron electrostatic interactions are not as extensively screened as they are in anatase due to the absence of intervening atoms in the open structure of the c channels, as can be seen in Figure 3. However, as was observed for rutile, we find that the activation energy converges to that of the isolated case within about four hops. It should be noted that the MA potential underestimates the dielectric
constants of both rutile and anatase and that this might affect the distance at which the lithium and electron polaron hopping parameters converge to the isolated case. However, a similar conclusion was reached for charge transfer in FeO,68 whereby an oxygen vacancy was found to have an effect on the rate of charge transfer within a radius of about 7-8 Å, despite the fact that FeO was calculated to have a higher dielectric constant that rutile and anatase (∼12 for FeO compared to ∼3 for rutile and anatase), which was also in better agreement with experiment90 (∼27) than in the case of rutile and anatase. The reorganization energy also shows the same behavior as in rutile. As mentioned in the Computational Methods section, VAB values have been calculated and published for a series of symmetric electron transfers in bulk rutile and anatase.31 To begin to evaluate the effects of the presence of a lithium ion on the electronic coupling matrix elements, we calculated VAB for the first edge-sharing transfer away from the lithium ion in rutile. For this purpose, we excised a small cluster composed of two edge-sharing TiO6 octahedra and a lithium ion and capped it with hydrogen atoms following an approach used in previous publications.31,68,76 As a reference point, we computed VAB for the isolated edge-sharing electron transfer (i.e., without the lithium ion present) and obtained a value of 0.26 eV. Deskins
21006
J. Phys. Chem. C, Vol. 113, No. 49, 2009
Kerisit et al.
TABLE 4: Activation Energies and Jump Rates for Elementary Processes of the Coupled Lithium-Electron Diffusion in Rutile and Anatase electron/ion
direction
νn (s-1)
∆G* (eV)
k (s-1)
electron ion electron ion electron ion
first forward first forward first backward first backward second forward second forward
Rutile 2.4 × 1013 2.9 × 1012 2.4 × 1013 2.9 × 1012 2.4 × 1013 2.9 × 1012
0.163 0.287 0.059 0.027 0.083 0.128
4.4 × 1010 4.3 × 107 2.5 × 1012 1.0 × 1012 9.7 × 1011 2.0 × 1010
electron ion electron ion electron ion
first forward first forward first backward first backward second forward second forward
Anatase 2.1 × 1013 2.9 × 1012 2.1 × 1013 2.9 × 1012 2.1 × 1013 2.9 × 1012
0.279 0.534 0.030 0.355 0.146 0.403
4.4 × 108 3.0 × 103 6.7 × 1012 3.1 × 106 7.5 × 1010 4.8 × 105
and Dupuis31 showed, for the same electron transfer, that VAB converges to 0.20 eV with increasing cluster size. We will therefore apply a correction factor of 0.76 to account for the smaller cluster size used in this work. The transition state configuration was obtained using the linear synchronous transit approach84 and the corrected difference between the electronic coupling at the transition state and reactants configurations was found to be 0.12 eV. Our calculations therefore predict a reduction of 60% in VAB in the presence of a lithium ion. The calculated rate constants (see Table 4) allow us to determine the most likely mechanism for coupled lithium/ electron diffusion in the two TiO2 polymorphs. Starting from a lithium/electron pair in a nearest-neighbor configuration, we concentrate on the first hop that takes the two species apart, whether made by a lithium ion or an electron polaron. From this new position, we evaluate the likelihood for increasing the separation distance between the two species compared to that of a backward transfer that would bring the two species back to a nearest-neighbor configuration. In rutile, although the activation energy for lithium diffusion is similar to the diabatic activation energy for electron transfer, the strong electronic coupling in the c direction means that the rate of electron polaron hopping is much faster than that of lithium diffusion. The electron is therefore much more likely to move first. From there, the next three hops (electron backward, electron second forward, and lithium backward) have essentially the same rate. This suggests a mechanism whereby the electron hops back and forth until a backward hop by the lithium ion takes place resulting in the net ambipolar diffusion of the pair. In other words, the electrons can be thought as dragging the lithium ions through the lattice. In anatase, a similar picture emerges except that the larger difference between the lithium and electron backward hops will translate into a greater amount of back and forth hops by the electron before the lithium ion undergoes a successful hop (i.e., the dragging effect is weaker). In both cases, the calculated rate constants suggest that the electron polarons remain highly mobile even in the presence of lithium ions. Finally, the binding energy of a lithium ion/electron polaron pair was calculated from the energy difference between the isolated species and the lithium/electron pair in a nearestneighbor configuration. The binding energy was found to be large in both materials with a stronger binding in rutile than in anatase: 0.45 and 0.28 eV, respectively. The higher binding energy in rutile is likely due to the closer approach distance achieved in rutile with ) 2.30 Å in rutile versus 2.63 Å in anatase. For comparison, Maxisch et al.67 calculated
a value of 0.37 eV for the binding energy between a lithium ion and an electron polaron in FePO4. The binding energies in rutile and anatase suggest that, for low lithium mole fractions, the lithium ions and electron polarons are expected to be found in close proximity and to diffuse in pairs. Conclusions The dynamics of coupled lithium/electron polaron diffusion in rutile and anatase were studied using molecular dynamics simulations and a potential shell model. When available, the model and simulation techniques used in this work yielded good agreement with experimental data and electronic structure calculations. Simulations of an isolated lithium ion in rutile predicted an activation energy of 0.05 eV for hopping along the c channels between the octahedral sites, in agreement with the DFT calculations of Koudriachova et al.47-49 for Li0.5TiO2. In anatase, the activation energies for intra- and interoctahedron lithium hopping were found to be consistent with the NMR data of Wagemaker et al.20,58 Electron polaron hopping rates modeled as electron transfer reactions between edge-sharing and corner-sharing TiO6 octahedra were calculated for both polymorphs using published electron transfer parameters.31,32 The rate of polaron hopping along the c direction was found to be particularly fast due to a high electron coupling in this direction. Overall, electron polarons were predicted to exhibit a higher diffusivity than lithium ions in both polymorphs. However, calculations of the binding energies between lithium ions and electron polarons revealed a strong tendency to form pairs, which should have a significant effect on their diffusion. Indeed, a detailed analysis of the individual lithium hops and electron polaron hops that may take place starting from a lithium-electron nearest-neighbor configuration showed a strongly coupled ambipolar process. The calculated rates of electron transfer suggest that the electron polarons remain fairly mobile and will generally move first, which, in turn, leads to the creation of an instantaneous driving force for lithium hopping. At low lithium mole fractions, the electron polarons can therefore be thought as dragging the lithium ions throughout the lattice. It is our hope that detailed atomistic-level studies of coupled elementary processes of diffusion such as the one presented in this paper will lead to the development of a microscopic theory of ambipolar diffusion in electrode materials. Acknowledgment. The research described in this paper was conducted under the Laboratory Directed Research and Development Program at Pacific Northwest National Laboratory, a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy under Contract No. DE-AC0576RL01830. The computer simulations were performed in part using the Molecular Science Computing Facility in the William R. Wiley Environmental Molecular Sciences Laboratory, a national scientific user facility sponsored by the U.S. Department of Energy’s Office of Biological and Environmental Research and located at Pacific Northwest National Laboratory (PNNL). References and Notes (1) Jiang, C.; Honma, I.; Kudo, T.; Zhou, H. Electrochem. Solid-State Lett. 2007, 10, A127. (2) Hu, Y.-S.; Kienle, L.; Guo, Y.-G.; Maier, J. AdV. Mater. 2006, 18, 1421. (3) Baudrin, E.; Cassaignon, S.; Koelsch, M.; Jolivet, J.-P.; Dupont, L.; Tarascon, J.-M. Electrochem. Commun. 2007, 9, 337. (4) Johnson, O. W. Phys. ReV. 1964, 136, 284. (5) Kavan, L.; Fattakhova, D.; Krtil, P. J. Electrochem. Soc. 1999, 146, 1375.
Coupled Lithium/Electron Diffusion in TiO2 Polymorphs (6) Murphy, D. W.; Di Salvo, F. J.; Carides, J. N.; Waszczak, J. V. Mater. Res. Bull. 1978, 13, 1395. (7) Borghols, W. J. H.; Wagemaker, M.; Lafont, U.; Keller, E. M.; Mulder, F. M. Chem. Mater. 2008, 20, 2949. (8) Zachau-Christiansen, B.; West, K.; Jacobsen, T.; Atlung, S. Solid State Ionics 1988, 28-30, 1176. (9) Macklin, W. J.; Neat, R. J. Solid State Ionics 1992, 53-56, 694. (10) Reddy, M. A.; Kishore, M. S.; Pralong, V.; Caignaert, V.; Varadaraju, U. V.; Raveau, B. Electrochem. Commun. 2006, 8, 1299. (11) Wang, D.; Choi, D.; Yang, Z.; Viswanathan, V. V.; Nie, Z.; Wang, C.; Song, Y.; Zhang, J.-G.; Liu, J. Chem. Mater. 2008, 20, 3435. (12) Qiao, H.; Wang, Y.; Xiao, L.; Zhang, L. Electrochem. Commun. 2008, 10, 1280. (13) Whittingham, M. S.; Dines, M. B. J. Electrochem. Soc. 1977, 124, 1387. (14) Bonino, F.; Busani, L.; Lazzari, M.; Manstretta, M.; Rivolta, B. J. Power Sources 1981, 6, 261. (15) Ohzuku, T.; Takehara, Z.; Yoshizawa, S. Electrochim. Acta 1979, 24, 219. (16) Ohzuku, T.; Hirai, T. Electrochim. Acta 1982, 27, 1263. (17) Ohzuku, T.; Kodama, T.; Hirai, T. J. Power Sources 1985, 14, 153. (18) van de Krol, R.; Goossens, A.; Meulenkamp, E. A. J. Electrochem. Soc. 1999, 146, 3150. (19) Koudriachova, M. V.; de Leeuw, S. W.; Harrison, N. M. Phys. ReV. B 2004, 69, 054106. (20) Wagemaker, M.; Van de Krol, R.; Kentgens, A. P. M.; Van Well, A. A.; Mulder, F. M. J. Am. Chem. Soc. 2001, 123, 11454. (21) Cava, R. J.; Murphy, D. W.; Zahurak, S.; Santoro, A.; Roth, R. S. J. Solid. State Chem. 1984, 53, 64. (22) Wagemaker, M.; Kentgens, A. P. M.; Mulder, F. M. Nature 2002, 418, 397. (23) Wagemaker, M.; Borghols, W. J. H.; van Eck, E. R. H.; Kentgens, A. P. M.; Kearley, G. J.; Mulder, F. M. Chem.sEur. J. 2007, 13, 2023. (24) Wagemaker, M.; Borghols, W. J. H.; Mulder, F. M. J. Am. Chem. Soc. 2007, 129, 4323. (25) Sudant, G.; Baudrin, E.; Larcher, D.; Tarascon, J.-M. J. Mater. Chem. 2005, 15, 1263. (26) Austin, I. G.; Mott, N. F. AdV. Phys. 1969, 18, 41. (27) Holstein, T. Ann. Phys. 1959, 8, 325. (28) Holstein, T. Ann. Phys. 1959, 8, 343. (29) Friedman, L.; Holstein, T. Ann. Phys. 1963, 21, 494. (30) Emin, D.; Holstein, T. Ann. Phys. 1969, 53, 439. (31) Deskins, N. A.; Dupuis, M. Phys. ReV. B 2007, 75, 195212. (32) Kerisit, S.; Deskins, N. A.; Rosso, K. M.; Dupuis, M. J. Phys. Chem. C 2008, 112, 7678. (33) Nowotny, J.; Radecka, M.; Rekas, M. J. Phys. Chem. Solids 1997, 58, 927. (34) Yagi, E.; Hasiguti, R. R.; Aono, M. Phys. ReV. B 1996, 54, 7945. (35) So¨dergren, S.; Siegbahn, H. J. Phys. Chem. B 1997, 101, 3087. (36) Henningsson, A.; Rensmo, H.; Sandell, A.; Siegbahn, H.; So¨dergren, S.; Lindstro¨m, H.; Hagfeldt, A. J. Chem. Phys. 2003, 118, 5607. (37) Richter, J. H.; Henningsson, A.; Karlsson, P. G.; Andersson, M. P.; Uvdal, P.; Siegbahn, H.; Sandell, A. Phys. ReV. B 2005, 71, 235418. (38) Lunell, S.; Stashans, A.; Ojama¨e, L.; Lindstro¨m, H.; Hagfeldt, A. J. Am. Chem. Soc. 1997, 119, 7374. (39) Persson, P.; Bergstro¨m, R.; Ojama¨e, L.; Lunell, S. AdV. Quantum Chem. 2002, 41, 203. (40) Fu, Z.-W.; Qin, Q.-Z. J. Phys. Chem. B 2000, 104, 5505. (41) Ottaviani, M.; Panero, S.; Morzilli, S.; Scrosati, B.; Lazzari, M. Solid State Ionics 1986, 20, 197. (42) Kanamura, K.; Yuasa, K.; Takehara, Z. J. Power Sources 1987, 20, 127. (43) Lindstro¨m, H.; So¨dergren, S.; Solbrand, A.; Rensmo, H.; Hjelm, J.; Hagfeldt, A.; Lindquist, S.-E. J. Phys. Chem. B 1997, 101, 7710. (44) Lindstro¨m, H.; So¨dergren, S.; Solbrand, A.; Rensmo, H.; Hjelm, J.; Hagfeldt, A.; Lindquist, S.-E. J. Phys. Chem. B 1997, 101, 7717. (45) van de Krol, R.; Goossens, A.; Schoonman, J. J. Phys. Chem. B 1999, 103, 7151. (46) Kavan, L.; Gra¨tzel, M.; Gilbert, S. E.; Klemenz, C.; Scheel, H. J. J. Am. Chem. Soc. 1996, 118, 6716. (47) Koudriachova, M. V.; Harrison, N. M.; de Leeuw, S. W. Phys. ReV. Lett. 2001, 86, 1275. (48) Koudriachova, M. V.; Harrison, N. M.; de Leeuw, S. W. Phys. ReV. B 2002, 65, 235423.
J. Phys. Chem. C, Vol. 113, No. 49, 2009 21007 (49) Koudriachova, M. V.; Harrison, N. M.; de Leeuw, S. W. Solid State Ionics 2003, 157, 35. (50) Hengerer, R.; Kavan, L.; Krtil, P.; Gra¨tzel, M. J. Electrochem. Soc. 2000, 147, 1467. (51) Koudriachova, M. V.; Harrison, N. M.; de Leeuw, S. W. Solid State Ionics 2002, 152-153, 189. (52) Koudriachova, M. V.; Harrison, N. M.; de Leeuw, S. W. Comput. Mater. Sci. 2002, 24, 235. (53) Koudriachova, M. V.; de Leeuw, S. W.; Harrison, N. M. Chem. Phys. Lett. 2003, 371, 150. (54) Koudriachova, M. V.; Harrison, N. M.; de Leeuw, S. W. Solid State Ionics 2004, 175, 829. (55) Mackrodt, W. C. J. Solid State Chem. 1999, 142, 428. (56) Stashans, A.; Lunell, S.; Bergstro¨m, R.; Hagfeldt, A.; Lindquist, S.-E. Phys. ReV. B 1996, 53, 159. (57) Tielens, F.; Calatayud, M.; Beltran, A.; Minot, C.; Andres, J. J. Electroanal. Chem. 2005, 581, 216. (58) Wagemaker, M.; Kearley, G. J.; van Well, A. A.; Mutka, H.; Mulder, F. M. J. Am. Chem. Soc. 2003, 125, 840. (59) Gligor, F.; De Leeuw, S. W. Solid State Ionics 2006, 177, 2741. (60) Kingsbury, P. I.; Ohlsen, W. D.; Johnson, O. W. Phys. ReV. 1968, 175, 1099. (61) Ajayi, O. B.; Nagel, L. E.; Raistrick, I. D.; Huggins, R. A. J. Phys. Chem. Solids 1976, 37, 167. (62) Streitz, F. H.; Mintmire, J. W. J. Adhes. Sci. Technol. 1994, 8, 853. (63) Olson, C. L.; Nelson, J.; Islam, M. S. J. Phys. Chem. B 2006, 110, 9995. (64) Matsui, M.; Akaogi, M. Mol. Simul. 1991, 6, 239. (65) Ellis, B.; Perry, L. K.; Ryan, D. H.; Nazar, L. F. J. Am. Chem. Soc. 2006, 128, 11416. (66) Ellis, B.; Subramanya Herle, P.; Rho, Y.-H.; Nazar, L. F.; Dunlap, R.; Perry, L. K.; Ryan, D. H. Faraday Discuss. 2007, 134, 119. (67) Maxisch, T.; Zhou, F.; Ceder, G. Phys. ReV. B 2006, 73, 104301. (68) Kerisit, S.; Rosso, K. M. J. Chem. Phys. 2005, 123, 224712. (69) Kerisit, S.; Rosso, K. M. Geochim. Cosmochim. Acta 2006, 70, 1888. (70) Swamy, V.; Gale, J. D.; Dubrovinsky, L. S. J. Phys. Chem. Solids 2001, 62, 887. (71) Dubrovinskaia, N. A.; Dubrovinsky, L. S.; Ahuja, R.; Prokopenko, V. B.; Dmitriev, V.; Weber, H.-P.; Osorio-Guillen, J. M.; Johansson, B. Phys. ReV. Lett. 2001, 87, 275501. (72) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Oxford University Press: Oxford, UK, 1954. (73) Dick, B. G.; Overhauser, A. W. Phys. ReV. 1958, 112, 90. (74) Kerisit, S.; Rosso, K. M.; Yang, Z.; Liu, J. In preparation. (75) Iordanova, N.; Dupuis, M.; Rosso, K. M. J. Chem. Phys. 2005, 122, 144305. (76) Rosso, K. M.; Smith, D. M. A.; Dupuis, M. J. Chem. Phys. 2003, 118, 6455. (77) Marcus, R. A. J. Chem. Phys. 1956, 24, 966. (78) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265. (79) Brunschwig, B. S.; Logan, J.; Newton, M. D.; Sutin, N. J. Am. Chem. Soc. 1980, 102, 5798. (80) Warshel, A. J. Phys. Chem. 1982, 86, 2218. (81) Hwang, J. K.; Warshel, A. J. Am. Chem. Soc. 1987, 109, 715. (82) King, G.; Warshel, A. J. Chem. Phys. 1990, 93, 8682. (83) Kuharski, R. A.; Bader, J. S.; Chandler, D.; Sprik, M.; Klein, M. L.; Impey, R. W. J. Chem. Phys. 1988, 89, 3248. (84) Farazdel, A.; Dupuis, M.; Clementi, E.; Aviram, A. J. Am. Chem. Soc. 1990, 112, 4206. (85) Smith, W.; Forester, T. R. DL_POLY is a package of molecular simulation routines, copyright The Council for the Central Laboratory of the Research Councils, Daresbury Laboratory at Daresbury, Nr. Warrington., 1996. (86) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (87) Melchionna, S.; Ciccotti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533. (88) Ewald, P. P. Ann. Phys. 1921, 64, 253. (89) Mitchell, P. J.; Fincham, D. J. Phys.: Condens. Matter 1993, 5, 1031. (90) Seagle, C. T.; Zhang, W.; Heinz, D. L.; Liu, Z. Phys. ReV. B 2009, 79, 014104.
JP9064517
