Nucleation in Systems with Multiple Stationary States BYA. NITZAN,P. ORTOLEVA AND J. Ross Department of Chemistry Massachusetts Institute of Technology Cambridge Massachusetts 021 39 Received 26th July 1974 We consider a reaction diffusion system far from equilibrium which has multiple stationary states (phases) for given ranges of external constraints. If two stable phases are put in contact then in general one phase annihilates the other and in that process there occurs a single front propagation (soliton). We investigate the macroscopic dynamics of the front structure and velocity for two model systems analytically and numerically and for general reaction-diffusion systems by a suitable perturb- ation method. The vanishing of the soliton velocity establishes the analogue of the Maxwell con- struction used in equilibrium thermodynamics.The problem of nucleation of one phase imbedded in another is studied by a stochastic theory. We show that if the reaction dynamics is derived from a generalized potential function then the macroscopic steady states are extrema of the probability distribution. We use this result to obtain an expression for the critical radius of a nucleating phase and confirm the prediction of the stochastic theory by numerical solution of the deterministic macro- scopic kinetics for a model system. 1. INTRODUCTION Chemical reaction mechanisms with macroscopic rate laws of sufficient non- linearity in systems maintained far from equilibrium may have multiple stationary stable states for given external constraints.We refer to each such state as a phase and transitions between phases are possible. The analogy of the theory of phase transitions and critical phenomena to transitions between stable states and critical (marginal stability) points has been discussed in some detail both at the macro- scopic and statistical mechanical level of anal~sis.~-~ In this paper we investigate the nucleation of one phase within another phase as well as the conditions of coexistence of two phases (the analogue of the Maxwell construction). These problems have been considered by Kobatake who showed for a particular case the similarity of the behaviour of the generalized entropy production in transitions between stable branches of steady states and the Gibbs free energy in equilibrium phase transitions ; by Schlogl,* who took into account reaction and diffusion and treated coexistence of phases and the analogue of the Maxwell construction for one variable systems with the help of a mechanical analogy; and by Nicolis Malek-Mansour van Nypzlseer and Kitahara,lO who analyzed the onset of instability as a nucleation process derived a ncn-linear master equation for that purpose and applied that result to some ex-amples.We approach these problems from two diff'erent points of view. First in section 2 we investigate the macroscopic dynamics of two stable stationary states of semi-infinite extent at given external constraints placed in contact with each other. Except for one value of the constraints one phase annihilates the other and in that process there occurs a single-front propagating wave a soliton.We investigate the behaviour of the concentration profile and velocity of such a wave for two model systems analytically and numerically and show that the values of the external con- straints for zero soliton velocity establish the Maxwell construction. 241 242 NUCLEATION Further we study these quantities for general reaction-diffusion systems by a suitable perturbation method and find that the soliton velocity vanishes linearly in the deviation of the external constraint from its value at the Maxwell construction. Next in section 3 we discuss some aspects of a stochastic theory of instability phenomena on the basis of a Fokker-Planck equation assumed for a reaction-diffusion system.That equation can be solved for the steady state probability distri- bution for which the reaction dynamics is obtainable from a generalized potential function. In that case we show that the macroscopic steady states are extrema of the probability distribution maxima for stable steady states. We apply this theory to analyze the behaviour of a system consisting initially of a nucleus of one phase embed- ded in an infinite bulk of the other phase ; we do so in a manner similar to the conven- tional treatment of nucleation in first order phase transitions.". l2 We thus obtain an expression for the critical radius for such a nucleus above which the nucleus grows and below which it disappears in time. This prediction agrees with numerical solutions of the deterministic diffusion-reaction equation for a simple model system.2. MACROSCOPIC DYNAMICAL THEORY We consider a diffusion-reaction system for which the niacroscopic deterministic equation of motion is The symbols denote $ a column vector of concentrations and possibly other state variables such as temperature; Y spatial coordinates; t time; B a matrix of diffusion coefficients ; F,the variations due to chemical reaction ; and X the set of external constraints (boundary concentrations of species light intensity in an illurn-inated ~ystern,~ etc.). We assune that 1. the diffusion matrix is constant and sym- metric and the diffusion process is stable in that in the absence of reaction eqn (2.1) with F = 0) the steady state solution is stable; 2.P is analytic in 9 and 1. We consider systems for which there exist two stable stationary states (and one unstable stationary state) in a given range of k. Thus if in a one-dimensional system (co- ordinate x) we prepare two systems under identical external constraints one system in one stable state (one phase) and the other system in the other stable state (the other phase) then the boundary at x = 0 between the two systems will move in one direction or another depending on X. There occurs a single front propagation,13 a soliton in which one phase annihilates the other. We shall find the condition for coexistence of two such systems (phases) to be zero soliton velocity or a standing front at a given X. We therefore seek wave solutions of (2.1) which are most readily studied in a reference frame moving with the front.W7edefine the phase 4 = x-ct (2.2) where the soliton velocity o is a function of k. With (2.2) the partial differential eqn (2.1) is converted into an ordinary differential equation which is a nonlinear eigenvalue equation for u at a given X. Coexistence occurs when u = 0 at h = AM. The value of hMis obtained from eqn (2.4) which may be inrer-preted as an eigenvalue equation for AM (2.4) A. NITZAN P. ORTOLEVA AND J. ROSS We first consider a soluble model system and then present an analysis for more general {D F}systems. The model system consists of a single species with third-order kinetics which obey the equation a2t,b all/-= D 7-q(t,b-a)(t,b-b)(ll/-c).at dx-Homogeneous steady states occur at $* = a,b c where these parameters are arbitrary functions of 2.. We take b 2 c > a q > 0,and thus only states a and b are stable to small homogeneous perturbations. If we define a new variable u such that *-a UE- 6-a' then we obtain dU a2u -= D-2-q(b-a)2u(u-l) at ax This is identical in form to a model considered by Montro11.12 In terms of our notation his solutions for soliton fronts are = -(b-a). (2y The front velocity is given by v = (+qD)*(a+b-2~). (2.9) Coexistence occurs at u = 0 i.e. c = +(a+b) at 2 = AM (2.10) which is obvious from the symmetry of the trilinear kinetics for which the two stable states a and b are located symmetrically about the unstable state c.It is interesting to take a particular case and analyze it further. Let F = -q@($2-2$ + 1 /A). (2.1 1) Then v = (+qD)3[3(1-1 /A)+-I] and = (q/2D)+[1+(1-1 /A)+] from(2.8,9). The uniform steady states and the soliton velocity are shown in fig. I(a) and (b). Co-existence is found at AM = 918. With this we obtain the structure of the coexistence region to be $coex .( x = $( +exp(kPMMX))-(2.12) where & = 3(8q/D)+;the structure is shown in fig. 2. We note from fig. l(b) that the velocity passes through the origin at AM with a finite but non-zero slope l(du/dL)AMl# co. This general property provides a useful aid in determining AM experimentally. Next we consider the general system (2.1) by analyzing the solutions of the front eqn (2.3)in the vicinity of AM.To do so we introduce a parameter A to be defined at a later stage in the development which measures the deviation of the system from coexistence conditions. In general A is a function of A 2 = &A" (2.13) n=O 244 NUCLEATION 0 5 10 A FIG.1 .-Steady states and soliton velocity versus the parameter A for the model system eqn (2.1 1). and we expect that A approaches zero proportional to the factor (A-AM) some positive power. Furthermore 1,6 and v depend on A raised to 03 (2.14) rm v = C vnAn. (2.15) n=O Such expansions have been employed elsewhere for the study of chemical waves.lS On substitution of the expansions (2.13-15) into (2.3) we recover to zero'th order in A FIG.2.4tructure of the coexistence region versus distance eqn (2.12) for the model system eqn (2.11).A. NITZAN P. ORTOLEVA AND J. ROSS 24s eqn (2.4) with the boundary conditions of the system being one stable stationary state at x = -a,and the other stationary stable state at x = +a. We also have by definition in this order uo = u(R& = 0. To first order in the parameter A we find (2.16) with the operator 2’ defined by d2 9 = D-+R (2.17a) d242 (aF/a$) (2.17b) and the subscript 0 following differentiation implies evaluation in lowest order. We note that differentiation of (2.4) with respect to 4 yields (2.18) We assume that d$/d@ is the only eigenfunction of 9 with zero eigenvalue. P-l(d$/d#) does not exist and hence the solution of (2.8) for t+h1cannot be found for arbitrary (vl R1).The problem is overcome by a proper choice of these coefficients. We then define a vector 0 such that 9+0= 0 (2.19) where P+is the adjoint operator of 9.If we define an inner product (LIP)= c d4 44 (2.20) 1 -53 where the sum on i runs over all species then from (2.17,20) L?+ is given by d2 9+= D+-+a+. (2.21) d+2 (By assumption we have D+ = Dt Onsager’s relations). We can be assured of the solubility of the first order eqn (2.16) if we can choose v1 and A1 such that the r.h.s. of (2.16) is orthogonal to the zero eigenvalue function d$,/d4. Thus we obtain by use of the adjoint eigenfunction to d$o/d4 that is 0,see (2.19) that u,(@I(d$o/d4)) +M@l(aFlwo)= 0. (2.22) The coefficient of v may be chosen to be unity by proper normalization of 0.Hence if the coefficient of A1 is finite then this relation can be satisfied for finite vl R1. Finally (aF/d;l) approaches constant values as 1&-+ co. A necessary condition for al solubility is 1 d$lO(4)l < 00 in order to find R1 vl # 0. For systems where -0 R+ = i-2 9is self adjoint and therefore 0= d$o/d&. Since d$o/d$ is localized to the region of the soliton front it satisfies the necessary condition of solubility. However we have not investigated this condition for systems with R+# s1. It is sufficient to all orders to choose Rn32 = 0 ; thus with the convenient choice R1 = 1 we have A R-RM. (2.23) To first order the perturbation solution is b Y-Y 246 NUCLEATION + ' VI4M -~(ol(aF/a~)M)/(ol(d~O/d#))](~-~~) ' (2.25) For the one variable system (2.5-12) we have verified these relations using the equation 0 = dt,bo/d4.I I --4.0 5.0 6.0 So FIG.3.-Steady state concentration and soliton velocity versus Sofor the model system described in Appendix A. For curve (a) the diffusion coefficients Ds Dp are both 1 and for curve (b)DS = 1 and Dp = 0.05. 0.0 8.0 16.0 X FIG.4.-Soliton structure against distance at different times for the model system Appendix A. The diffusion coefficients are Ds= 1 Dp = 0.05. The external concentration S" = 5. A. NITZAN P. ORTOLEVA AND J. ROSS Finally we note that the structure and velocity of the soliton depends in general not only on the reaction term P but also on the diffusion matrix D.If D is equal and diagonal then it may always be eliminated from the equation for the soliton front by a suitable scaling of the length 4-+(JD)4'. However for more complex diffusion matrices this simple scaling no longer holds. In Appendix A we consider a two species system with rate mechanism & S"$S+P+P" (2.26) where So,Poare maintained constant and the activity of the enzyme 8depends on P. In fig. 3 we show the steady state concentration P* as a function of s"and the results of the numerical integration of the partial differential equations yielding the variation of the soliton velocity against the concentration of s" for various choices of diffusion coefficients of S and P. For this system the value of Soat coexistence (AM of the general theory) does not appear to depend on the choice of diffusion coefficients to within the accuracy of the numerical computations; the results of the stochastic theory section 3 indicate no dependence on that choice for systems for which the kinetics is derivable from a generalized potential function.In fig. 4 we show the soliton stucture for the same system (Appendix A) for the choice of diffusion coefficients Ds = 1 Dp = 0.05 and So = 5. 3. STOCHASTIC THEORY Our treatment is based on a model in which to the set of reaction diffusion rate equations we add a stochastic source term to account for fluctuations in the system. The set of concentrations at each spatial point becomes then our set of stochastic variables. It is convenient to make this set discrete by dividing the system into cells of homogeneous composition.We thus consider the system of volume V to be partitioned into N cells each of volume ij = V/N. The macroscopic reaction kinetics is again given by with a similar equation for $n in the nth cell. The rate of change of numbers of particles of any species x, = fiq,,,due to diffusion is where the sum extends over the 2d nearest neighbour cells (d is the dimensionality of the system). The factor of Ui-2/d comes from multiplying by the surface area ijl-l'd of the cell and dividing by the characteristic length of a cell such that ($n+i-$,,)/fi*/d is an approximation to the concentration gradient. Now we write a Langevin equation for this system as (3.3a) or (3.3b) For simplicity we take the stochastic term to be Gaussian and thus to satisfy the usual relations (A,(?))= 0and (&,,(~')h,Jt)) = 6,p6,,,$(t-t') (a,/3 denote different chemical NUCLEATION components).Eqn (3.3b)is equivalent to the following Fokker-Planck equation for the probability density of finding the concentration distribution $,, P[Y], in the system where Mn(Y](with components MaJ3!]) is given by 2d and where Y denotes the “ supervector ” *n (n = 1 . . . N) (= (It = 1 . . . N a goes over the different chemical species)). In (3.4) K was taken for simplicity as a scalar matrix K = KI where I is the unit matrix. The integration of this many- variable Fokker-Planck equation for steady state conditions (Plat) = 0 is possible provided the reaction kinetics are such that the curl of’M is zero,I6 aMan -aMpn (3.6) 89s wa This condition should be modified if D is not a symmetric matrix or if K is not a scalar matrix.If (3.6) holds then we obtain P[W] = Cexp(-U[Y]) (3.7) where C is a normalization constant. For simplicity we assume now that we have a one-component system for which (3.6)always holds. The results for a more general case with a symmetric diffusion matrix are given in Appendix C. We write U[+]as which in the continuum limit becomes The function G[$(v)]is defined by (3.10) which in case of many chemical species should be understood as a line integral in $ space. Condition (3.6) ensures that this line integral is well defined and does not depend on the integration path.It is interesting to note that (3.9) has the same form as the free energy in the Landau-Ginzburg theory (there G($) is usually taken as a quartic form). Also critical phenomena have been discussed with Hamiltonians of that form.17 Having arrived at an expression for the steady state probability distribution P[$] we show next that extrema in that distribution are given by the steady states of the deterministic equation. Differentiation of (3.7) yields (3.11) A. NITZAN P. ORTOLEVA AND J. ROSS 249 or hp= dr{F[$]G$-DV$ VS$) (3.12) which on partial integration becomes K dr(F[1,b]+DV~$}6$. (3.13) For P to be an extremum for the arbitrary variation S$ we must have the integrand vanish ; that is the macroscopic condition for steady states.It can be further shown that stable steady states correspond to maxima of the probability distributions (see Appendix B). The system may be in a single stable stationary state (although there exist other such states at the same external constraints) and in that case the system is homogeneous i.e. F[$] = 0. However two (or more) stable stationary states may coexist and the distribution of concentrations is the solution of DV'++F[$] = 0 for given boundary conditions. Note that this is exactly the equation for the zero velocity soliton structure (eqn (2.3) with u = 0). The stochastic theory for the probability distribution of a system in a steady state may be applied to the problem of nucleation. Consider a spherical nucleus of radius R of one stable stationary state labelled phase B immersed in an infinite bulk of another stable stationary state labelled phase A both under the same external constraints.Thus we have F[@J= F[i,bB]= 0 for the individual homogeneous phase. For phase B immersed in A however we may inquire about the expected stability of the spherical nucleus that is the nucleus is expected to evolve in the direc- tion of increasing the probability of the overall structure (nucleus +bulk). The Iogarithm of the ratio of the probability distribution for B phase immersed in A phase, PBA,to that of pure A phase PAis where I is the width of the interface region and where we have approximated the gradient by ($B-$A)/Z. In eqn (3.14) we have used the following expression for the surface area S and the volume V of a d-dimensional sphere (ris the r-function) The derivative (with respect to R) of the r.h.s.of (3.14),set equal to zero gives the radius of the sphere of phase B for which the 1.h.s. of (3.14) is a minimum that is (for a given I) (3.15) If G[#,] > G[@*] then there exists a distribution (characterized by R,) of maximum "free energy " (minimum probability) such that for R > R the nucleus will grow (as this is the direction of increasing the probability) and for R < R it will disappear in time. If G[I,~~] < G[$A] then no positive R exists. Finally if G[$J = G[t,bA]then R = a,and phases A and B coexist with a planar (zero curvature) surface of contact. This condition is clearly the analogue of the Maxwell construction in equilibrium 250 NUCLEATION phase transitions.Furthermore we expect that this condition is equivalent to the one of zero soliton velocity and have confirmed the expectation for the example given by (2.5). In a many component curl-free system with a symmetric non-scalar diffusion matrix the result (see eqn (C.l) Appendix C) shows that the Maxwell construction is again given by the condition G[31g]= G[$J. We stress that the arguments presented here are based on probabilistic consider- ations especially on the postulate that a system will evolve in the direction of increasing the probability. The result (3.15) is similar in many respects to that for the radius of critical nucleus in equilibrium first order phase transitions.' * l2 The contribution of the diffusion process to the " free energy " U eqn (3.9) takes here the place of the surface tension in the equilibrium nucleation treatment.The appearance of the parameter I io eqn (3.15) takes into account the non-zero thickness of the interphase layer. It should be remembered that I is not an independent parameter but is deter- mined by the external parameters. Also the simplifying assumption of linear varia- tion of the concentration t,b across the interphase layer limits the applicability of the result (3.15). We compare the results of the stochastic theory to the solution of the deterministic dynamical equations. To this end we numerically integrated the one variable partial differential equation !?!! = DV2t+b-J13+bt,b+~ (3.16) at in which the cubic polynomial F($) = -G3+b$+c plays the role of a chemical source term.For b > 0 and Jcl < (2b/3)(,/b/3) this equation has two stable and one I .€ 4 O.( -1.t 1 t=O I J 0.0 4.5 91 1.0 4.5 9.0 r r FIG.5.-Evolution of a nucleus as a function of radial coordinate for various times as obtained from numerical solution of eqn (3.16) with b = 3 c = 1.S and D = 1. Figures at left and right are for initial radius smaller and larger than & respectively. A. NITZAN P. ORTOLEVA AND J. ROSS 25 1 unstable homogeneous steady states. Denoting the two stable homogeneous steady states by (phase A) and t,kB (phase B) we put our initial condition to be a sphere of radius R centred around the origin and consisting of phase B imbedded in a spherical region of phase A which in the computation process has to be taken finite but with radius R1much larger than R.Since the initial condition is spherically symmetric and the coefficients of (3.16) are constants it is sufficient to consider the radial coordinate only. The boundary conditions are taken as aIrl//ar = 0 at r = 0 and at the larger bounding sphere. Fig. 5 shows the time and space dependent solution to eqn (3.16) in three dimen- sions for the choice of coefficients D = 1 b = 3 c = 1.5 and for two different choices of the initial radius R R < R and R > R,. The nucleus shrinks and disappears for R < R,and grows in time for R > R,. It is interesting to note that the growth process (for R > R,)accelerates for small R but the velocity approaches a constant for larger radii as the front becomes more planar.We calculate the velocity of a planar front (from eqn (2.9) using the steady state solutions t,ka = -1.38 t,hb = 1.94 (stable) and Irl/c = -0.558 (unstable)) to be 1.2 while from fig. 5b the velocity for the largest radius obtained is 1.0. TABLECRITICAL RADIUS FOR NUCLEATION FOR THE SYSTEM GIVEN BY EQN (3.16) Rccalculated from eqn (3.15) Rc from numerical computations with 1 = 2 -C 2 dimensions 3 dimensions 2 dimensions 3 dimensions 0.50 2.6 < R < 3.2 3.4 0.75 3.5 < R < 5.5 4.6 1.oo 1.5 < R < 1.75 2.7 < R < 3.7 1.7 3.4 1.25 2.1 c R < 2.3 2.6 1SO 0.85 < R < 1.0 1.65 < R < 1.95 1.05 2.1 1.75 0.8 < R < 1.0 1.35 -c R < 1.5 0.9 1.75 1.875 1.2 < R < 1.4 1.6 2.0 R < 0.5 R < 0.8 0.7 1.35 In table 1 we present the results of numerical computation of R,.as well as results based on eqn (3.19 for b = 3 and for different values of c in two and three dimensions.The calculations based on eqn (3.15) were done with I = 2 which is an estimate based on the numerical results such as given by fig. 5b. The agreement between the values of R,as obtained from (3.15) and from numerical integration of (3.16) is fairly good for all the values of the parameter c that were tried excluding the bifurcation value c = 2. (Our failure to provide a more accurate numerical answer for R near the bifurcation point is caused by the slow time evolution of the system near that point which makes the numerical integration very expensive.) This agreement is remark- able as the two ways of obtaining R,are different; one is based on solutions of the time-dependent deterministic equation and the other on the steady state probability distribution.To conclude this section we should note that even though the dynamics of nucleus growth as given by eqn (3.16) or fig. 5 give important information on the rate of nucleation of one phase within another the rate of this process depends also on the rate of formation of nuclei by spontaneous fluctuations which is not discussed in the present work. We thank John M. Deutch for helpful discussions. This work was supported in part by the National Science Foundation and Project SQUID Office of Naval Research.252 NUCLEATION APPENDIX A For the reaction mechanism (2.26) the reaction-diffusion equations are as -= 82s -2D +(SO-S)-G(P)S at ax ap a2p --D -+(Po-P)+d(P)S at ax2 with all rate coefficients taken to be unity. The assumed form for the dependence of 8 on P is P2 (A.3) &(P) = 1+P+P2 for which chemical examples are available.'* The steady state solutions of (A.l,2) are P* = 0; t(So-1)+$[(So-1)2-8J* 64.4) which are shown in fig. 3a. APPENDIX B Here we show for the probability distribution given by eqn (3.7) and (3.9) that stationary points which correspond to stable macroscopic steady states are local maxima of the distribution or minima of the potential U[$] given by (see eqn (C.l)) U[$l = -1lK 1dr{G[$(r)l-4D :(V$m2}.(B.1) Let the steady state variable $*(r) satisfy the equation F[$*(r)]+DV2$*(r) = 0 (B.2) and let $(r) = $*(r)+a$(r). We assume that \II*(r)corresponds to a stable steady state of the system. The first order (in S$) term of U then vanishes and the second order term is given by a2u = -3 JdrslC/(n+ Dv2)S$ where i2 = (iW/a$)+*. We have used the identity JdrVWeV V= -J drWV2V+1 WneVVdo (B.4) and the fact that $*(r) and \II(r)must satisfy the same boundary conditions and there- fore = 0. In eqn (B.4) E is the boundary of the system do is an element of this boundary and n is a unit vector perpendicular to the boundary. Since we assume that $*(r) is a stable steady state the operator +DV2 must be negative definite and therefore d2 U in eqn (B.3) is positive.This establishes the fact that U[$] is a mini-mum (and P[$]is a maximum) for $ = $*. APPENDIX C Here we write the modified relations of section 3 for a many variable system and for a general symmetric diffusion matrix. The potential U[$] (eqn (3.9)) is given by wwl = -1 lq dww91-3D (V$(rN2} (C.1) A. NITZAN P. ORTOLEVA AND J. ROSS where The expression for the critical radius R is now replaced by N. Rashevsky Mathematical Biophysics (Dover New York 1960) Vol. I. B. B. Edelstein,J. Theoret. Bid 1970 29 57. P. GIansdorfT and I. Prigogine Thermodvnamic Theory of Structure Stability and Fluctuations (Wiley-Interscience New York,1971). G. Nicolis Adv. Chem. Phys.1971 19 209. A. Nitzan and J. Ross J. Chem. Phys. 1973 59 241. H. Hahn P. Ortoleva and J. Ross J. Theoret. Biol. 1973 41 503. ’Y. Kobatake Physica 1970,48 301. F. Schlogl 2.Phys. 1972 253 147. A. Nitzan P. Ortoleva J. Deutch and J. Ross J. Chem. Phys. 1974,61 1056 see also other references therein. lo G. Nicolis M. Malek-Mansour A. van Nypolseer and K. Kitahara private communication (preprint). L. D. Landau and E. M. Lifshitz Statistical Physics (Addison-Wesley Reading Massachusetts 1958.) l2 see for example B. Widom in Phase Transitions and Critical Phenomena Vol. 2 ed. C.Domb and M. S.Green (Academic Press New York 1972). l3 R.Davidson Methods in Nonlinear Plasma Theory (Academic Press New York 1972). l4 E. W. Montroll in Statistical Mechanics ed.S. A. Rice K. F. Freed and J. C. Light (Univ. of Chicago Press 1972) p. 69. P. Ortoleva and J. Ross J. Chem. Phys. 1974 60,5090. l6 R. L. Stratonovich Topics in the Theory of Random Noise (Gordon and Breach New York 1963). see for example S. Ma Rev. Mod. Phys. 1973 45 589. l8 H. Hahn A. Nitzan P. Ortoleva and J. Ross Pruc. Nat. Acad. Sci. (submitted).