Theory of anomalous voltage–discharge behavior for topotactic intercalation

Solid State Ionics 127 (2000) 163–168 www.elsevier.com / locate / ssi

Theory of anomalous voltage–discharge behavior for topotactic intercalation V.I. Kalikmanov*, M.V. Koudriachova, S.W. de Leeuw Department of Applied Physics, Delft University of Technology, Lorentzweg 1, 2628 CJ Delft, The Netherlands Received 6 April 1999; received in revised form 16 July 1999; accepted 26 July 1999

Abstract Anomalous voltage–discharge behavior observed in topotactic intercalation materials is explained in terms of a lattice–gas model for interstitial sites with short-range nearest-neighbor repulsive interactions. For sufficiently strong interactions a second order phase transition results in the formation of fine structure of intercalated ions at certain critical concentrations. As a result a voltage–discharge curve shows a step-like behavior with nonanalytic features at the critical concentrations.  2000 Elsevier Science B.V. All rights reserved. Keywords: Intercalation; Lithium ion battery; Lattice–gas models; Phase transitions; Open circuit voltage PACS: 64.60.Cn; 84.60.Hs; 05.70.Fh

An intercalation compound represents a matrix of host atoms and a certain amount of guest ions penetrated into it and occupying interstitial sites. A typical example of such a material is Li c MO 2 where lithium is intercalated into a metal–oxide matrix; c is a number of intercalated ions per metal atom M. In a number of materials c can be between 0 and 1 without substantial changes to the structure of the host matrix. This so called topotactic intercalation mechanism is the basis for application of intercalation compounds as an electrode in a rechargeble battery (see e.g. Refs. [1,2] and references therein). During a discharge cycle Li ions diffuse through electrolyte from anode to cathode while the electron

*Corresponding author. Tel.: 131-15-278-5596; fax: 131-15278-6081. E-mail address: [email protected] (V.I. Kalikmanov)

transport through electrolyte is completely prohibited. Guest ions accumulated in the host matrix during the discharge cycle are removed from it upon charging. In an ideal cell these two processes are reversible, and a large voltage difference between the electrodes can be achieved. Material research in this field is aimed at development of electrodes which would be able to provide high voltage in combination with low molecular weight. One of the most important characteristics of a battery is a voltage–discharge curve, also called the open circuit voltage (OCV), describing the equilibrium voltage difference V between electrodes as a function of the amount of intercalated material. If ze e is a charge of the intercalating ion, where e e is an electron charge (z 5 1 for a Li ion), then a simple energy balance gives: 2 ze eV 5 m 2 m0

0167-2738 / 00 / $ – see front matter  2000 Elsevier Science B.V. All rights reserved. PII: S0167-2738( 99 )00261-1

164

V.I. Kalikmanov et al. / Solid State Ionics 127 (2000) 163 – 168

where m is a chemical potential of an intercalated ion in a host matrix (cathode) and m0 is its value in the pure metal (anode). Hence, the OCV is fully characterized by dependence of m on concentration: 1 1 V(c) 5 2 ] m (c) 1 ] m0 ze e ze e

‘lattice’ refer to the lattice of interstitials). Then the energy of an arbitrary configuration of sites reads: E( p1 , . . . , pN ) 5 u

O pp i

j

(2)

(i, j )

(1)

It was observed experimentally that at certain conditions the OCV possesses an anomalous step-like behavior while in other cases it is a smooth function. In this paper we present a model of topotactic intercalation which explains these features. An interstitial ion represents a positively charged point defect in the host matrix. Compensating electron reduces the metal ion preserving electroneutrality. Coulombic repulsion between interstitial ions is screened by the host matrix. The nature of the metal and the strength of its M 41 /M 31 redox couple determines to a large extent the screening radius. It is expected that in the most of the intercalation materials screening is quite substantial [2,3] so that repulsive interaction between interstitials can be considered short-ranged. It is plausible therefore to simplify the description of a topotactic intercalation compound by replacing the two-component system ‘host matrix1guest ions’ by an effective one-component medium of particles, occupying interstitial sites, with short range repulsive interaction u . 0. In view of the preceding discussion u should be regarded as an effective interaction potential which includes effects due to the presence of the host matrix. In our model u is an input parameter to be provided by experiment or ab initio computer simulations; typically, u . 2 4 6k B T. In a simplified description only repulsion between nearest neighbors can be taken into account. These assumptions lead naturally to formulation of a lattice–gas model for the interstitials with nearest neighbor repulsive interactions. About the lattice we assume that it is bichromatic, i.e. it can be partitioned into two interpenetrating sublattices so that nearest neighbors belong to different sublattices (this partitioning is possible e.g. for the square lattice, the simple cubic- and body-centered cubic lattice). Frustration effects (occurring in the triangular and face centered cubic lattice) are not considered. With each site i (i 5 1, . . . , N) we associate a variable pi which is equal to unity if the site is occupied and zero otherwise (in what follows the terms ‘site’ and

The round brackets in the sum denote summation over nearest neighbors with each pair counted only once; let q be the coordination number of the lattice. This expression implies that all sites are equivalent. The number of occupied sites is n(h pi j) 5 o Ni 51 pi . It is more convenient to measure concentration of intercalated particles in terms of the average fraction of occupied sites: knl x5] N where k l denotes the thermal average. Since interactions are purely repulsive the equilibrium configuration results from the competition between energetic and entropic terms in the free energy. If repulsion is weak the entropic contribution prevails for all fractions, which means that ions occupy lattice sites at random. If repulsion is strong enough random occupation is favourable for fractions smaller than some critical x (1) c . Beyond x c the energy contribution becomes comparable to the entropic one which results in an disorder–order transition. Particles occupy preferably one of the sublattices up to x 5 1 / 2 when preferable occupation of the second sublattice starts. At large concentrations, x . x (c2 ) , in view of a large number of mutually repelling particles it becomes favourable to place them again at random in order to maximize entropy. So, at x (2) the order– c disorder transition takes place. The particle–hole ( 1) symmetry of the model yields x (2) c 5 1 2 x c . The phase transition is of second order and involves spontaneous symmetry breaking. The order parameter is the difference in occupation of sublattices. Since n can vary it is convenient to work in the grand canonical ensemble with the fixed values of the chemical potential m, the total number of sites N and the temperature T. The average number of occupied sites is given by the thermodynamic relationship [4]: 1 ≠V x 5 2 ] ]u N,T N ≠m

(3)

V.I. Kalikmanov et al. / Solid State Ionics 127 (2000) 163 – 168

where V( m,N,T ) is the grand potential. Instead of searching for m (x) we will study the inverse function x( m ) described by Eq. (3). The phase diagram can be obtained by studying the isothermal compressibility: 1 ≠x x 5 ]2 ] T x ≠m Its divergence signals about a phase transition. An important characteristic of the OCV is the derivative: dx ≠x 2 ] 5 ze e ] dV ≠m proportional to x 2 x . A phase transition results in T nonanalyticity of the OCV at corresponding critical concentrations. It is convenient to replace the lattice gas by an equivalent magnetic system representing a fully occupied lattice of N spins s i which can be either 1 1 (‘up’), or 2 1 (‘down’) according to s i 5 2pi 2 1 [5]. Taking the thermal average we obtain: 1 x 5 ]s1 1 Md 2

(4)

Thus, the fraction of the occupied sites in the original lattice gas is ultimately related to the magnetization per spin M ; ks i l of the corresponding magnetic system. The grand potential becomes:

m qu V 5 2 N ] 2 ] 1 ^Is 2 8

S

D

m (x 5 1 / 2) 5 qu / 2

(7)

The order–disorder transition in the lattice gas has its analogue in the magnetic system. The antiferromagnetic nature of spin–spin interactions implies that for coupling larger than a critical one in weak magnetic fields the structure corresponding to the stable equilibrium is antiferromagnetic with two sublattices: one with spins ‘up’ and the other with spins ‘down’ [6]. In strong enough fields the antiferromagnetic structure becomes unstable and the energetically favourable configuration corresponds to parallel orientation of spins of the both sublattices along the field (paramagnetic state). The phase transition (occurrence and disappearance of the sublattices) takes place when the energy in the external field becomes comparable with the interaction energy. We analyze the lattice gas behavior in terms of the equivalent antiferromagnetic Ising system within the Bethe–Peierls approximation [7,8] which takes into account for each spin its interactions with the nearest neighbors explicitly and with all the rest – by means of averaging. We introduce two elementary q-clusters, each containing a tagged spin s 0 and s˜ 0 in its center belonging to the different sublattices (see Fig. 1). Let us discuss first the cluster shown on the left side of Fig. 1. Its Hamiltonian reads:

(5)

where ^Is is the Helmholtz free energy of an antiferromagnetic Ising model with the coupling u / 4 in an external field: m qu H5]2] (6) 2 4

165

q

b* q 5 K

q

O s s 2 hs 2 h O s 0 j

j 51

0

1

j

j 51

where b 5 1 /(k B T ), K 5 b u / 4, h 5 b H and H1 5

The lattice gas compressibility is proportional to the magnetic susceptibility of the Ising model: 1 dM x 2x 5 ] ] T 4 dH Magnetization is an odd function of the field implying that x(2H ) 5 1 2 x(H ), so that in the lattice–gas terms:

m (x) 1 m (1 2 x) 5 qu for all x providing that for all temperatures:

Fig. 1. q-Clusters for two sublattices. Spins belonging to the ‘ferromagnetic’ sublattice (‘up-spins’) are denoted by open circles. Spins belonging to the ‘antiferromagnetic’ sublattice (‘downspins’) are denoted by crossed circles. External field h is along the positive z axis; directions of ‘interaction fields’ h 1 and h˜ 1 are shown by means of arrows.

V.I. Kalikmanov et al. / Solid State Ionics 127 (2000) 163 – 168

166

h 1 k B T is the (average) field acting on the q neighbors of s 0 . The partition function reads: q

Zq 5 e hse 2K e h 1 1 e K e 2h 1d 1 e 2hse K e h 1 1 e 2K e 2h 1d

F

dg ≠ d ≠ ] 5 2 2g ] 1 ]1 ] dh ≠g dg ≠g1

G

q

The average values of s 0 and s 1 are: ≠ 1 ≠ m 0 ; ks 0 l 5 ] ln Zq , m 1 ; ks 1 l 5 ] ] ln Zq ≠h q ≠h 1

The only term which may cause singularity is dg1 / dg describing the change in h 1 with the external field h (see Fig. 1). Performing routine algebra we find that the singularity condition is: zg 21 2 bg1 1 z 5 0, where b 5 q(1 2 z 2 ) 2 2

Due to symmetry the similar relations hold for the second cluster: ≠ 1 ≠ ks˜ 0 l 5 ] ln Z˜ q , ks˜ 1 l 5 ] ] ln Z˜ q ≠h q ≠h˜ 1

The unknown fields h 1 and h˜ 1 are found from the consistency conditions setting the equivalence between sites of the same sublattice: ks 0 l 5 ks˜ 1 l

(8)

ks 1 l 5 ks˜ 0 l

(9)

Representing the system as a combination of two sublattices we have counted each spin twice, therefore, the total magnetization is M 5 ]12 (m 0 1 m 1 ). Let us define a reduced temperature as:

(12)

Positive roots exist if b . 2z which implies that z , z c 5 (q 2 2) /q. In terms of the coupling constant it means that transition takes place if:

S

1 q K . K Bc 5 ] ln ]] 2 q22

D

Note, that t Bc 5 1 /K Bc coincides with the critical temperature for a ferromagnetic Ising model in a zero field within the Bethe approximation [5]. For t , t cB the solutions of Eq. (12) read: ]]] 1 g1,c 5 ] f b6Œb 2 2 4z 2 g 2z Critical values of the lattice–gas parameters are expressed via g1,c as: z 1 g1,c x c 5 ]]]]] z 1 2g1,c 1 zg 21,c

1 4k B T t ; ] 5 ]] K u and introduce convenient notations: z 5 e 22K , g 5 ˜ e 22h , g1 5 e 22h 1 , g˜ 1 5 e 22h 1 . In the disordered state both sublattices are identical, g1 5 g˜ 1 , and the chemical potential reads:

S

1 1 zg1 m (x) 5 2 k B T ln g1 1 k B T(q 2 1) ln ]] z 1 g1 qu 1] 2

D (10)

where ]]]]]]] 1 g1 (x) 5 ] f 1 2 2x 1œ(1 2 2x)2 1 4z 2 x(1 2 x) g 2xz (11) (note, that at high temperatures m ¯ k B T ln[x /(1 2 x)] independent of the interaction and the type of the lattice). The lattice gas compressibility reads: ~ 0 1m ~ 1 ] where the upper dot dex 2 k B Tx 5 (1 / 8)[m T notes the full derivative with respect to h:

S

D

z 1 g1,c ut ut qu mc 5 2 ] ln g1,c 2 ](q 2 1) ln ]]] 1 ] 4 4 1 1 zg1,c 2 We can obtain the low temperature approximation by expanding g1,c in powers of z 2 ; e 24K . In the 1) leading approximation we have: g (1,c 5 (q 2 2) /z, (2) g 1,c 5 z /(q 2 2) which implies that: 1 1 ) x (c1 ) 5 ], x (2 , z2 < 1 c 512] q q Thus, the critical concentration is finite even at t → 0. This is in agreement with a general argument that since only the nearest neighbor interactions are present, the energy contribution to m becomes important if on average a particle has at least one occupied neighboring site. In Fig. 2 the phase diagram of the square lattice gas (q 5 4) is shown on the x 2 t-plane (upper panel) and m 2 t-plane (lower panel). Areas bounded by the critical curves and the

V.I. Kalikmanov et al. / Solid State Ionics 127 (2000) 163 – 168

167

Fig. 3. OCV for the square lattice and various reduced temperatures. Solid lines: theory. For t , t Bc each curve has two kinks at the critical concentrations; for t . t Bc OCV is analytical. Closed circles, squares and diamonds: Monte Carlo simulations.

Fig. 2. Lattice–gas phase diagram for the square lattice. (a) x 2 t-plane, (b) m 2 t-plane. Critical temperature t Bc ¯ 2.88. Solid lines: theory. Open circles: simulations. Domain between the dashed lines corresponds to a typical range of interaction strength u . 2 4 6k B T.

y-axis correspond to the ordered state. The critical temperature t Bc (q 5 4) ¯ 2.88 (the exact Onsager’s result is: t c ¯ 2.27 [5]). Open circles are results of Monte Carlo simulation of the square lattice gas; phase transition points were estimated from the staggered susceptibility data using the method described in [9] (details of the simulation technique will be presented elsewhere [10]). Domain between the dashed lines corresponds to a typical range of interaction strength u . 2 4 6 k B T. Given a number of nearest neighbors, interaction potential and the absolute temperature we calculate the OCV and its inverse derivative 2 dx / dV for the square lattice (Figs. 3 and 4). If the reduced temperature is higher than the critical, t . t Bc , the lattice gas

Fig. 4. The inverse derivative of the OCV, 2 dx / dV, vs. concentration for the square lattice and various reduced temperatures in the Bethe approximation. For t , t Bc each curve has two singularity points (jumps) at the critical concentrations corresponding to the second order phase transition.

is in the disordered state for all x and m (x) is a smooth function given by Eqs. (10) and (11). For t , t Bc the disordered phase with m (x) given by Eqs. (10) and (11) exists only at the small and large concentrations, i.e. for x , x c( 1 ) and x . x c(2 ) . In the intermediate range x (c1 ) , x , x (c2 ) corresponding to m c(1) , m , m c( 2) (see Fig. 2) the lattice gas is in the ordered state, for which instead of searching for m (x) we calculate the inverse function x 5 x( m ) by gradu-

168

V.I. Kalikmanov et al. / Solid State Ionics 127 (2000) 163 – 168

ally varying the chemical potential. The set of consistency Eqs. (8) and (9) is solved iteratively. To find g1 and g˜ 1 for an arbitrary point (t, m ) of the parameter space we start from the known solution at (t 0 5 0, m0 5 qu / 2) where all sites of the first sublattice are occupied while all sites of the second sublattice are free (recall that m0 5 qu / 2 corresponds to H 5 0 according to Eq. (6)). This implies that g1 (t 0 , m0 ) 5 1 `, g˜ 1 (t 0 , m0 ) 5 0. The temperature interval (t 0 ,t) is divided into small steps Dt and the solution for the point (t, m0 ) is found by iterations in Dt starting from (t 0 , m0 ). When this is done the chemical potential interval ( m0 , m ) is divided into small steps Dm and iterations are performed in Dm at the given t starting from the point (t, m0 ). At temperatures below the critical, t , t Bc , each OCV curve in Fig. 3 has two singularity points (kinks) in which the second order phase transition takes place. In the vicinity of x 5 1 / 2 OCV rapidly decreases though remaining analytical. The lower the temperature the more pronounced is the voltage drop; at t 5 0 the OCV becomes a step function at x 5 1 / 2. It is important to stress that the step-like behavior at x 5 1 / 2 at non-zero temperatures is a consequence of the phase transitions taking place at x (c1 ) and x (c2 ) but not at the point x 5 1 / 2 itself. For t . t Bc OCV is analytical for all concentrations. Monte Carlo results are in good agreement with the theory. Calculations of the inverse derivative of the OCV, 2 dx / dV, are presented in Fig. 4. At high temperatures 2 dx / dV is an analytical function of x with a single maximum at x 5 1 / 2. For lower t, but still beyond t Bc , it remains analytical but the central maximum becomes minimum while two other maxima appear, which are symmetric about x 5 1 / 2.

Below t Bc the function retains both maxima becoming discontinuous at concentrations corresponding to the phase transition. As a concluding remark note that by its nature the lattice–gas model discussed in the present paper depends on a single structural parameter: the coordination number of the interstitial lattice q, originating from the structural properties of the host matrix. However, peculiarities of the host matrix manifest themselves also in the effective interaction potential u. Besides, one has to bear in mind that the results obtained refer exclusively to the bichromatic lattices. To illustrate this, mention as an example that the model is valid for the simple cubic lattice with q 5 6 but is not applicable for the triangular lattice with the same q 5 6, since the latter is frustrated.

References [1] A.D. Robertson, A.R. West, A.G. Ritchie, Solid State Ionics 104 (1997) 1. [2] M.K. Aydinol, A.F. Kohan, G. Ceder, K. Cho, J. Joannopoulus, Phys. Rev. B 56 (1997) 1354. [3] A. Van der Ven, M.K. Aydinol, G. Ceder, J. Electrochem. Soc. 145 (1998) 2149. [4] L.D. Landau, E.M. Lifshitz, Statistical Physics, Pergamon Press, Oxford, 1969. [5] R.J. Baxter, Exactly Solved Models in Statistical Mechanics, Academic Press, London, 1982. [6] L.D. Landau, E.M. Lifshitz, Electrodynamics of Continuous Media, Pergamon Press, Oxford, 1984. [7] H. Bethe, Proc. R. Soc. 150 (1935) 552. [8] R. Peierls, Proc. R. Soc. 154 (1936) 207. ¨ Physica A 251 (1998) 224. [9] J.R. Heringa, H.W.J. Blote, [10] M.V. Koudriachova et al. (manuscript in preparation).