A Theoretical Study on Turing Patterns in Electrochemical Systems

J. Phys. Chem. B 2000, 104, 6081-6090

6081

A Theoretical Study on Turing Patterns in Electrochemical Systems Nadia Mazouz and Katharina Krischer* Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, 14195 Berlin ReceiVed: January 14, 2000; In Final Form: April 1, 2000

We present theoretical studies on pattern formation in electrochemical systems with an S-shaped current potential curve (S-NDR systems) under potentiostatic control. Linear stability analysis and simulations of the reaction-migration equation give evidence that stationary patterns with a defined wavelength exist in a large parameter range. As it is the case for Turing structures, the patterns form due to an interplay of shortrange activation and long-range inhibition. It is shown that the constraint on the ratio of the diffusion constants of activator and inhibitor in reaction-diffusion equations transforms into a condition involving diffusion and migration lengths of the system. This condition is fulfilled in practically all electrochemical systems. The experimental parameters under which the patterns form should be readily accessible.

Introduction In his pioneering 1952 paper,1 Turing predicted the existence of regular stationary concentration patterns in some chemical systems maintained far from the thermodynamic equilibrium, a phenomenon that exerts a continuing fascination. The novel and most astonishing aspect of the mechanism was that diffusion, that is, a transport process that was known to smooth any concentration gradient, can destabilize a homogeneous steady state that is stable in the absence of diffusion. To form this nonequilibrium structure, certain requirements on the reaction kinetics and the transport properties have to be met. The kinetics has to be of the well-known activator-inhibitor type,2,3 and indeed such kinetics has been found in a large variety of chemical systems [see e.g., refs 4 and 5]. Given such an activator-inhibitor system, a Turing structure can develop if the inhibitor diffuses faster than the activator. In chemical systems, however, this latter requirement turned out to be an obstacle for the experimental verification because the diffusion constants of activator and inhibitor, mostly ions, are typically very similar. Thus, despite the comparatively large number of chemical reactions known to obey activator-inhibitor kinetics, the experimental search for Turing patterns was successful only 40 years after the derivation of the theoretical concept.6,7 In most electrochemical systems that exhibit typical nonlinear behavior, the interfacial potential of the working electrode acts as autocatalytic variable (see e.g., the recent reviews8,9). Here, the source of the instability lies in a negative impedance characteristic of a faradaic process. In general, two cases are distinguished. In the first case, the stationary polarization curve exhibits a region of negative differential resistance (NDR). In the second one, the latter is hidden in the steady-state currentpotential curve by a further potential dependent and inherently slow process (HNDR, hidden negative differential resistance). But, independent of whether the NDR is visible in the stationary current-potential curve or not, it causes an instability that can be traced back to an N-shaped current-potential characteristic. The coupling of such a characteristic to a slow, often nonelectric process brings about oscillations, or more generally, temporal * Corresponding author. E-mail: [email protected]. URL: http://w3.rz-berlin.mpg.de/∼nadia/.

dynamics typical for an activator-inhibitor system where the potential takes on the role of the activator. The inhibitor is then a chemical quantity, that is, the concentration of the electroactive species or the coverage of the electrode by an adsorbate. Spatial inhomogeneities in activator and inhibitor induce migration and diffusion fluxes, respectively, which in turn lead to the occurrence of electrochemical pattern formation. Since the typical time scale associated with migration is much shorter than the one associated with diffusion, we are dealing with a fast activator and a slow inhibitor transport process. The spatio-temporal dynamics is dominated by the fast migration, that is, by fast spatial coupling of the activator, and the potential waves observed in the systems can be understood in terms of reactionmigration systems. The main difference from a reactiondiffusion system with fast activator diffusion is that, in contrast to diffusion, migration represents a long-range spatial coupling.10,11 The latter fact has been shown to be responsible for the acceleration of wave phenomena characteristically found in experiments.12-15 Thus, according to the current knowledge of pattern formation in electrochemical systems, the existence of stationary patterns in electrochemistry that emerge in a Turing-like mechanism seems to be very unlikely. The scope of this paper is to demonstrate that this conjecture is not true. On the contrary, we demonstrate that there exists a class of electrochemical systems in which a Turing-like instability destabilizes the homogeneous steady state in a large parameter region. Furthermore, we show that in contrast to chemical systems, in electrochemical systems the conditions on the parameters governing the spatial coupling are practically always met. Hence, whenever the requirement on the reaction kinetics is fulfilled, in this class of electrochemical systems stationary patterns are just as likely to occur as temporal oscillations in systems with an N-shaped current-potential curve. The condition on the reaction kinetics is that the current-potential characteristic possesses a bistable region whereby the autocatalytic branch must possess a negative differential resistance.16 The easiest, though not the only, realization of this mechanism gives rise to an S-shaped current potential characteristic. Zn deposition represents an early experimental example in which an S-shaped current-potential characteristic was found.17,18 Obviously, in

10.1021/jp000203+ CCC: $19.00 © 2000 American Chemical Society Published on Web 06/07/2000

6082 J. Phys. Chem. B, Vol. 104, No. 25, 2000

Mazouz and Krischer

these systems, the current is not a unique function of the potential, and the autocatalysis responsible for the multistability must stem from a variable different from the interfacial potential. This variable takes on the role of the activator, and the interfacial potential acts like an inhibitor. To distinguish the two classes of systems with a negative differential resistance (NDR) in the current-potential curve, we use in the following the terms N-NDR (HN-NDR) and S-NDR according to the form of the current-potential characteristic. After reviewing the general reaction-migration equation governing the potential evolution in a spatially distributed electrochemical system, we discuss in the following first the homogeneous dynamics of an electrochemical system with any S-shaped current-potential characteristic and then present a prototype electrochemical model possessing such a behavior. Then, the linear stability of the homogeneous state with respect to spatial perturbations is analyzed, again, first for the general case, that is, without specifying a reaction mechanism. Subsequently, bifurcation diagrams and simulations of the model equations are presented. Finally, the specific “electrochemical” features of the models and of the stationary patterns are discussed.

direction). (Note that eq 1 does not contain a spatial operator that depends on x.) Finally, we have to relate the double layer potential, φDL, and the potential in the electrolyte, φ. Setting the potential at the reference electrode (z ) 0) to zero,19 under potentiostatic operation conditions (which are exclusively considered in this paper) this relation is given by

U ) φDL + φ|z)-w

where U is the externally applied potential. As pointed out in the Introduction, we considered systems that possess an S-shaped current-potential curve. In this case, the state of the system can never be expressed as a function of the double-layer potential only. Rather, the multivalued region of the current-potential curve is caused by some autocatalytic kinetics of a second variable u (or the interplay of several variables20). At infinite conductivity (or zero resistance), the temporal evolution of the electrochemical system is then described by the dynamics of this second variable only, which, in general terms, can be expressed as

∂u ) f(u,φDL) + D∇2u ∂t

Model and Linear Stability Analysis The central variable in electrochemical pattern formation is the interfacial potential, φDL. The basic equation governing the spatio-temporal evolution of φDL results from the charge balance at a given point of the electrode:

|

∂φDL ∂φ ) - ireac(u,φDL) - σ ∂t ∂z z)-w

C0

(1)

Here, C0 is the specific double-layer capacity (which is taken to be time-independent), and t is time. ireac denotes any faradaic current density and u a variable on which the current depends, such as the concentration of an electroactive species or the coverage of the electrode with an adsorbate. σ is the specific conductivity of the electrolyte divided by the electrode area. z is the spatial direction perpendicular to the electrode, where z ) -w is at the electrode, or more precisely, at the end of the electric double layer from where on the electrolyte is electroneutral. φ stands for the electric potential in the electrolyte outside the electric double layer. Thus, the last term on the righthand side (rhs) of eq 1 describes the migration current density “entering” the double layer at a given point of the electrode. It is balanced by the faradaic current density and the capacitive current density, the first term on the rhs and the term on the lhs, respectively. For a discussion of this equation in the context of pattern formation see refs 9-11. To solve eq 1, it is necessary to determine the electric field at the double layer perpendicular to the electrode, which is possible once the electric potential φ inside the electrolyte is known. Due to the electroneutrality of the electrolyte, φ is very well approximated by Laplace’s equation

∆φ ) 0

(1a)

where ∆ is the Laplace operator. For any spatially inhomogeneous potential distribution at the double layer, the solution of eq 1a has electric field components parallel to the electrode, which thus induce migration currents parallel to the electrode. These enable a spatial “communication” or coupling between different locations of the electrode and are thus responsible for the observed pattern formation along the electrode (our x

(1b)

(2)

Here, f is a nonlinear function that is autocatalytic in u. Furthermore, it has been assumed that the only relevant transport process parallel to the electrode is diffusion (either surface diffusion or diffusion in the Helmholtz layer), and D is the respective diffusion constant. All possible transport processes perpendicular to the electrode are lumped in the function f. In the following, we discuss nonequilibrium patterns in S-NDR systems that are described by Equations 1-2. Homogeneous System. First consider the homogeneous case, in which eqs 1 and 2 reduce to

du ) f(u,φDL) dt C0

dφDL σ ) -ireac(u,φDL) + (U - φDL) dt β

(3) (4)

where β is a geometric factor such that the cell resistance for an electrode with unit area is given by β/σ. The second term on the rhs of eq 4 reflects that, for a uniform distribution of the double-layer potential, the migration current density at every point of the interface (the second term on the rhs of eq 1) can be expressed as the potential drop between working and reference electrode divided by the cell resistance. Systems described by eqs 3 and 4 may undergo a Hopf- or a saddlenode bifurcation as a parameter is varied. They thus exhibit bistability and oscillations for some region in parameter space. The parameter region in the U-σ plane, in which the system possesses multiple steady states, can be easily rationalized when plotting the polarization curve together with the so-called loadline, the second term on the rhs of eq 4, in one plot (Figure 1a). The intersections of both curves are stationary states. As long as the cell resistance is smaller than the absolute value of the (zero frequency) faradaic impedance, [(direac/dφDL)ss]-1 at the intersection ss, that is, as long as

| ||

direac σ > β dφDL

(5)

ss

the steady state is unstable and two further intersections exist, corresponding to stable stationary states for most parameter

Turing Patterns in Electrochemical Systems

J. Phys. Chem. B, Vol. 104, No. 25, 2000 6083 values. When the resistance becomes larger than the faradaic impedance at the inflection point of the polarization curve, the curves intersect in only one point and the system possesses a single steady state. Consequently, measurements of the current as a function of the applied potential U exhibit a hysteresis for small values of β/σ, whereas the S unfolds with increasing electrolyte resistance, leading to a continuous I/U diagram with a positive slope over the whole range (Figure 1b). The bistable region in the U-(β/σ) plane is indicated by the solid line in Figure 1c. It is broadest at vanishing electrolyte resistance, becoming smaller and being shifted to larger values of the external voltage with increasing resistance. Note that condition (5) for the existence of multiple steady states is exactly opposite to the corresponding one for a system with an N-shaped current potential curve, and all the other mentioned properties are complementary as well:21 N-NDR systems are bistable for infinite resistance and become monostable with decreasing resistance, with a correspondingly inflected bistable region in the U-(β/σ) plane. If the autocatalytic variable u changes on a faster time scale than the double layer potential, eqs 3 and 4 possess limit cycle solutions. The oscillations exist in most of the parameter region in which the current-potential curve and load line have only one intersection on the unstable (autocatalytic) branch of the polarization curve. In the U-(β/σ) parameter plane, the Hopf bifurcation occurs along a V-shaped region whose tip is close to the tip of the bistable region and which opens toward larger values of resistance and voltage (Figure 1c, dashed line). To be able to carry out numerical simulations, we have to introduce a specific model for an electrochemical system with an S-shaped current potential curve. The choice of the equations below was motivated by experiments that are currently being carried out in our laboratory. The model describes an electrontransfer reaction that is inhibited by an adsorbate that undergoes a first-order phase transition when the potential is varied. Examples for such adsorbates are thymine, uracil, or camphor on Au single-crystal electrodes.22-25 A detailed analysis of the dynamic behavior of the model can be found in ref 26. We take f to be

f(u,φDL) ) kadcA(1 - u)ew˜ (φDL)e-gu/(2RT) - kdue-w˜ (φDL)egu/(2RT) (6) with

1 1 (7) w˜ (φDL) ) - (C0 - C1) φ2DL + C1φDLφsh 2 NmaxkBT

(

Figure 1. (a) Schematic of an S-shaped current-potential characteristic (dashed line) and a load line (solid line). (b) Current-voltage curve resulting if the cell resistance is such that the current-potential characteristic and load line have a single intersection in the whole potential range (as it is the case in (a)). (c) General locations of saddlenode bifurcations (solid line) and Hopf bifurcations (dashed line) in the resistance R (R ) σ/β)/voltage U parameter plane for an S-NDR system.

)

f then models adsorption and desorption of a neutral (organic) molecule with attractive interactions in the Bragg-Williams approximation following Frumkin’s early suggestions.27 u can thus be interpreted as a number of adsorbed molecules (or atoms), where Nmax denotes the maximum number of molecules that can adsorb on a unit area and kB the Boltzmann constant. The two terms in w˜ (φDL) account for the work expended in the adsorption process due to a decrease in the dielectric constant and to an increase in the double layer thickness upon adsorption on one hand, and the energy associated with the dipole moment of the molecules adsorbed, on the other hand. (φDL is measured versus the point of zero charge.) The specific capacities of the bare electrode and the adsorbate-covered electrode are denoted by C0 and C1, respectively. g represents an interaction energy of the molecules. For g < -4, f describes a multivalued, S- or Z-shaped isotherm where the two outer branches of the isotherm correspond to a condensed and a dilute phase of the adsorbate.

6084 J. Phys. Chem. B, Vol. 104, No. 25, 2000

Mazouz and Krischer specific reactions that give rise to the S-shaped polarization curvesthe linear stability of a homogeneous state in the spatially extended system described by eqs 1 and 2. Therefore, we first have to specify boundary conditions (b.c.) to eq 1a that depend on the cell geometry. Consider a simple geometry that already served as a caricature of an electrochemical experiment in the context of pattern formation in N-NDR systems:10,11,28,29 The electrode is modeled as a one-dimensional domain with periodic b.c., an idealized description of a ring electrode whose diameter is large compared to the thickness of the ring. The electrolyte is modeled on a cylindrical surface with a constant potential at some distance from the electrode, at which the reference electrode is located. This choice of b.c. allows us to solve Laplace’s equation analytically. Before we do so, it is advantageous to introduce the dimensionless variables

Figure 2. (a) Coverage of a neutral species, u, as a function of the electrode potential as described by eq 6. The multivalued region results from the attractive interactions between the adsorbed molecules. The isotherm was calculated with the following dimensionless parameters (see eq 32): γ ) -4, p ) 0.5, ν ) 2, C ) 0.1. (b) Current-potential characteristic in the absence (dashed line) and in the presence of the adsorbate (solid line) according to eq 8 (eq 33). In the latter case, the characteristic possesses an S- and a Z-shaped region.

It, thus, implies the existence of a thermodynamic phase transition. A typical coverage/potential dependence is displayed in Figure 2a. Suppose that the adsorbate itself does not react electrochemically but that it inhibits the electron transfer of a faradaic reaction that proceeds via Butler-Volmer kinetics on the free surface. Then (assuming a large enough overvoltage such that the backward reaction can be neglected), the reaction current density is modeled by

ireac ) (1 - u)nFcrkreaceδ(φDL-V) with δ )

RnF RT

(8)

where n is the number of electrons per molecule transferred during the reaction and R the transfer coefficient. cr stands for the bulk concentration of the electroactive species, kreac for the rate constant, and V for the equilibrium potential of the reaction. In the two regions in which the isotherm is multivalued, the stationary current-potential curve possesses the form of an S or Z, respectively (Figure 2b). For the specific situation depicted above, the double layer capacitance depends on the coverage u. This gives rise to an additional term in the evolution equation of the double layer potential. As it is the aim of this paper to discuss the phenomena that can be expected in a general S-NDR system and this complication does not seem to be generic, we neglected the coverage dependence of C0 in our model and approximated the evolution equation for φDL by

C0

dφDL σ ) -ireac(u,φDL) + (U - φDL) dt β

(9)

However, we confirmed that the results presented here are qualitatively not affected by this simplification. Of course, quantitatively the location of the bifurcations change to some extent, especially the parameter range in which oscillations occur becomes considerably smaller. Spatially Extended System. Our next objective is to investigatesfor a general case, that is, independent of the

xf

2πx RnF u z 4π2D t f t, fu , zf , φ f φDL, L w RT DL u0 L2 (10)

where the dimensionless quantities are given on the rhs of the expressions. x is the direction parallel to the electrode and L the length of the electrode in physical space. In general, u0 will be a bulk concentration or the maximum number of adsorbed molecules per unit area, depending on whether u represents the concentration of a solution species or the number of adsorbed molecules per unit surface. The transformation of all potentials (U, V, φ) in this paper is analogous to the one of the doublelayer potential φDL given above. With the definitions given in eq 10, Laplace’s equation reads

∂ 2φ ∂2φ ) - β2 2 , x  [0,2π), z  [0,-1] 2 ∂z ∂x

(11)

where x is the direction parallel to the electrode, and β ) 2πw/ L. z ) -1 denotes a location at the electrode and z ) 0, correspondingly, one at the equipotential plane. In our geometry, it is subject to the following boundary conditions:

φ(x + 2π,z,t) ) φ(x,z,t) φ(x,z ) 0,t) ) 0 φ(x,z ) -1,t) ) φ0

(12)

The solution of Laplace’s equation is a Fourier series parallel to the electrode that decays into the electrolyte for increasing values of z according to ∞

φ(x,z,t) )

∑[(c˜ n(t)cos(nx) + d˜ n(t)sin(nx))sinh(nβz)] + n)1 c˜ 0(t)z (13)

where n is the wavenumber. Finally, we translate eqs 1 and 2 in dimensionless form:

∂u ∂2u ) µf(u,φDL) + 2 ∂t ∂x

( |

(14)

∂φDL d d ∂φ ) κg(u,φDL) + (U - φDL) + φ|z)-1 ∂t β β ∂z z)-1 (15)

)

Here, β has the same meaning as in eq 11, and the parameters µ, κ and d relate to the original problem in the following way:

Turing Patterns in Electrochemical Systems

µ)

Rn2F2L2kre-V L2 Lσ , κ ) (16) k , d) u 2 RTC04πD 2πDC0 4π D

ku and kr represent characteristic rates involved in the reaction kinetics of u and the electron-transfer kinetics, respectively and have the unit [s-1] and [mol/(s cm2)]. Equation 15 has been written such that the last term on the right-hand side vanishes for a homogeneous potential distribution (cf. Equation 3 and ref 11). This term describes the spatial coupling through cross migration currents in the electrolyte, and it is thus the analogous term to the diffusion term in eq 14. Knowing the solution of Laplace’s equation for our geometry (eq 13), it is easy to see that the functions

φˆ ) (c˜ n(t)cos(nx) + d˜ n(t)sin(nx))sinh(nβz)

(17)

solve eq 11 and are eigenfunctions of the equation

∂ φˆ ) mφˆ ∂z

(18)

with m ) nβ coth(nβz). For our periodic boundary conditions, the eigenfunctions uˆ of the spatial eigenvalue problem

∂2uˆ ) -n2uˆ ∂x2

(19)

uˆ ) ancos(nx) + bnsin(nx), n ) 0,1,2,...

(20)

are Fourier functions

When realizing, finally, that after expanding the double-layer potential in a Fourier series ∞

φDL(x,t) )

∑[cn(t) cos(nx) + dn(t)sin(nx)] + c0(t) n)1

(21)

the potentiostatic control condition

(22)

connects the Fourier coefficients in the expansion of φDL (eq 21) and of φ|z)-1 in the following way:

c0 ) U - c˜ 0 cn ) -sinh(nβ)c˜n, dn ) -sinh(nβ)d˜ n

(23)

we can derive the linearized problem. Let (uss,φss DL) be a homogeneous steady-state such that

)|

( ) )( ) ( ) ( )

(

(uss,φss DL)

u φDL -

u u n2 0 ˜ u ) J˜ss φ -D φDL 0 dk(n) φDL DL

Perturbing this homogeneous state by small local perturbations, we set

)

(25)

Then, for small values of u˜ and φ˜ DL, eqs 14 and 15 become

(26)

with

k(n) ) ncoth(nβ) - β-1 ˜Jss is the Jacobian matrix of the homogeneous system evaluated at steady state. The stability of a homogeneous steady state is determined by the real part of the eigenvalues λn of the characteristic polynomial

˜ (n)| ) 0 |J˜ss - λn˜I - D

(27)

Here, ˜I stands for the identity matrix. If the real part of all eigenvalues Re(λn) is negative for any n, then the perturbations decay exponentially in time, and the homogeneous stationary solution is asymptotically stable. If, on the other hand, one eigenvalue has a positive real part for any wavenumber n, perturbations with that wavenumber grow exponentially in time, and the homogeneous system is unstable. Then, the system evolves to a new patterned stationary state. The homogeneous state will loose stability if Re(λn) ) 0. To change sign for the real part of λn for any n * 0, the determinant ∆ and the first derivative d∆/dn of the linear operator L must simultaneously be 0, where L is given by

L)

(

µfu - n2 µfφDL κgg

κgφDL -

d - dk(n) β

)|

(28) (uss,φss DL)

Let us consider under which conditions ∆ ) 0 can be fulfilled. First, recall that a prerequisite for the occurrence of an instability is that the homogeneous steady state lies on the autocatalytic branch of f, that is, fu > 0. Furthermore, we require that the homogeneous steady state is asymptotically stable and thus that

(29)

where Tr and Det denote trace and determinant of J˜ss, respectively. Consequently, a Turing-like instability will only occur if the stationary state of the system lies on the NDR branch of the S-shaped current potential curve and if that stationary state is not a saddle point, as for most parameter values within the bistable region. With this in mind, one can show that a necessary condition for a spatial symmetry breaking is that

d>

d ss ss ss f(uss,φss DL) ) 0 and κg(u ,φDL) + (U - φDL) ) 0 (24) β

( )(

(

µfu µfφDL ) κg d u κgφDL β

Tr ˜Jss < 0 and DetJ˜˜ss > 0

U ) φDL ) φ|z)-1

u˜ u - uss ) φ˜ DL φDL - φss DL

( ) ∂u ∂t ∂φDL ∂t

J. Phys. Chem. B, Vol. 104, No. 25, 2000 6085

n2 n2 ) k(n) ncoth(nβ) - (1/β)

(30)

This inequality represents the electrochemical analogue to the requirement that for a RD system to undergo a Turing bifurcation, the ratio of the diffusion constants of inhibitor and activator has to be larger than one.2 In contrast to the latter condition, inequality (30) depends on two quantities, the wavenumber n and the aspect ratio of the cell, β, which is a measure of the nonlocality of the coupling.11 Let us investigate how these parameters affect the necessary condition on d. The rhs of the inequality (30) is always larger than 1 for n g 1 irrespective of the value of β. Moreover, it is a monotonically increasing function of n. This means that the lower limit for d

6086 J. Phys. Chem. B, Vol. 104, No. 25, 2000

Mazouz and Krischer

is always larger than one, and it is larger the larger the critical wavenumber is. The implications of the value of β on the minimal d can be seen best if we consider two limiting cases. For large aspect ratios β, k(n) f n, whereas for small values of β, we are close to the diffusive limit in which k(n) f (β/3) n2.11 Thus. inequality (30) becomes

d > n for β f ∞ and d >

3 for β f 0 β

(31)

Hence, in the diffusive limit the condition becomes independent of n, as one would expect. However, because the condition only holds if β , 1, for decreasing coupling range the necessary value of d becomes larger. This is a consequence of the fact that in this region the “electrochemical diffusion constant”, which can be identified with Del ) σβ/3, decreases with decreasing β. As a rule of thumb, we can say that the parameter region in which the homogeneous state is unstable with respect to spatial perturbations becomes smaller for increasingly localized coupling. Remember, however, that changing β also changes the homogeneous steady state, which makes things more complicated than it appears from inequality (30) alone. Inequality (30) is necessary but not sufficient for a positive real part of some λn. In principle, one can derive further conditions that restrict possible parameter combinations. However, these conditions are rather complex; they involve, for instance, in the long-range limit, a cubic equation in n. We refrain from giving the formulas here. Instead, we demonstrate typical properties of eqs 14 and 15 with numerical results. Therefore we set

f ) (1 - u)e(w(φDL)-γu) - pue(-w(φDL) + γu) g ) -(1 - u)e

φDL

(32) (33)

with

w(φDL) ) -ν(1 - C)φ2DL, ν ) γ)

g 2RT

C0R2T

, C)

2NmaxR n F kB 2 2 2

p)

C1 C0

kd kadcA

f and g are the dimensionless forms of eqs 6 and 8, where the second term in w(φDL) ( eq 7), that is, the change of the energy of the double-layer condenser owing to the dipole moments of the molecules adsorbed, has been neglected. The numerical integration of eqs 1 and 2 were done employing the pseudospectral method as described in refs 11 and 28. The integration of the resulting set of ordinary differential equations for the temporal evolution of the coefficients of the Fourier functions in eqs 20 and 21 was performed with the Livermore solver for ordinary differential equations (LSODE).30 The bifurcation diagrams were calculated with AUTO.31 Results Consider a state of the system that is close to the cusp but outside the bistable region (cf. Figure 1c). Typical dispersion relations, that is, curves of the real part of the eigenvalue λ(n) of eq 27 as a function of the wavenumber n, are displayed in Figure 3 for three different values of d. In the bottom curve, all eigenvalues are negative, and thus the homogeneous state is asymptotically stable with respect to any small perturbation. When increasing d, the maximum of the curves is shifted toward

Figure 3. Dispersion relation displaying the growth rate of the perturbations, λ(n), versus their wavenumber n for three homogeneous stationary states of eqs 14 and 15 close to the cusp (cf. Figure 1c or 4) with f and g as defined in eqs 32 and 33. From bottom to top: d ) 150, 185, 213; U ) 1.5; remaining parameter values as in Figure 4.

more positive values, crossing zero at a critical value dcr that is close to the middle curve. For d somewhat larger than dcr (top curve), some eigenvalues possess a positive real part, and any perturbation of the homogeneous state with the corresponding wavenumbers will grow, the system eventually settling down in a patterned state. Thus, at dcr the homogeneous state loses stability in a bifurcation that, although the underlying transport process is not described mathematically by a diffusion term, bearssas we will discuss belowsthe characteristics of a Turing bifurcation. We therefore call it an electrochemical Turing bifurcation. A Turing structure arises from a balance between local activation processes due to the autocatalytic reaction of the activator and long-range inhibition provided by fast diffusion, or more generally a fast transport process of the inhibitor. In our system (eqs 14 and 15), d is a measure of the ratio of typical time scales involved in the communication across the electrode by migration and diffusion, and one could expect that a homogeneous state is unstable whenever d > dcr, as is true for RD systems. However, as already discussed above, a peculiarity of electrochemical systems is that the homogeneous steady state also depends on the conductivity that defines the time scale of migration. Consequently, in electrochemical systems, there is not only a lower bound for d but also an upper bound above whichsdue to the properties of the homogeneous dynamicss the homogeneous state becomes asymptotically stable again. The set of electrochemical Turing bifurcations in the U/d-1 parameter plane is mapped out in Figure 4, together with the location of the saddle node bifurcations of the homogeneous system. Note that d plays in the homogeneous system the role of the specific conductivity of the electrolyte (cf. Figure 1c). The set of Turing bifurcations loops around the cusp and marks the region where no homogeneous stable solutions exist at all. Its location is similar to the location of the Hopf bifurcation set (which exists for larger values of µ or smaller values of κ). Contrary to the Hopf bifurcation set, the Turing-unstable region is bounded at large values of U and d-1. This stems from the weaker longrange inhibition, that is, slower migration, for small values of d.

Turing Patterns in Electrochemical Systems

J. Phys. Chem. B, Vol. 104, No. 25, 2000 6087

Figure 4. Two-parameter bifurcation diagram displaying the locations of the saddle node (solid line) and Turing bifurcations (dot-dashed line). The curve connecting the square symbols gives the approximate location of a saddle-node bifurcation of nontrivial stationary solutions that coexist with the homogeneous state. Their existence region extends into the Turing-unstable region. This second boundary is not shown. Parameter values: γ ) -4, ν ) 2, p ) 0.5, C ) 0.1. β ) 6.28, µ ) 10, κ ) 10.

Figure 6. (a) Asymmetric stationary stable potential profile coexisting with the homogeneous solution at d ) 100, U ) 1.75; remaining parameters as in Figure 4. (b) Asymmetric stationary stable potential profile coexisting with a symmetric one that emerges in a supercritical Turing instability. Parameter values: d ) 100, U ) 1.8; remaining parameters as in Figure 4. Figure 5. Stationary potential patterns for d ) 50 and U ) (from top to bottom) 2.7 (homogeneous solution on the cathodic side), 2.8, 2.9, 3.1, 3.3, 3.4, 3.5, 3.6 (homogeneous solution on the anodic side); remaining parameter values as in Figure 4.

Spatial profiles of the stationary solutions at d-1 ) 0.02 for different values of U are reproduced in Figure 5. The two straight lines indicate the homogeneous potential values just outside the Turing-unstable region. When crossing the Turing bifurcation, sinusoidal period-2 patterns emerge with small amplitude. A further change in U leads first to an increase in the amplitude of the pattern, which then decreases again close to the second Turing bifurcation, whereas the homogeneous component of the solution monotonically shifts toward the value of the stable solution beyond the second Turing bifurcation. At lower values of d-1, a stable inhomogeneous asymmetric potential distribution exists already outside the Turing-unstable region and coexists with the stable homogeneous solution (Figure 6a). It must be born in a saddle-node bifurcation whose approximate location is also depicted in Figure 4. The existence range of the asymmetric solution extends into the Turingunstable region, where it then coexists with the more sinusoidal pattern that is born in the Turing bifurcation (Figure 6b). Above, we discussed that upon increasing the time scale of the activator, which corresponds in our case to a larger value of µ or, equivalently, to put it in chemical terms, increasing the adsorption rate of the reaction inhibiting species, the homogeneous state becomes unstable with respect to homogeneous perturbations in a Hopf bifurcation. Then, the emerging limit

Figure 7. Gray scale representation of the evolution of the coverage in a Turing- and Hopf-unstable region for two different initial conditions: left, the homogeneous steady state was slightly perturbed in the first three modes; right, random initial conditions. The nonlinear gray scale is displayed at the bottom of the figure. Parameter values: µ ) 100, d ) 100, U ) 2; remaining parameters as in Figure 4.

cycle coexists with the stationary pattern in a comparatively large region in parameter space which, in a first approximation, is the intersection of the Turing- and Hopf-unstable areas. Such coexisting solutions are displayed in Figure 7, and the locations of the Hopf and Turing bifurcations in the U/d-1 plane are shown in Figure 8a. As is evident in the blowup, at large values of d (small d-1), the homogeneous solution becomes unstable first in a spatial mode and then in the homogeneous mode. For lower values of d, the bifurcation lines cross and the homogeneous oscillations exist prior to the stationary patterns. The critical wavenumber changes along the Turing bifurcation set, as shown in Figure 8b. Here, the upper line of the loop

6088 J. Phys. Chem. B, Vol. 104, No. 25, 2000

Figure 8. (a) Two-parameter bifurcation diagram displaying the locations of the saddle node (solid line), Turing (dot-dashed line), and Hopf (dashed line) bifurcations. The inset shows a magnification of the diagram close to the origin. µ ) 100, remaining parameters as in Figure 4. (b) Critical wavenumber along the Turing bifurcation as a function of U. The upper part of the curve corresponds also to the upper part of the Turing set shown in (a).

Figure 9. Gray scale representation of the long-term behavior of the coverage (top) and the potential (bottom) in the parameter range of mixed-mode oscillations. The calculation was started on the homogeneous stationary state with a small perturbation in the forth mode. Parameter values: µ ) 1000, d ) 100, U ) 2; remaining parameters as in Figure 4.

corresponds also to the upper line in Figure 8a and vice versa, such that the larger wavenumber is observed at smaller values of d. An increase of µ with all other parameters, as in Figure 7, gives rise to the spatiotemporal dynamics shown in Figure 9. The Turing pattern is not stationary anymore but oscillates in time. This mixed-mode oscillation coexists with the homoge-

Mazouz and Krischer

Figure 10. (a) Critical value of d as a function of the system length. (b) Dependence of the wavelength of the stationary structure on the system length.

neous limit cycle. Different types of mixed-mode oscillations close to a Turing bifurcation have been reported in the literature and are compiled in ref 32. The perhaps most striking property of Turing structures is that the wavelength of the pattern does not depend on the system length but is completely determined by the rate constants of the homogeneous system. In our electrochemical model, it was not possible to define one scale parameter, as in RD-systems. Rather, L enters, to different powers, in all four parameters, µ, κ, d, and β. To test for the influence of the system length on the wavelength of the spatial structures, we took the length dependence of three of the parameters explicitly into account and continued the Turing bifurcation with L as bifurcation parameter and dcr as adjustable quantity. The outcome of such a calculation is shown in Figure 10, where the dependence of dcr as well as of the critical wavelength are displayed versus the system length L. dcr linearly increases with L and reflects the functional relation of both parameters (eq 16).33 In accordance with this linear relation between dcr and L, we observe that for a given steady state a change of the length of the system, L, does not affect the wavelength of the nonequilibrium structure. It thus justifies the notation electrochemical Turing structure for the stationary patterns we discuss in this paper. Discussion In the preceding sections, we considered a class of electrochemical systems, so-called S-NDR systems, which are

Turing Patterns in Electrochemical Systems characterized by a multivalued current potential characteristic whereby the slope of the unstable (autocatalytic) branch of the characteristic has to be negative (see footnote 3). We presented a linear stability analysis and simulations of an electrochemical model that can be regarded as a generic model for an S-NDR system. Both approaches gave evidence that in these systems, stationary nonequilibrium structures exist due to an interplay of short-term activation and long-range inhibition. Their wavelength is defined by the parameters that determine the local dynamics. As these are just the properties of Turing structures, we can conclude that, indeed, Turing patterns can also be expected in electrochemical systems. The next question is, given such an S-NDR system, how likely is it that Turing patterns form, or, in other words, how large is the parameter range in which the stationary patterns exist, and how easily can this range be accessed experimentally? First, there are two conditions on the local dynamics, namely, (a) that the intersection between load line and current-potential characteristic has to be on the NDR branch of the latter, and (b) that the absolute value of the slope of the polarization curve at the fixed point is larger than the slope of the load line in the I/φDL plane. It is easy to see that these requirements do not constitute a restriction. A change of the externally applied voltage U leads to a parallel shift of the load line. Thus, the first inequality will always be fulfilled in some voltage interval. The condition on the cell resistance is similar of the one necessary to observe bistable behavior in N-NDR systems. A low electrolyte conductivity and a large current, that is, a high concentration of the electroactive species and a larger electrode area, are favorable here. It seems to be always possible to adjust the concentrations of base electrolyte or the electroactive species accordingly. Whether or not the local dynamics supports pattern formation is apparent from a cyclic voltammogram. The current-voltage dependence should be unique with a sigmoidal form, as depicted in Figure 1c. Whether the sigmoidal form in fact arises from an S-shaped current-potential curve can easily be tested. An increase in conductivity should recover the multivalued region and hence the bistability in the cyclic voltammogram. Now, we are left with the last condition (inequality (30)), which is a mathematical expression for the mechanistic requirement of a fast long-range inhibition. With the definition of d (eq 16), we can estimate typical experimental values of d. Taking as a lower bound of the conductivity σ ) 1.7 × 10-5 Ω cm-1, which corresponds to a 10-4 M solution of a 1:1 salt, and setting L ) 1 mm, D ) 10-5 cm2/s, and C0 ) 20 µF/cm2, d is of the order of 103. Considering the two limiting cases of diffusive (β f 0) and long range (β f ∞) coupling, inequality (30) becomes d > 3/β and d > n, respectively. (Remember that β is the aspect ratio of electrode circumference and 2π times the distance between working and reference electrode.) Thus, for the diffusive coupling, inequality (30) is fulfilled as long as β is larger than 3 × 10-3, which is true for any experimental situation. In the long-range limit, d > n, which holds for the parameters given above as long as the characteristic wavelength of the pattern L/n is larger than 1 µm. It is not possible to derive a useful analytical expression for the critical wavelength. However, our simulations suggest that, in a typical experimental situation, one should be well above this threshold: Assuming room temperature and the same values for C0, D, and L as in the estimation above, µ and κ are given by µ ) 2.5 × 10-5 kad and κ ) 5 × 106 kr exp [-V], respectively. Thus, in the simulations shown in Figures 4 and 8, kad ) 106 or 107 [s-1], kr exp [V] ) 0.2 × 10-5 [mol s-1 cm-2], and the wavelength of

J. Phys. Chem. B, Vol. 104, No. 25, 2000 6089 the corresponding patterns amounts to 500 µm, more than 2 orders of magnitude above the lower limit for the wavelength. In most experimental situations, electrode size and electrolyte conductivity will be larger than in the estimation given here, so that the threshold is still 1 to 2 orders of magnitude larger. Even when taking into account that the kinetic constants kad and kr can vary largely, it is evident that inequality (30) will be fulfilled for many experimental situations. This can be traced back to the fact that the characteristic length of migration, that is, of the inhibitor transport process, is much larger than the one of diffusion, that is, of the activator transport process. Therefore, Turing structures will form more easily in electrochemical S-NDR systems than in chemical systems. Next, let us compare how likely it is to observe oscillations compared to Turing structures in S-NDR systems. For oscillations to occur, it is necessary that the time scale of the activator, which is in our case a chemical species, is larger than the time scale of the inhibitor, the double-layer potential. This constraint is not very likely to be met, as can already be seen from a comparison with N-NDR systems. Due to the interchanged roles of potential and chemical variables, this condition is again just opposite to that in N-NDR systems, where the time scales were shown to be mostly favorable for oscillatory behavior.8,9 In terms of our model parameters, oscillations can be expected when µ is at least of the order of κ. As we have seen, such a restriction on the time scales does not exist for the formation of Turing patterns. The latter also occurs when µ < κ, making in S-NDR systems the experimental occurrence of a Turing bifurcation much more probable than the observation of oscillations. Stationary patterns in electrochemical systems were observed during the reduction of persulfate, a well-known N-NDR system.35 It is important to distinguish between the mechanism giving rise to those structures and the Turing patterns discussed in this paper. In the persulfate experiment, a Haber Luggin capillary was used, an electrode arrangement that, in connection with the potentiostatic operation mode of the system, induces a global coupling.35-37 The global coupling, in turn, destabilizes the homogeneous stationary state. There are two essential differences between the two types of instabilities. First, from a dynamic point of view, a Turing structure possesses an intrinsic wavelength, independent of the system size. In contrast, a destabilizing global coupling leads to the formation of two spatial domains, that is, the wavelength of the pattern is identical to the system’s length. Second, from a physical point of view, the transport processes involved in the destabilization are qualitatively different: The transport process associated with the Turing instability is migration of ions in the electrolyte. In contrast, the global coupling results from the external circuitry; thus, if one wants to identify at all the corresponding transport process here, it is the motion of electrons in the metal. At last, consider again function f (eq 5), which describes adsorption and desorption of an adsorbate with attractive interactions. One should be aware that it is a drastic simplification to couple this function with eq 2, where the communication of the adsorbate species over the surface is described by a simple diffusion term. Rather, the “diffusion constant” is proportional to ∂2A/∂u2, where A is the free energy density of the uniform adsorbate. This quantity will be negative for a uniform steady state within the spinodal envelop,38 opening another mechanism that can destabilize the uniform steady state. The corresponding patterns are typically of nanometer scale.39,40 Here, we disregarded this effect because it is our aim to demonstrate under which conditions a migration-driven instability can occur in

6090 J. Phys. Chem. B, Vol. 104, No. 25, 2000 S-NDR systems. In a system in which a phase transition of the adsorbate layer plays a role, two mechanisms inducing an instability of the homogeneous system will be present, and consequently the spatio-temporal dynamics is more complicated than in both subsystems. This fact alone already suggests that the influence of each of the two mechanisms be first investigated separately. Concluding Remarks S-NDR systems represent a further class of electrochemical systems in which the interaction between a negative differential resistance caused by the electrode kinetics and the electric control of the experiment brings about spatial and temporal instabilities. Many of them are not obvious at all from measurements of averaged quantities, such as the current density. This underlines the necessity to take into account the possibility of a spatially inhomogeneous interface when interpreting experimental data. A prerequisite, therefore, is a thorough understanding of pattern formation in electrochemical systems, which has only begun to be developed. Thus, such studies are not only fruitful in the context of nonlinear dynamics but can also help to understand the interfacial behavior in different situations. Acknowledgment. This work is part of project B4 of the Sfb 555, which is supported by the Deutsche Forschungsgemeinschaft. References and Notes (1) Turing, A. M. Philos. Trans. R. Soc. London, Ser. B 1952, 237, 37. (2) Murray, J. D. Mathematical Biology; Springer: Berlin, 1990. (3) An activator acts autocatalytic for its own production and also produces the inhibitor. The inhibitor, on the other hand, consumes the activator and decays with a certain rate. (4) Field, R. J.; Burger, M. Oscillations and TraVeling WaVes in Chemical Systems; Wiley: New York, 1985. (5) Epstein, I. R.; Pojman, J. A. An Introduction to Nonlinear Chemical Dynamics; Oxford University Press: New York, 1998. (6) Castets, V.; Dulos, E.; Boissonade, J.; De Kepper, P. Phys. ReV. Lett. 1990, 64, 2953. (7) Boissonade, J.; Dulos, E.; De Kepper, P. In Chemical WaVes and Patterns; Kapral, R., Showalter, K., Eds.; Kluwer Academic: Dodrecht, 1995. (8) Koper, M. T. M. In AdVances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; Wiley: New York, 1996; Vol. 92. (9) Krischer, K. In Modern Aspects of Electrochemistry; Conway, B. E., Bockris, J. O., White, R., Eds.; Kluwer Academic/Plenum: New York, 1999; Vol. 32.

Mazouz and Krischer (10) Fla¨tgen, G.; Krischer, K. J. Chem. Phys. 1995, 103, 5428. (11) Mazouz, N.; Fla¨tgen, G.; Krischer, K. Phys. ReV. E 1997, 55, 2260. (12) Fla¨tgen, G.; Krischer, K. Phys. ReV. E 1995, 51, 3997. (13) Fla¨tgen, G.; Krischer, K.; Ertl, G. Z. Naturforsch. 1995, 50a, 1097. (14) Otterstedt, R.; Plath, P. J.; Jaeger, N. I.; Sayer, J. C.; Hudson, J. L. Chem. Eng. Sci. 1996, 51, 1747. (15) Otterstedt, R.; Plath, P. J.; Jaeger, N. I.; Hudson, J. L. Phys. ReV. E 1996, 54, 3744. (16) An example of a bistable polarization curve which does not fall into this category is a system with a Z-shaped polarization curve. (17) Lejay, M. M. E.; Wiart, R. C. R. Acad. Sci. Paris 1973, 277C, 833. (18) Epelboin, I.; Ksouri, M.; Lejay, E.; Wiart, R. Electrochim. Acta 1975, 20, 603. (19) It is assumed that the reference electrode is at the height of the counter electrode, which provides an equipotential plane. Thus, φ(z)0))0 is equivalent to setting the potential at the equipotential plane to 0. (20) We consider here only systems whose “chemical part” can be expressed mathematically by one variable. (21) Koper, M. T. M. Electrochim. Acta 1992, 37, 1771. (22) Roelfs, B.; Bunge, E.; Schroeter, C.; Solomun, T.; Meyer, H.; Nichols, R. J.; Baumga¨rtel, H. J. Phys. Chem. B 1997, 101, 754. (23) Ho¨lzle, M. H.; Wandlowski, T.; Kolb, D. M. Surf. Sci. 1995, 335, 281. (24) Dretschkow, T.; Dakkouri, A. S.; Wandlowski, T. Langmuir 1997, 13, 2843. (25) Striegler, H. Ph.D. Thesis, Universita¨t Ulm, 1998. (26) Mazouz, N. Ph.D. Thesis, FU Berlin, 1999. (27) Frumkin, A. N.; Damaskin, B. B. In Modern Aspects of Electrochemistry; Bockris, J. O.; Conway, B. E., Eds.; Buttlerworths: London, 1964. (28) Mazouz, N.; Krischer, K.; Fla¨tgen, G.; Ertl, G. J. Phys. Chem. 1997, 101, 2403. (29) Mazouz, N.; Fla¨tgen, G.; Krischer, K.; Kevrekidis, I. G. J. Electrochem. Soc. 1998, 145, 2404. (30) Hindmarsh, A. C. ACM-SIGNUM Newsl. 1980, 15, 10, (lsode is available from the netlib library ([email protected])). (31) Doedel, E. J.; Kevernez, J. P. AUTO: A program for continuation and bifurcation problems in ordinary differential equations; California Institute of Technology, Applied Mathematics Report, 1986. (32) deWit, A. In AdVances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; Wiley: New York, 1999; Vol. 109. (33) Note that in our model a change of L does not change the homogeneous solution. This is different when two-dimensional electrode geometries are considered [34]. (34) Christoph, J. Ph.D. Thesis, FU Berlin, 1999. (35) Grauel, P.; Christoph, J.; Fla¨tgen, G.; Krischer, K. J. Phys. Chem. B 1998, 102, 10264. (36) Global coupling means that a local change in the state of the system affects instantaneously the whole spatial domain with the same strength. (37) Christoph, J.; Otterstedt, R. D.; Eiswirth, M.; Jaeger, N. I.; Hudson, J. L. J. Chem. Phys. 1999, 110, 8614. (38) Kevrekidis, I. G.; Aris, R.; Schmidt, L. D. In I. Chem. E. Symp. Series 87; Pergamon: Elmsford, NY, 1984. (39) Hildebrand, M.; Mikhailov, A. J. Phys. Chem. 1996, 100, 19089. (40) Hildebrand, M.; Mikhailov, A.; Ertl, G. Phys. ReV. E 1998, 68, 5483.