Numerical inversion of Laplace transforms.

Journal of

Electroanalytical Chemistry Journal of Electroanalytical Chemistry 608 (2007) 37–46 www.elsevier.com/locate/jelechem

Numerical inversion of Laplace transforms. A useful tool for evaluation of chemical diffusion coefficients in ion-insertion electrodes investigated by PITT C. Montella *, R. Michel, J.P. Diard E´cole Nationale Supe´rieure d’E´lectrochimie et d’E´lectrome´tallurgie de Grenoble, Laboratoire d’E´lectrochimie et de Physicochimie des Mate´riaux et Interfaces, UMR 5631 CNRS-INPG-UJF, Domaine Universitaire, B.P. 75, 38402 Saint Martin d’He`res, France Received 2 February 2007; received in revised form 25 April 2007; accepted 8 May 2007 Available online 24 May 2007

Abstract This work starts with PITT experimental data reported by Xia et al. from thin-film LiCoO2 electrodes prepared on a Si substrate by pulse laser deposition. In this article, we show that (i) the value of chemical diffusion coefficient of lithium-ion evaluated by graphical processing techniques is often erroneous, (ii) non-linear curve fitting is a better approach for accurate determination of thermodynamic, kinetic and diffusion parameters of ion-insertion electrodes from PITT data, (iii) the concept of anomalous diffusion can be applied to ion-insertion electrodes investigated by PITT, (iv) numerical inversion of Laplace transforms is a useful tool to derive the relevant model for representation of PITT data and (v) implementation of the whole calculation procedure becomes an easy task in the environment of Mathematica or other equivalent softwares.  2007 Published by Elsevier B.V. Keywords: PITT; Diffusion coefficient; Intercalation; Numerical inversion of Laplace transforms; Non-linear fit

1. Introduction Since the pioneering works of Weppner and Huggins [1] and Wen et al. [2], the so-called potentiostatic intermittent titration technique (PITT) has been widely employed for investigation of ion-insertion electrodes (IIEs) involved in electrochemical batteries and electrochromic devices. A large bibliography can be found in the recent articles by Deiss [3] and Markevich et al. [4]. PITT consists in the application of a staircase potential signal to IIEs and measurement of the resulting current transients. At the nth step of the input signal the electrode is perturbed by imposing a potential step from the initial potential En1 to the applied potential En = En1 ± DEn, where n = 1,2, . . . is the running index. The current transient In(t), due to potential stepping, is measured until the *

Corresponding author. Tel.: +33 4 76826526; fax: +33 4 76826630. E-mail address: [email protected] (C. Montella).

0022-0728/$ - see front matter  2007 Published by Elsevier B.V. doi:10.1016/j.jelechem.2007.05.011

current approaches zero because of restricted (finite-space) diffusion conditions in the host material, and hence full equilibrium conditions are reached at the applied potential. The incremental charge passed after the application of potential perturbation is calculated R 1 by integration of current with respect to time, DQn ¼ 0 I n ðsÞds; and the differential intercalation capacitance of host material is obtained, from the incremental charge and the potential step amplitude, as Cint,n = DQn/DEn. Thermodynamic informations about insertion host materials are directly accessible from the experimental dependence, Cint versus electrode potential, while diffusion and kinetic parameters of IIEs are generally determined by means of theoretical analysis of current transients through graphical [1,2,5–7] and/or numerical [8,9] processing techniques based upon the analytical solution of Fick’s equations of diffusion under the relevant initial and boundary conditions in the host material. The main aims of this work are the following: (i) to point out some problems relative to graphical processing

38

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

Nomenclature Abbreviations ADI anomalous diffusion impedance GS Gaver–Stehfest algorithm working in standard (double) precision GSmp GS algorithm working in multi-precision computing environment IIE ion-insertion electrode NILT numerical inversion of Laplace transforms NL-Fit non-linear fit procedure PITT potentiostatic intermittent titration technique SEI solid electrolyte interphase Roman A c Cint D E F I(t) I(0) I0 Ik J

letters electroactive surface area concentration differential intercalation capacitance chemical diffusion coefficient electrode potential Faraday’s constant Faradaic current initial current level I0 = I(0) kth value of measured current diffusion flux

of PITT experimental data when the influence of heterogeneous kinetics and Ohmic drop is disregarded, (ii) to promote the use of non-linear fit (NL-Fit) procedures for accurate determination of thermodynamic, kinetic and diffusion parameters of IIEs from PITT data, (iii) to propose a set of four graphical representations of experimental data and theoretical curves in order to judge easily the goodness of fit, (iv) to apply the concept of anomalous diffusion to PITT data, (v) to prove the usefulness of numerical inversion of Laplace transforms (NILT) to derive the relevant model employed in the NL-Fit procedure, and finally (vi) to show that implementation of both NILT and NL-Fit techniques becomes an easy task in a high-level language for scientific computing like MathematicaTM [10] or other equivalent softwares. 2. Problems relative to graphical processing of PITT data Both theoretical and experimental problems are related to the application of PITT to insertion electrodes. Such problems were thoroughly examined in the recent papers by Deiss [3], Markevich et al. [4], Han et al. [11] and Montella [12]. In this article we are mainly concerned with the substantial error in the determination of chemical diffusion coefficients resulting from the approximation of diffusioncontrolled processes. The problem is well illustrated by taking the example of lithium-ion insertion in LiCoO2 film electrodes. Among the

L M*(s) nG Rd Rext s t Vm x x, y z Z(s) Zd(s)

diffusion length dimensionless mass-transfer function number of Gaver functionals in GS algorithm diffusion resistance sum of external resistances Laplace variable time molar volume Abscissa insertion level charge of mobile ion electrode impedance diffusion impedance

Greek letters bk kth positive root of Eq. (5) c fractional order of derivation DE amplitude of potential step DQ incremental charge erel relative residual K dimensionless parameter in Eq. (6) sd diffusion time constant xk kth weight factor in the NL-Fit procedure

large number of references pertaining to this oxide, we selected the article by Xia et al. [13] because of (i) the use of a pure (binder-free) and compact (without any cracks and pinholes) material prepared on a Si substrate by pulse laser deposition, (ii) the experimental data collected in the composition range where LixCoO2 exists as a single phase, (iii) linear diffusion of lithium-ion in a thin-film planar electrode, and finally (iv) the use of small potential steps (DEn = ±10 mV). Assuming very fast heterogeneous kinetics and negligible Ohmic drop, the potentiostatic current transient is described by the well-known series [1,2]:   1 2 DQ X 2 p t IðtÞ ¼ 2 exp ð2k  1Þ ð1Þ sd k¼1 4sd where the running index ‘n = 1, 2, . . . ’ is omitted in In (t), DQn and sd,n for the sake of simplification of notation. The same applies in what follows. In Eq. (1), sd = L2/D is the time constant for linear diffusion, L is the film thickness and D is the chemical diffusion coefficient of guest ion in the host material. An exponential decay of diffusion current with respect to time is predicted from the leading term of the above series within the long-time approximation (t  sd). Then, the diffusion coefficient can be calculated from the slope of the asymptotic straight line in the plot of lnjI(t)j versus time: D¼

4L2 d ln jIðtÞj t; d p2

ðt  sd Þ

ð2Þ

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

The above equation was used by Cao and Prakash [14] to evaluate the chemical diffusion coefficient of lithiumion in LiMn2 O4 spinel thin-film electrodes. The same equation was employed by Das et al. [15] for the analysis of nano-crystalline lithium manganate thin-films. A similar approach was used by Rho et al. to evaluate the chemical diffusion coefficient of lithium-ion in LiCoO2 [16], Li4Ti5O12 [17] and LiMn2O4 [18] film electrodes prepared by spin-coating and sol–gel methods. Hydrogen absorption and diffusion in vanadium [19] and gadolinium [20] thinfilms was investigated by Di Vece et al. through the same equation. On the other hand, both slope and intercept of the asymptotic straight line were used by Levi et al. [21] for interpretation of the chronoamperometric responses of thin-film lithium-insertion V2O5 electrodes. Finally, Xia et al. [13] employed Eq. (2) to investigate the potential dependence of the chemical diffusion coefficient of lithiumion in LiCoO2. By way of illustration, the values L = 3 · 105 cm and D = 3 · 1012 cm2 s1 were reported by these authors when the electrode potential is stepped from 3.98 to 3.97 V vs. Li/Li+. From the above values, we get the diffusion time-constant, sd = L2/D = 300 s. Experimental data scanned and discretized from Fig. 10 in reference [13] are presented by dots in Fig. 1 of this article. The interpolated curve has been plotted using cubic spline interpolation provided by the ‘Interpolation’ command of Mathematica. The incremental charge is calculated by integration of current with respect to time, using the ‘Integrate’ command of the same software, as DQ = 7.9 · 104 C. Then, the theoretical variation of current versus time can be plotted from Eq. (1), setting sd = 300 s and DQ = 7.9 · 104 C. The current transient is compared to the experimental data of Xia et al. in Fig. 2 of this article. A large discrepancy in value is observed between the experimental and theoretical data, thus indicating that the value of chemical diffusion coefficient evaluated from Eq. (2) is erroneous. Indeed, due account of both diffusional processes in the electrode bulk and kinetic limitations on its surface is an important fea-

39

Fig. 2. Potentiostatic current transient of LixCoO2 film electrode caused by a potential step from 3.98 V to 3.97 V vs. Li/Li+(x  0.6), applied at t = 0. The experimental data (dots) of Xia et al. [13] are compared to the theoretical prediction (solid line) from Eq. (1).

ture of electrochemical insertion/extraction into/from host materials [8,22–24]. 3. Theoretical expression of potentiostatic current transient The theory of PITT is now well established [8,22–24], at least for thin planar films of pure (binder-free), singlephase and ideal materials (without structural defects). The model covers the cases where finite-space diffusion in the electrode, slow heterogeneous kinetics on its surface, Ohmic drop effects and ionic transfer/transport through the passive surface layer (so-called solid electrolyte interphase (SEI)) are taken into consideration. Provided phase transformation is disregarded (singlephase materials and/or potential domains) and the amplitude of potential step is small enough that linearity conditions are experimentally satisfied, the electrode response is given by [22]: IðtÞ ¼ DEL1

1 sZðsÞ

ð3Þ

where s is the Laplace complex variable, Z is the IIE impedance measured at the same initial potential and L1 denotes the inverse Laplace transformation. Provided the double-layer, passive film and contact capacitances are negligible, compared to Cint, the IIE impedance takes on the expression reported in the second entry to Table 1. Then, the current transient derived from Eq. (3) is described by the infinite series [22,23]:   1 X K 2 t IðtÞ ¼ 2Ið0Þ exp b ð4Þ k 2 2 sd k¼1 K þ K þ bk

Fig. 1. Plot of current versus time. The experimental data (dots) have been scanned and discretized from Fig. 10 of reference [13]. Noisy data measured at t > 400 s have been removed. The solid line is the interpolated curve plotted using Mathematica.

where the initial current level is I(0) = DE/Rext = KDQ/sd, the coefficient bk is the kth positive root of the transcendent equation: b tan b  K ¼ 0 and the dimensionless parameter:

ð5Þ

40

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

Table 1 Theoretical expressions of the dimensionless mass-transfer function M*(s) and the electrode impedance Z(s) pertaining to IIEs M*(s) pffiffiffiffiffi Coth ssd pffiffiffiffiffi ss d

As above Coth½ðssd Þc=2  ðssd Þc=2 Coth½ðssd Þc=2  ðssd Þ1c=2

K

Z(s)

I(t)

Refs.

1

RdM*(s)

Eq. (1)

[1,2]

Finite

Rext[1 + KM *(s) ]

Eq. (4)

[22,23]

As above

As above

–

[30,31]

As above

As above

Eq. (13)

[31]

Restricted (finite-space) diffusion of guest ions in a thin planar film of pure (binder-free) insertion host material is assumed. The double-layer, surface layer (SEI) and contact capacitances are disregarded in Z(s).

K¼

Rd Rext

ð6Þ

is the quotient of the diffusion resistance of guest ions by the sum of external resistances, i.e. the sum of Ohmic resistances, charge-transfer resistance, contact resistances, etc. It is noticed here that Eq. (4) does not apply at very short times where the Faradaic current is corrupted by the double-layer, passive film and contact capacitances. An exponential decay of Faradaic current with respect to time is predicted from Eq. (4) within the long-time approximation (t  sd), setting I(0) = KDQ/sd [23], as:   DQ K2 2 t exp b1 ð7Þ IðtÞ ¼ 2 ; ðt  sd Þ sd K2 þ K þ b21 sd where the first positive root, b1, of Eq. (5) can be well pffiffiffiffi approximated by a rational function with respect to K [12]. Both slope (Sl) and intercept (Int) of the asymptotic straight line in the plot of lnjI(t)j versus time are then required for evaluation of K, sd and therefore the chemical diffusion coefficient of guest ions from PITT data [9,23]: ! b21 jDQj K2 Sl ¼  ; Int ¼ ln 2 ð8Þ sd K2 þ K þ b21 sd As indicated in a previous work [23], the chemical diffusion coefficient of guest ions is underestimated when Eq. (2) is used in place of Eq. (8). In particular, a substantial error in the evaluation of D is predicted at small values of K. In addition, Eqs. (7) and (8) are only valid in the long-time region (t  sd). Alternatively, the more general Eqs. (4) and (5) can be used to fit PITT experimental data in the whole time domain of interest. 4. Non-linear curve fitting NL-Fit procedures, based on Eqs. (4) and (5), were proposed in previous works focused on hydrogen absorption/ desorption into/from Pd and Pd-alloy foils [25], as well as lithium-ion insertion/extraction into/from thin-films of carbon and transition metal oxides [8,9,23]. The series in Eq. (4) can be limited to its first terms provided K 6 10 [9], and the bk,s calculated from Eq. (5) by any numerical

method for solving transcendent equations, e.g. the ‘FindRoot’ procedure of Mathematica. The ‘NonlinearRegress’ procedure of the same software is applied to the experimental data of Xia et al. [13] in Fig. 3 of this article. Due to monotonous decrease of jI(t)j with respect to time during PITT experiments, under restricted (finite-space) diffusion conditions in insertion electrodes, the influence of long-time data is very low when the current is directly used in the NL-Fit procedure. The usual approach is then to use weight factors (xk) in the fitting procedure in order to minimize the sum of square of weighted residuals. Setting xk = 1/Ik, where Ik is the kth value of measured current after the application of potential perturbation, the influence of the experimental data measured in the long-time region is increased. Alternatively, logjI(t)j could be used in the NL-Fit procedure, rather than I(t). As a matter of fact, the best-fitting parameters obtained by the two methods are nearly the same. Hence, only the first procedure has been used in Fig. 3. Four graphical representations of the experimental data (dots) and the best-fitting curves (solid lines) are presented in Fig. 3, in order to judge easily the goodness of fit. The current transient is plotted first in Fig. 3a, whilepthe ffi Cottrell function representation of PITT data, IðtÞ t versus log t, initially proposed by Levi and Aurbach [26] and widely used in the electrochemical literature focused on IIEs, is shown in Fig. 3b. The semi-logarithmic plot in Fig. 3c is a usual representation of PITT data since the pioneering works of Weppner and Huggins [1] and Wen et al. [2]. Finally, the logarithmic current is plotted versus logarithmic time in Fig. 3d. This representation is essentially employed in the articles by Shin and Pyun [27]. The four graphical representations, used together in Fig. 3, can be proposed as a standard representation for easy (visual) comparison of PITT experimental and theoretical data. By visual inspection of Fig. 3, and more especially of Fig. 3b, it is observed that the series in Eq. (4) is only a first (rough) approximation of the experimental data collected from LixCoO2 film electrodes under the operating conditions of reference [13]. 5. Anomalous diffusion in IIEs Close examination of the electrochemical impedance spectra (EIS) plotted by Xia et al. in Fig. 7 of reference [13] shows that the behavior of LixCoO2 electrodes is not purely capacitive in the low-frequency (LF) range, in contrast to the theoretical prediction from the second entry to Table 1. Indeed, a steep sloping line is observed in the LF domain, rather than a vertical straight line. In a previous paper [28] some of us derived the general formulation of linear diffusion impedance as: Z d ðsÞ ¼ Rd M  ðsÞ M*(s)

ð9Þ

is the relevant expression of dimensionless where mass-transfer function (DMTF) pertaining to the electrochemical system under consideration. A lot of limiting

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

41

I ( t ) t μA s1 2

I ( t)

μA

9

6

3

0 0

200

30

15

0 0

400

1

3

2

3

log (t s )

t s -5

lo g I ( t) A

-5

lo g I ( t) A

2

-6

-6

-7

-7 0

200

400

0

1

log (t s )

t s

Fig. 3. The experimental data (dots), from Figs. 1 and 2, are compared to the best-fitting curves (solid lines) obtained from Eq. (4), using the four graphical representations: (a) current versus time, (b) Cottrell function plot, (c) semi-logarithmic plot and (d) log–log plot. The ‘NonlinearRegress’ procedure of Mathematica has been employed together with the Levenberg–Marquardt algorithm for minimization of the sum of square of weighted residuals (xk = 1/Ik). The best-fitting parameters are I0 = 11.8 lA, K = 1.7 and sd = 114 s.

expressions of DMTFs were examined in the above article. From Eq. (9), the DMTF is equal to the quotient Zd(s)/Rd. In terms of equivalent circuit [29], this impedance can be represented by a uniformly distributed structure terminated by an appropriately chosen impedance which is related to the DMTF under consideration. From Eqs. (3), (6) and (9), the theoretical response of IIEs to a small potential step is given by the following equation where I(0) = DE/Rext: IðtÞ ¼ Ið0ÞL1 fs½1 þ KM  ðsÞg1

ð10Þ

For example, Eq. (4) can be obtained directly from Eq. (10) and the second entry to Table 1 [22,23]. Unfortunately, Eq. (4) fails to predict accurately the potentiostatic response of LixCoO2 film electrode under the operating conditions of reference [13]. Hence, a description of transport mechanisms which differ from ordinary diffusion is required. We suggested the DMTF reported in the third entry to Table 1, on a heuristic basis, for investigation of hydrogen absorption/desorption processes into/from HxNb2O5 thinfilm electrodes by EIS [30]. However, Bisquert and Compte [31] showed that the main problem relative to its application to diffusion in solid materials is that the physical system is then in a situation where the number of diffusing

particles is not conserved. Hence, the above DMTF is disregarded in what follows. Using a theoretical approach to anomalous diffusion impedance (ADI) that employs fractional calculus, several models were presented by Bisquert and Compte [31] depending on the microscopic picture of diffusion processes in solid materials. Anomalous diffusion is characterized by a mean squared displacement of the diffusing particles that does not follow the ordinary linear law Ær2æ  t, but, more generally, has a power law dependence on time [31]. The DMTF reported in the fourth entry to Table 1 comes from the substitution of the generalized constitutive equation of diffusion flux [31]: J ¼ D

o1c oc ; ot1c ox

c61

ð11Þ

into the usual mass balance equation, which leads to the generalized diffusion equation [31]: oc c o2 c ¼ D ; otc ox2

c61

ð12Þ

where the diffusion constant, D/(cm2 sc), has no more conventional units, c is the fractional order of derivation operator and the other symbols have their usual meaning

42

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

and units. Eqs. (11) and (12) can be derived as a limit from a stochastic scheme (continuous time random walk) in the article by Compte and Metzler [32]. The time domain model used hereafter for representation of PITT data is derived, from Eq. (10) and the fourth entry to Table 1, as: h i391 82 c=2 < = Coth ðssd Þ 5 IðtÞ ¼ Ið0ÞL1 s41 þ K ; c61 1c=2 : ; ðssd Þ ð13Þ where I(0) = DE/Rext and K = Rd/Rext, with: Rd ¼

L2=c1 zFAD1=c jdc=dEj

ð14Þ

and sd ¼ ðL2 =DÞ

1=c

ð15Þ

The above formulations of diffusion resistance and diffusion time constant keep the dimensions in order, Rd/X and sd/s, with D/(cm2 sc). In addition, z is the charge of mobile ion (e.g. z = 1 for lithium-ion insertion/extraction processes), F is the Faraday’s constant, A is the electroactive surface area between the electrolyte and the host material, and dc/dE is the insertion-isotherm slope, i.e. the reciprocal of slope of the coulometric titration curve of host material, dc=dE ¼ V 1 m ðdy=dEÞ, where Vm and y denote, respectively, the molar volume and the insertion level. Setting c = 1 in Eqs. (13)–(15), the usual formulations apply and the ideal behavior of IIEs is recovered from Eq. (4). In contrast, the inverse Laplace transform is unknown when c < 1. Hence, there is a lack of mathematical model for representation of PITT experimental data in this case. Fortunately, efficient algorithms are available from the specialized literature for numerical inversion of Laplace transforms. 6. Usefulness of numerical inversion of Laplace transforms It is surprising how little attention this approach has received in the literature relative to electrochemical kinetics, although a large number of engineering application papers have been published in the last decade. A bibliography of several hundred such papers is available on the WEB [33]. As we are concerned here with electrochemical systems, the following references should be quoted. Floriano et al. [34] employed the Crump method [35] for simulating and fitting the oxygen reduction reaction at a Nafion-covered rotating-disk electrode investigated by the chronoamperometry technique. Ogata et al. [36] used the so-called fast numerical inversion of Laplace transformation, where exp(st) in the Bromwich inversion integral is replaced by an approximating function, for solving the diffusion equations relative to hydrogen permeation in a multilayered membrane. Later, Barsoukov et al. [37] adapted the above method to prediction of the time domain charge

curves resulting from lithium-ion intercalation into mesocarbon–microbeads based anodes. Finally, the Stehfest method [38] was employed in the recent and attractive textbook by Honeychurch [39] to plot concentration profiles of electroactive soluble species, caused by potential stepping. In this article we use the Stehfest method, also called Gaver–Stehfest (GS) algorithm [40]. A Mathematica package called ‘NumericalInversion’ is available from the MathSource library on the WEB [41]. This package provides five methods (the Stehfest method being included) for numerical inversion of Laplace transforms (NILT), along with references to the original mathematical literature. Following Abate and Valko´ [40], the GS algorithm contains only one parameter nG that is the number of Gaver functionals (approximants) used in the calculation, while n = 2nG is the number of terms used in the Salzer summation to accelerate convergence [42]. It should be noted that the GS algorithm results from a finite linear combination of values of the Laplace transform to be inverted. Its implementation requires only two lines programming when the inverse Laplace transform is evaluated with standard (double) precision computing. As soon as an efficient approximation of the IIE response to a small potential step is obtained by NILT from Eq. (13), it can be used as a model for representation of PITT experimental data. Then, the ‘NonlinearRegress’ procedure of Mathematica is well suited to estimation of best-fitting parameters from the above model. 7. Using the whole ADI–NILT–NL-fit procedure The whole procedure (GS approximation of Eq. (13) and NL-Fit of experimental data) requires only a few lines programming in Mathematica or equivalent scientific languages. It is applied to the experimental data of Xia et al. [13] in Fig. 4 of this article. A good fit of PITT data is observed, whatever the graphical representation is employed in the figure. The estimates of parameters, I0, K, sd and c, are reported in Table 2, together with the 95% confidence intervals. In addition, it is observed in Fig. 5 that the relative deviation of experimental data from the proposed model is nearly randomly distributed around zero, thus validating the statistical analysis. Indeed, the mean value of relative residual is Æerelæ = 2.8 · 104, while the mean value of jerelj is equal to 7.8 · 103. From the numerical value of diffusion time constant in Table 2, and the film thickness given by Xia et al. [13], L = 3 · 105 cm, the diffusion constant of lithium-ion in LixCoO2 is evaluated as D ¼ L2 =scd ¼ 8:7  1 011 cm2 sc at 3.98 V vs. Li/Li+ where LixCoO2 (x  0.6) exists as a rhombohedral single-phase [43,44]. On the other hand, an apparent value of chemical diffusion coefficient, with usual units, can be calculated as Dap = L2/sd = 2.2 · 1011 cm2 s1, which is approximately one order of magnitude larger than the value calculated by Xia et al. from Eq. (2) above.

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

43

Fig. 4. The experimental data (dots), from Figs. 1–3, are compared to the best-fitting curves (solid lines) obtained from the GS approximation (nG = 8) of Eq. (13). The ‘NonlinearRegress’ procedure of Mathematica has been employed together with the Levenberg–Marquardt algorithm for minimization of the sum of square of weighted residuals (xk = 1/Ik). The best-fitting parameters are reported in Table 2.

Table 2 Estimates of parameters and 95% confidence intervals (CI) obtained from the experimental data of Xia et al. [13] and the PITT model derived by GS approximation (nG = 8) of Eq. (13) Parameter

Estimate

95% CI

I0/lA K sd/s c

9.33 0.45 41.0 0.63

9.17–9.48 0.38–0.53 35.2–46.8 0.57–0.69

The ‘NonlinearRegress’ procedure of Mathematica has been employed together with the Levenberg–Marquardt algorithm for minimization of the sum of square of weighted residuals (xk = 1/Ik).

Moreover, from the best-fitting parameters in Table 2, we get the values Rext = DE/I0  1072 X and Rd = KRext  486 X. The differential intercalation capacitance is obtained as Cint = sd/Rd = 8.4 · 102 F. It can be compared to the value Cint = DQ/DE = 7.9 · 102 F where the incremental charge is calculated by integration of current with respect to time as DQ = 7.9 · 104 C, as well as to the incremental capacity measured by Xia et al. from the charge–discharge curves of LixCoO2 electrode. Cint  0.18 F is obtained at 3.98 V vs. Li/Li+ from Fig. 4 in reference [13]. Clearly this value is overestimated by a factor  2, probably because

Fig. 5. Relative deviation of the experimental data of Xia et al. [13] from the proposed model, plotted versus time.

the charge–discharge curves were plotted by these authors far from equilibrium conditions. The slope of the coulometric titration curve of LixCoO2 could also be estimated from the formulation Cint = zFALjdc/dEj, but, unfortunately, the electroactive surface area is not indicated in the article of Xia et al. However, the values Ajdc/dEj = 2.9 · 102 mol cm1 V1 and Ajdx/dEj = 0.57 cm2 V1 are obtained from Cint, using Vm = 19.56 cm3 mol1 [13].

44

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

Hence, accurate values of thermodynamic, kinetic and diffusion parameters relative to lithium-ion insertion in LixCoO2 can be obtained from the proposed model by using the GS approximation of Eq. (13), coupled to a NL-Fit procedure. 8. Testing the accuracy of numerical inversion with mathematica The main advantage of GS algorithm is the simple formulation of the Gaver–Stehfest inversion formula. In addition, the above algorithm is directly compatible with the ‘NonlinearRegress’ procedure of Mathematica. Indeed, the statistical analysis of numerical results (estimates of parameters, confidence intervals, estimated variance, asymptotic correlation matrix, Fit curvature Table) are directly accessible from the NL-Fit procedure. However, because of the presence of the binomial coefficients in the weights, the GS formula tends to require high system precision in order to yield good accuracy in the calculations [45]. In other words, the process is numerically unstable in a fixed-precision computing environment because of the propagation of round-off errors when the number of Gaver functionals is increased in the calculation. A significant improvement of NILT methods was suggested recently by Abate and Valko´ [40]. Their key idea was to use multi-precision computing in the environment of Mathematica. Following the suggestion of these authors, a home-made version of GS algorithm, called GSmp below, was written in Mathematica language. As recommended by Abate and Valko´, the precision was set at 2.2 nG, where nG is the number of Gaver functionals used in the calculations. We tested the performance of GSmp algorithm for the computation of ‘good’ transforms. Transforms are said to be ‘good’ by Abate and Valko´ if the transforms have all their singularities on the real axis to the left of s = a, with a P 0, and the functions are infinitely differentiable for all t > 0. We observed that about 0.9 nG significant digits are produced by the GSmp algorithm with such transforms, in agreement with the predictions by Abate and Valko´ [40]. Hence, the home-made GSmp algorithm performs well. Note however that the Laplace transforms reported in Table 1 are not ‘good’, because they have singularities out of the real axis. Then, the number of significant digits of the numerical inverse transforms may be not proportional to nG [40]. Moreover, we observed that the number of significant digits decreases at constant nG and rising time. The GSmp algorithm cannot be employed directly in the NL-Fit procedure because of the approximate numbers (with decimal point) used to estimate the model parameters. However, GSmp can be used to control the accuracy of GS algorithm and therefore to validate, or not, the results of NL-Fits. This approach is illustrated in Fig. 6 where the relative deviation of PITT experimental data from the proposed model (dots) is compared to the relative

accuracy of GS approximation (solid line), which is the relative deviation of GS from GSmp output. nG = 8 and nG = 20 have been used in Fig. 6 for the GS and GSmp algorithms, respectively. In addition, the best-fitting parameters, as well as the list of measurement times have been rationalized before computing in order to yield the required precision with the GSmp algorithm. A logarithmic scale is used on Y-axis in Fig. 6. It is observed that the accuracy of GS algorithm is sufficient

Fig. 6. The relative residual jerelj (dots) is plotted versus time, from Fig. 5, and compared to the relative accuracy of GS approximation (solid line), using a logarithmic scale. nG = 8 and nG = 20 have been used in the GS and GSmp procedures, respectively.

Fig. 7. Same figure captions than in Fig. 6, except for the values, nG = 6 (a) and 10 (b), used in the GS procedure.

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

compared to the residuals of the NL-Fit procedure. Indeed, the residuals are two or more orders of magnitude larger than the GS accuracy, depending on the time value. Hence, the standard GS method performs well under the operating conditions in this work. Note that n = 2nG = 16 terms in the Gaver–Stehfest inversion formula is the optimal number of terms that is recommended in the mathematical literature. The problem related to the value of nG is illustrated by the examples presented in Figs. 6 and 7. As a matter of fact, the performance (accuracy) of GS algorithm increases with nG, up to the optimal value nG = 8 or 9, and then it decreases rapidly at higher values of nG. The ‘noise’ generated by the GS algorithm in Figs. 6 and 7b is probably due to the propagation of round-off errors, while the GSmp algorithm remains stable irrespective of the value of nG. 9. Conclusion Starting with a set of PITT experimental data from a pure (binder-free), compact (crack-free) and single-phase LixCoO2 film electrode prepared and investigated by Xia et al. [13], we have shown in this article that:  the numerical value of chemical diffusion coefficient of lithium ion calculated from Eq. (2) above is erroneous because of the influence of heterogeneous kinetics and Ohmic drop effects,  a non-linear curve fitting procedure is required to compare PITT experimental and theoretical data in the whole time domain of interest, and hence to estimate correctly the values of kinetic and diffusion parameters from ion-insertion electrodes,  the set of four graphical representations of PITT data, used together in Figs. 3 and 4 of this article, is well suited to judge visually the goodness of fit,  the theoretical expression of potentiostatic current transient derived in previous works [22–24] is only a first (rough) approximation of PITT experimental data under the operating conditions in reference [13],  a more efficient model for representation of PITT data is available from Eq. (13) above, using the concept of anomalous diffusion [31] in IIEs, at least for thin-film materials investigated in single-phase potential domains,  numerical inversion of Laplace transforms is required to derive the relevant model in the time domain from Eq. (13),  the Gaver–Stehfest (GS) algorithm, which is available from the specialized literature, performs well under the operating conditions in this work. In addition, it can be easily coupled to standard NL-Fit procedures through a few lines programming,  the performance of GS algorithm can be controlled, a posteriori, by comparison to the home-made GSmp algorithm working in the multi-precision computing environment of Mathematica,

45

 finally, the proposed model in Eq. (13) fairly describes the experimental data of Xia et al. [13]. Using the same procedure, we tried successfully to fit other PITT data from the electrochemical literature focused on IIEs, although the obtained results are not presented here for shortening reasons. Such results will be analyzed elsewhere [46]. The new approach proposed and experimented in this work for estimation of thermodynamic, kinetic and diffusion parameters of ion-insertion electrodes from PITT experimental data can be easily adapted to other intermittent titration techniques (e.g. GITT), other IIEs (electrode or particle geometry, diffusion length distribution, diffusion trapping processes, multilayered materials and composite electrodes could be investigated in the time domain through NILTs provided the relevant expression of IIE impedance is known a priori), as well as to other electrochemical systems/ methods. In addition, different approaches can be envisaged to fit PITT and GITT experimental data, because different models are available to fit the impedance spectra of IIEs. Discrimination of such models in the time domain will be greatly facilitated by using numerical inversion of Laplace transforms. All these points will be developed in forthcoming articles. References [1] W. Weppner, R.A. Huggins, Ann. Rev. Mater. Sci. 8 (1978) 269. [2] C.J. Wen, B.A. Boukamp, R.A. Huggins, W. Weppner, J. Electrochem. Soc. 126 (1979) 2258. [3] E. Deiss, Electrochim. Acta 47 (2002) 4027. [4] E. Markevich, M.D. Levi, D. Aurbach, J. Electroanal. Chem. 580 (2005) 231. [5] M.A. Vorotyntsev, M.D. Levi, D. Aurbach, J. Electroanal. Chem. 572 (2004) 299. [6] M.D. Levi, E. Markevich, D. Aurbach, J. Phys. Chem. B 109 (2005) 7420. [7] M.D. Levi, R. Demadrille, A. Pron, M.A. Vorotyntsev, Y. Gofer, D. Aurbach, J. Electrochem. Soc. 152 (2005) E61. [8] A.V. Churikov, A.V. Ivanischev, Electrochim. Acta 48 (2003) 3677. [9] C. Montella, J. Electroanal. Chem., submitted for publication. [10] S. Wolfram, The Mathematica Book, Wolfram Media/Cambridge University Press, Cambridge, 1999, Version 4. [11] B.C. Han, A. Van der Ven, D. Morgan, G. Ceder, Electrochim. Acta 49 (2004) 4691. [12] C. Montella, Electrochim. Acta 51 (2006) 3102. [13] Hui Xia, Li Lu, G. Ceder, J. Power Sources 159 (2006) 1422. [14] F. Cao, J. Prakash, Electrochim. Acta 47 (2002) 1607. [15] S.R. Das, S.B. Majumder, R.S. Katiyar, J. Power Sources 139 (2005) 261. [16] Y.H. Rho, K. Kanamura, J. Electrochem. Soc. 151 (2004) A1406. [17] Y.H. Rho, K. Kanamura, J. Solid State Chem. 177 (2004) 2094. [18] Y.H. Rho, K. Dokko, K. Kanamura, J. Power Sources 157 (2006) 471. [19] M. Di Vece, A. Remhof, J.J. Kelly, Electrochem. Commun. 6 (2004) 17. [20] M. Di Vece, I. Swart, J.J. Kelly, J. Appl. Phys. 94 (2003) 4659. [21] M.D. Levi, Z. Lu, D. Aurbach, Solid State Ionics 143 (2001) 309. [22] J.S. Chen, J.-P. Diard, R. Durand, C. Montella, J. Electroanal. Chem. 406 (1996) 1. [23] C. Montella, J. Electroanal. Chem. 518 (2002) 61.

46

C. Montella et al. / Journal of Electroanalytical Chemistry 608 (2007) 37–46

[24] A.V. Churikov, Russ. J. Electrochem. (Elektrokhimiya) 38 (2002) 103. [25] J.S. Chen, Thesis, Grenoble, 1992. [26] M.D. Levi, D. Aurbach, J. Phys. Chem. B101 (1997) 4641. [27] H.-C. Shin, S.-I. Pyun, in: C.G. Vayenas, B.E. Conway, R.E. White (Eds.), Modern Aspects of Electrochemistry, vol. 36, Kluwer Academy Publishers/Plenum Publishers, New York, 2003, p. 255. [28] J.-P. Diard, B. Le Gorrec, C. Montella, J. Electroanal. Chem. 471 (1999) 126. [29] J. Bisquert, G. Garcia-Belmonte, P. Bueno, E. Longo, L.O.S. Bulhes, J. Electroanal. Chem. 452 (1998) 229. [30] R. Cabanel, G. Barral, J.-P. Diard, B. Le Gorrec, C. Montella, J. Appl. Electrochem. 23 (1993) 93. [31] J. Bisquert, A. Compte, J. Electroanal. Chem. 499 (2001) 112. [32] A. Compte, R. Metzler, J. Phys. A: Math. Gen. 30 (1997) 7277. [33] P.P. Valko´, B.L. Vojta, The List, 2001; http://www.pe.tamu.edu/ valko/. [34] J.B. Floriano, E.A. Ticianelli, E.R. Gonzalez, J. Electroanal. Chem. 367 (1994) 157.

[35] K.S. Crump, J. Assoc. Comput. Mach. 23 (1976) 86. [36] Y. Ogata, T. Sakka, M. Iwasaki, J. Appl. Electrochem. 25 (1995) 41. [37] E. Barsoukov, J. Hyun Kim, J. Hun Kim, C.O. Yoon, H. Lee, Solid State Ionics 116 (1999) 249. [38] H. Stehfest, Comm. ACM 13 (1970) 47 and 624. [39] M.J. Honeychurch, Simulating electrochemical reactions with Mathematica, IBNH, St Lucia (Australia); 2006. [40] J. Abate, P.P. Valko´, Int. J. Numer. Eng. 60 (2004) 979. [41] A. Mallet, Numerical Inversion of Laplace Transform, 2000, Mathematica ‘‘Add On’’; http://www.mathsource.com/Content/Enhancements/Algebraic/0210-968. [42] P.P. Valko´, J. Abate, Comput. Math. Appl. 48 (2004) 629. [43] Y.-I. Jang, B.J. Neudecker, N.J. Dudney, Electrochem. Solid-State Lett. 4 (2001) A74. [44] P.J. Bouwman, B.A. Boukamp, H.J.M. Bouwmeester, P.H.L. Notten, J. Electrochem. Soc. 149 (2002) A699. [45] J. Abate, W. Whitt, INFORMS J. Computing 18 (2006) 408. [46] C. Montella, R. Michel, J.-P. Diard, in preparation.