Anal. Chem. 2007, 79, 6341-6347
Time-Dependent Diffusion-Migration at Cylindrical and Spherical Microelectrodes: Steady- and Quasi-Steady-State Analytical Solution Can Be Used under Transient Conditions Oleksiy V. Klymenko,† Christian Amatore,*,‡ and Irina Svir*,†,‡
Mathematical and Computer Modelling Laboratory, Kharkov National University of Radioelectronics, 14 Lenin Avenue, Kharkov, 61166, Ukraine, and Departement de Chimie, Ecole Normale Supe´ rieure, UMR CNRS 8640 “PASTEUR”, 24 rue Lhomond, 75231 Paris Cedex 05, France
Mass transport at cylindrical and spherical microelectrodes involving diffusion and migration is analyzed by means of numerical simulation under transient conditions. The origin of the intrinsic difficulties encountered during the numerical solution of the diffusion-migration equations using implicit finite differences are outlined, especially for the particular case when the number of electrons transferred equals the charge number of the electroactive species. The numerical results for transient conditions have been compared to the general analytical solutions for the current enhancement or diminishment due to migration under steady- and quasi-steady-state conditions at 1D geometry microelectrodes (Amatore, C.; Fosset, B.; Bartelt, J.; Deakin, M. R.; Wightman, R. M. J. Electroanal. Chem. 1988, 256, 255-268). This yields that the analytical limiting currents are applicable, within experimental error, to the analysis of transient diffusion-migration current responses at microelectrodes of cylindrical and spherical geometries except extremely short times after the application of the potential step, i.e., when current measurements are anyway already corrupted by ohmic drop when the supporting electrolyte concentration is low. Also, this confirms that the current enhancements or diminishments due to migration are identical for both electrode geometries when steady- or quasi-steady states are approached and do not drastically differ even under transient regimes. Miniaturization of electrodes has prompted their use for many micro- and nanoapplications. Interestingly, miniaturizing the electrode sizes gives them two valuable additional properties by providing them extremely fast time responses (i.e., very small cell time constants, RC) and making them less sensitive to ohmic drop. This has been exploited with great advantages by physical electrochemists and by electroanalysts to probe * Corresponding authors. E-mail:
[email protected]; irina.svir@ kture.kharkov.ua;
[email protected]. † Kharkov National University of Radioelectronics. ‡ Ecole Normale Supe ´ rieure. 10.1021/ac0706168 CCC: $37.00 Published on Web 07/19/2007
© 2007 American Chemical Society
and study phenomena that were inaccessible to electrochemistry before the introduction of microelectrodes.1-3 These successes among electrochemists and the simultaneous introduction of microfabrications in the analytical community have helped the design of increasingly miniaturized analytical devices in which microelectrode-based detectors may be implemented to provide precise and suitable in situ measurements. This important drive has induced a change of the demand imposed on electrochemical measurements, which progressively shifted from straight obedience to electrochemical constraints to fulfillment of electrochemical measurements under specific, and generally nonelectrochemical, conditions. Among such unusual conditions, many regard electrochemical measurements under resistive conditions imposed by the analytical separating device and protocol, since this takes full advantage of the experimental immunity of microelectrodes toward ohmic drop. However, this experimental immunity does not mean that electrochemical measurements performed under such resistive conditions have the same analytical significance (i.e., a known relationship between current and concentration) as those performed under ideal electrochemical conditions. For some applications, this may be bypassed by relying on calibrations but this is not possible for all of them. Analytical relationships have been proposed4-6 for electrodes amenable to a 1D formulation7 yet only when steady- (spherical electrode) or quasi-steady- (cylindrical electrode) state conditions are achieved. Since the achievement of quasi-steady- and steadystate conditions may require long times, these may not be suitable stricto-sensu for many analytical applications especially those (1) Michael, A. C.; Wightman, R. M. In Microelectrodes in Laboratory Techniques in Electro-analytical Chemistry; Kissinger, P. T., Heinemann, W. R., Eds.; Marcel Dekker: New York, 1996; Chapter 12, pp 367-403. (2) Amatore, C. In Electrochemistry at Ultramicroelectrodes in Physical Electrochemistry: Principles, Methods and Applications; Rubinstein, I., Ed.; Marcel Dekker: New York, 1995; Chapter 4, pp 131-208. (3) Aoki, K. Electroanalysis 1993, 5, 627-639. (4) Amatore, C.; Deakin, M. R.; Wightman, R. M. J. Electroanal. Chem. 1987, 220, 49-63. (5) Amatore, C.; Fosset, B.; Bartelt, J.; Deakin, M. R.; Wightman, R. M. J. Electroanal. Chem. 1988, 256, 255-268. (6) Oldham, K. B. J. Electroanal Chem. 1992, 337, 91-126. (7) Amatore, C.; Fosset, B. Anal. Chem. 1996, 68, 4377-4388.
Analytical Chemistry, Vol. 79, No. 16, August 15, 2007 6341
involving measurements of fast biological responses or within microfluidic systems. One alternative is to rely on the full numerical solutions of the problem at hand, taking advantage of the performances of simulation procedures and availability of fast computers. Yet, even under strict and simple electrochemical conditions, such approaches have appeared extremely difficult compared to what occurs for pure diffusional conditions and have therefore been limited to a few simple cases.8-12 Indeed, the extremely tight coupling between the migration electrical fields and concentration profiles of ionic species creates a specific numerical difficulty, in the sense that any minute error on the overall balance of ions at a given space and time point may result in drastic variations of the migration field due to Poisson’s law.13 Such difficulties become apparent through the detailed analysis we present here. In fact, it stems from the fact that the general diffusion-migration mathematical problem, which conditions any simulation numerical algorithm, consists of coupled partial derivative equations, which describe the transport of molecules and ions, plus the algebraic equation giving the electrochemical charge balance, which closes the system. Overall, one obtains a differential-algebraic (DAE) equation system whose properties differ from those of most common diffusion-reaction models usually considered in electrochemistry, which are composed of parabolic partial derivative equations (PDE) only. The present work explains in detail this point and shows why it creates specific numerical difficulties when the electrochemical reaction interconverts a charged species and a neutral one. Nevertheless, its consequence is that one cannot expect to address easily such problems by relying onto classical numerical procedures. Conversely, an analytical solution bypasses intrinsically this difficulty since it does not require any discretized approximation. However, analytical solutions may be obtained only when steadyor quasi-steady states are achieved, i.e., are not expected to be valid under conditions that depart from them. It is therefore the second purpose of this work to examine whether or not the steady- and quasi-steady-state analytical solutions may be used within a given precision to account for the treatment of experimental transient measurements involving diffusion-migration transport. THEORY In the following, we consider a cell equipped with an electrode of cylindrical or spherical geometry so that the problem is amenable to a 1D formulation for which analytical solutions4-6 have been offered. We also consider the simplest situation involving a single electroactive species, A, undergoing the redox reaction: (8) Hyk, W.; Palys, M.; Stojek, Z. J. Electroanal. Chem. 1996, 415, 13-22. (9) Palys, M. J.; Stojek, Z. J. Electroanal. Chem. 2002, 534, 65-73. (10) Palys, M. J.; Sokolowska, H.; Stojek, Z. Electrochim. Acta 2004, 49, 37653774. (11) Ciszkowska, M.; Jaworski, A.; Osteryoung, J. G. J. Electroanal. Chem. 1997, 423, 95-101. (12) Dan, C.; Van, den Bossche, B.; Bortels, L.; Nelissen, G.; Deconinck, J. J. Electroanal. Chem. 2001, 505, 12-23. (13) Smith, C. P.; White, H. S. Anal. Chem. 1993, 65, 3343-3353.
6342
Analytical Chemistry, Vol. 79, No. 16, August 15, 2007
AzA + ne- f BzB
(1)
where n ) zA - zB is the number of electrons transferred and zA, zB are the respective charges of A and its product B. Within this framework, n is necessarily equal to (1. Indeed, should |n| be larger that unity, the system would necessarily involve conproportionation reactions taking place into the solution14,15 and could not be amenable to the straight formulation given in eq 1. In addition to species A and B, one must consider also the possible counterion P of A, and the constituent of the inert 1:1 electrolyte, M+ and X-. For simplification of the present analysis, we assume that P is a monoion species, which is a common experimental situation. Additionally, we assume that all diffusion coefficients are identical since it has been shown6 that, except under special conditions where large differences between diffusion coefficients are involved, no drastic changes occur compared to the most simple case of identical diffusion coefficients. Under such conditions, it has been shown previously (see eqs 27-30 in ref 5) that one may define a dimensionless current ratio, Ψ/γ, where γ is the supporting electrolyte concentration ratio to A, which under steady- (or quasi-steady-) state conditions gives the ratio between the current in the presence of migration to that measured at the same electrode when transport by diffusion only prevails, viz., when γ f ∞: Ψ/γ ) i(γ)/i(γ f ∞). In the presence of migration, the fluxes of any species at any time and any point in the solution or at the surface of a cylindrical or spherical electrode are given by
(
ji ) -D
)
∂ci F ∂φ + zi ci ∂r RT ∂r
(2)
where i ) A, B, P, M, X; ci is the concentration of ith species, zi its charge number, r the radial coordinate in the cylindrical or spherical coordinate system, φ the electric potential, and D the diffusion coefficient assumed equal for all species. The symbols F, R, and T have their usual meaning. Double layer effects are not considered here, and it is assumed that enough electrolyte is present for a conventional double layer to be constructed. In the following we assume also that the electrode radius r0 is much larger than the Debye length (viz., r0 is larger than a few nanometers), so that, beyond the double layer,13 the solution remains electrically neutral everywhere and at any time. Therefore, the closing equation for both the cylindrical and spherical cases is
∑ zc ) 0 i i
(3)
i
at any point where the diffusion-migration equations need to be solved. The mass transport equations for the electroactive and supporting electrolyte species may thus be written as (14) Amatore, C.; Bento, M. F.; Montenegro, M. I. Anal. Chem. 1995, 67, 28002811. (15) Amatore, C.; Paulson, S. C.; White, H. S. J. Electroanal. Chem. 1997, 439, 173-182.
cylinder:
∂ci 1 ∂ )(rj ) ∂t r ∂r i
(4a)
sphere:
∂ci 1 ∂ ) - 2 (r2ji) ∂t r ∂r
(4b)
Carrying out the differentiation in eqs 4a and 4b and rearranging the terms, we note that both equations for the cylindrical and spherical electrodes are particular cases of a more general PDE taking into account the diffusion and migration effects:
[
)]
(
∂ci ∂2ci λ ∂ci F ∂2φ λ ∂φ ∂ci ∂φ ci 2 + c i , )D 2 + + zi + ∂t r ∂r RT ∂r r ∂r ∂r ∂r ∂r i ) A, B, P, M, X (5) with λ ) 1 for the cylindrical and λ ) 2 for the spherical coordinate systems, respectively. We assume that initially the electrode potential is not applied so that all species are uniformly distributed within the solution volume and the initial conditions for the DAE system 3, 5 are
t ) 0,
cA ) c0,
r gr0:
cP ) |zA|c0,
cB ) 0,
cM ) cX ) γc0,
φ ) 0 (6)
The boundary conditions valid for both the cylindrical and spherical electrodes under mass transport limited conditions at any t > 0 are as follows:
electrode (r ) r0):
cA ) 0,
jA ) - jB,
∑ z c ) 0,
Table 1. Coordinate Transformations and Nonconstant Coefficients. geometry
λ
coordinate transformation
simulation interval
f (y)
cylindrical spherical
1 2
y ) lnR y ) 1 - 1/R
[0, ln(1+6xτmax)] [0, 1]
e-2y (1 - y)4
For the numerical solution of the above PDEs, the spatial coordinate was transformed for the cylindrical and spherical cases as shown in Table 1. Thus, a uniform grid in the transformed space provides a perfectly tailored nonuniform spatial grid in the real normalized space with grid step sizes increasing away from the electrode surface. This yields the mass transport equation in the form
[
i i
jP ) jM ) jX ) 0 (7a) c A ) c 0,
cP ) |zA|c0,
φ ) 0 (7b)
where c0 is the bulk concentration of the electroactive ion A. Dimensionless Model and Coordinate Transformations. First we introduce the following normalized coordinates and variables to cast the model into the more general dimensionless form:
r , r0
τ)
Dt , r02
Ci )
ψ)
ci F ,Φ) φ c0 RT
(8)
Consider the general eq 5 describing the diffusion-migration mass transport at cylindrical and spherical electrodes. Upon substitution of the dimensionless variables and parameters, eq 8, the latter becomes
)
∂Ci ∂2Ci λ ∂Ci ∂2Φ λ ∂Φ ∂Ci ∂Φ ) 2+ + z i Ci 2 + Ci + ∂τ R ∂R R ∂R ∂R ∂R ∂R ∂R
where i ) A, B, P, M, X.
)
y)0
(11)
where S is the electrode area (S ) 2πr0L for the cylindrical electrode, where L is its length and S ) 4πr02 for the spherical electrode). It is obvious that when the electroactive species A happens to be uncharged the respective mass transport equation and the expression for the current will exactly coincide with those for diffusion-only conditions. Therefore, we will not consider the case of zA ) 0 when comparing the time behaviors of diffusionmigration and diffusion-only systems. Note that the mass transport eqs 10 in the limit of large τ (steady- or quasi-steady-state) become
d2Ci
(
(
∂CA i ∂Φ + zACA ) nFSDAc0/r0 ∂y ∂y
cB ) 0,
cX ) cM ) γc0,
R)
(10)
where the function f (y) is given in Table 1 for each situation. The simulation interval is also indicated in Table 1 being timedependent for the cylindrical case (due to the absence of real steady-state regime) and constant for the spherical one. The electroneutrality relation eq 3 and the initial and boundary conditions eqs 6 and 7 change in an obvious way upon forcing the dimensionless variables and the pertinent coordinate transformation from Table 1, so they are not given here. The normalized electrochemical current is given by the same expression for both the cylindrical and spherical electrodes:
i
infinity (r f ∞):
)]
(
∂Ci ∂2Ci ∂2Φ ∂Ci ∂Φ + z C + ) f (y) i i ∂τ ∂y ∂y ∂y2 ∂y2
dy
2
(
+ zi Ci
)
d2Φ dCi dΦ + ) 0, dy dy dy2
i ) A, B, P, M, X
(12)
which is now independent of the microelectrode geometry. The previously related analytical solutions, viz., Ψ/γ analytical expressions,4-6 correspond to the solution of these equations using the parameters recalled in Table 1.
(9) NUMERICAL SOLUTION The mass transport eqs 10 were discretized using a three-timelevel implicit finite difference scheme with second-order apAnalytical Chemistry, Vol. 79, No. 16, August 15, 2007
6343
proximation in both space and time.16 The ensuing nonlinear algebraic equations in concentrations and the electrical potential were solved using the Newton-Raphson method with the generalized Thomas method17 employed for linear systems with banded matrices. As indicated in the introduction, the numerical solution of the linear systems at hand is complicated by the fact that when migration terms are involved the matrices are not diagonally dominant and are ill-conditioned even after the appropriate reordering of unknowns. This implies that the linear solver employed here (the generalized Thomas algorithm) may give divergent solutions when the condition number of the matrix, κ(M) ) |M| |M-1|, is greater than ∼5 × 106. The latter limit was devised by running numerical experiments with different parameter values. This way it was established that when κ(M) exceeded 5 × 106 the solutions to linear systems became oscillatory and the convergence of the Newton-Raphson method was destroyed. For lower values of κ(M), the generalized Thomas method yielded meaningful solutions showing better stability properties than those predicted by the theory, as often occurs in practice.17 In fact, even though the errors in the solution of linear system may be significant when κ(M) is large, the convergence of the NewtonRaphson method could always be achieved provided that κ(M) < 5 × 106, thus ensuring the overall convergence of the numerical simulation of the mass transport problem. However, neglecting this condition resulted into severe numerical problems. The computer programs for the numerical simulation were written in Borland Delphi 7 and executed on a PC equipped with the Intel Pentium 4 processor at 3 GHz and 512 MB of RAM. To ensure a numerical convergence of better than 0.1%, in all cases reported here, a grid with 1000 nodes in the spatial coordinate and 4000 nodes in time was implemented. RESULTS AND DISCUSSION Under the conditions when the migration effects are not suppressed (i.e., when γ is comparable to or smaller than ∼100), the electric current at (quasi) steady-state is expected to follow the analytical results given in refs 4 and 5. The latter describe the ratio Ψ/γ ) iγ/i∞ of the (quasi) steady-state current due to diffusion and migration (iγ) to that due to diffusion only (i∞) as a function of zA, n, and the support ratio γ. These functional dependences for different values of zA and n, i.e., for any system described in eq 1, are summarized in eqs 27-30 in ref 5. To analyze the numerically computed results according to the diffusion-migration model given by eqs 6, 7, 10, and 11, we first normalize the current transients with respect to the diffusion-only currents as ψ/γ(τ) ) ψγ(τ)/ψdif(τ). It is expected that as τ f ∞ the current ratio ψ/γ(τ) will tend to its limiting value ψ/γ depending on the values of the three parameters zA, n, and γ. The migration effects may lead to either enhancement (nzA > 0) of diminishment (nzA < 0) of the electric current with respect to the diffusion-only conditions, which results in values of ψ/γ greater or less than unity, respectively.5 The deviation of the normalized current from unity is necessarily more pronounced (16) Richtmyer, R. D.; Morton, K. W. Difference Methods for Initial-Value Problems, 2nd ed.; Interscience Publishers: New York, 1967. (17) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, UK, 1992.
6344
Analytical Chemistry, Vol. 79, No. 16, August 15, 2007
Figure 1. (a) Time dependence of the simulated current function ζ(τ) at cylindrical microelectrode for different values of γ and zA ) -1, n ) 1. Black dots mark the values of ζ at τ ) 0.1, which are given in text. (b) Relative deviation of ζ(τ) from unity for different times (cylindrical electrode).
at low values of γ corresponding to weakly supported electrolytes. Conversely, in the limit γ f ∞, the limiting current ψ/γ is unity for any values of zA and n. A common framework for the comparison between the (quasi)steady-state and transient currents may be accomplished by normalizing again by the analytical value Ψ/γ.5 Thus, we wish to analyze the behavior of the following time-dependent function:
ζ(τ) )
ψ/γ(τ) 1 ψγ(τ) ) Ψ/γ ψdif(τ) ψ/γ
(13)
It is clear that for any values of the parameters zA, n, and γ the limiting value of ζ(τ) as τ f ∞ is 1. However, the temporal evolution of ζ(τ) is expected to depend on the values of the abovementioned parameters as well as on the microelectrode geometry. Transient Diffusion-Migration Effects at the Cylindrical Microelectrode (zA * n). Figure 1a shows the simulated time dependence of ζ(τ) exemplified for the cylindrical electrode for different values of the support ratio, γ, and zA ) - 1, n ) 1. At times close to zero, ζ(τ) is already close to unity (see solid symbols in Figure 1a) and tends rigorously toward one when the polarization time increases so that the system approaches quasi-steady-
state. The amplitude of the initial ζ(τ) value depends on the value of γ in a way such that when γ f ∞ the overall deviation span of ζ(τ) vanishes due to the system approaching the diffusion-only behavior while the deviation tends to be highest when γ f 0. This trend is further clarified in Figure 1b, which shows the distribution of the relative deviation of ζ(τ) versus unity, , defined by
) |1 - ζ(τ)|
(14)
as a function of the support value γ for three different values of the dimensionless time. Increasing the value of γ leads to an exponential convergence toward the current obtained under diffusion-only conditions (for log γ > 0). Conversely, when γ tends to zero, the deviation levels off corresponding to the unsupported electrolyte limit γ f 0. The results shown in Figure 1a for γ ) 100 and γ ) 10-4 are representative of each of the above limiting behaviors since extremely minute changes in the ζ(τ) curves are observed upon further increasing or decreasing the value of γ, respectively. Thus, the largest values of the deviation corresponding to the values of ζ marked with black dots in Figure 1a at the shortest shown time τ ) 0.1 are as follows: γ ) 10-4, ) 4.94%; γ ) 0.1, ) 4.55%; γ ) 1, ) 2.65%; γ ) 10, ) 0.49%; γ ) 100, ) 0.05%. Since the largest deviation corresponds to the limit γ f 0, which is achieved at γ ) 10-4, we will restrict ourselves to the study of chronoamperometric curves simulated with γ ) 10-4 and different values of zA while keeping n ) 1. As recalled above, the latter is conditioned by the fact that the transfer of more than one electron requires the consideration of reproportionation reaction(s) and thus necessitates a different theoretical treatment.14,15 Note that the cases for n ) -1 are readily obtained by transposing zA into -zA. The current transients ζ(τ) have been simulated for the values of zA ranging between -4 and 4 (except zA ) 1; see Numerical Difficulties in Solving the Case zA ) n (|n| ) 1)). Figure 2a shows the log-log plot of against the normalized time τ, for the cylindrical electrode. This plot allows the delineation of time scales required for the ratio between the diffusion-migration controlled current versus its pure diffusional value at the same time to converge toward its steady-state value to within a given relative error. Horizontal lines in Figure 2a indicate the limits of 1, 2, 3, and 5% deviations. It should be noted that for most experimental measurements even a deviation of 5% is usually considered acceptable owing to many other unknowns and other uncertainty sources. In this connection, we can see from Figure 2a that this level of agreement between the steady-state and transient ratios of the diffusion-migration and diffusion-only currents is achieved for any charge number value zA within the dimensionless time of ∼0.1, i.e., under conditions where the planar transport behavior is significantly approached.18 It should be noted that we have not analyzed here the double layer effects, which will inevitably alter the current transients at extremely short times while the diffusion layer is still comparable in thickness with the diffuse double layer. Indeed, the results presented in Figure 2 should be used only for the values of τ when (18) Amatore, C.; Deakin, M. R.; Wightman, R. M. J. Electroanal. Chem. 1986, 206, 23-36.
Figure 2. Deviation of the normalized current from the steady-state limit as a function of time for different charge numbers zA (shown for each curve on the right). (a) Cylindrical electrode; (b) spherical electrode.
the condition δκ-1 . 1 is satisfied, where δ ∝ (Dt)1/2 is the diffusion layer thickness and κ-1 is the Debye length.19,20 For most electrochemical conditions, κ is less than a few nanometers as soon as a classical double layer may be built so that the approximation is correct whenever the electrode dimension exceeds a few tens of nanometers19 or the scan rate is below 100,000 V‚s-1.20 Table 2 shows the values of τ1%, τ2%, and τ3% corresponding to dimensionless times when the deviation becomes less than 1, 2, and 3%, respectively, for the considered range of zA. It is evident that experimentally insignificant deviations of less than 3% from the steady-state limit are achieved even at sufficiently short dimensionless times of ∼0.3 except the case of zA ) -1 for which the value of τ3% is ∼2. Experimentally, such conditions correspond to situations where transient current components still dominate over the quasi-steady-state limit.18 Transient Diffusion-Migration Effects at the Spherical Microelectrode (zA * n). The behavior of transient currents flowing at the spherical microelectrode is qualitatively similar to that at the cylindrical electrode. However, since the time required to approach the steady-state for the sphere is much shorter than (19) Norton, J. D.; White, H. S.; Feldberg, S. W. J. Phys. Chem. 1990, 94, 67726780. (20) Amatore, C.; Lefrou, C. J. Electroanal. Chem. 1990, 296, 335-358.
Analytical Chemistry, Vol. 79, No. 16, August 15, 2007
6345
zA
τ3%
τ2%
τ1%
-4 -3 -2 -1 2 3 4
