Phase ordering with a global conservation law: Ostwald ripening and coalescence Massimo Conti1, Baruch Meerson2, Avner Peleg2 and Pavel V. Sasorov3 1 Dipartimento di Matematica e Fisica, Universit`a di Camerino, and Istituto Nazionale di Fisica della Materia, 62032, Camerino, Italy 2The Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel 2 3Institute for Theoretical and Experimental Physics, Moscow, 117259, Russia 0 0 Globally conserved phase ordering dynamics is investigated in systems with short range corre- 2 lations at t = 0. A Ginzburg-Landau equation with a global conservation law is employed as the phase field model. The conditions are found underwhich the sharp-interface limit of this equation n isreducibletothearea-preservingmotion bycurvature. Numericalsimulations showthat,for both a J critical and off-critical quench, the equal time pair correlation function exhibits dynamic scaling, and the characteristic coarsening length obeys l(t) t1/2. For the critical quench, our results are 3 ∼ inexcellentagreementwithearlierresults. Foroff-criticalquench(Ostwaldripening)weinvestigate ] thedynamicsofthesizedistributionfunctionoftheminorityphasedomains. Thesimulationsshow h that,atlargetimes,thisdistributionfunctionhasaself-similarformwithgrowthexponent1/2. The c scaled distribution, however, strongly differs from the classical Wagner distribution. We attribute e thisdifferencetocoalescenceofdomains. AnewtheoryofOstwaldripeningisdevelopedthattakes m intoaccountbinarycoalescenceevents. Thetheoreticalscaleddistributionfunctionagreeswellwith - that obtained in the simulations. t a t s I. INTRODUCTION . t a m Phase ordering is emergence of order from disorder through domain growth and coarsening. The standard setting - when phase ordering occurs is a temperature quench from a high-temperature disordered phase into a two-phase or d a multiphase region. Phase ordering has been the subject of extensive research during the last two decades [1]. An n importantsimplifyingassumptioninphaseorderingtheoryisdynamicscaleinvariance. Accordingtothisassumption, o the coarsening system possesses, at late times, a single relevant dynamic length scale l(t) (the characteristic domain c [ size)whichgrowswithtimeasl(t) tα [1]. Itisbynowwellestablishedthatinsystemswithshortrangecorrelations α=1/2 for non-conserved (model∼A) dynamics, while α=1/3 for locally conserved (model B) dynamics . 2 Thereis,however,animportantadditionalcoarseningmechanism: globallyconserved phaseordering[2,3,4,5,6,7]. v Globally conserved dynamics can be thought of as model A dynamics constrained by global conservation of the 7 5 order parameter: for example, Ising model with fixed magnetization. This global conservation law is maintained 1 by an external field (for example, a magnetic field) which depends on time but is uniform in space. The globally- 9 conserved phase ordering is accessible in experiment. Consider the sublimation/deposition dynamics of a solid and 0 its vapor in a small closed vessel kept at a constant temperature below the melting point. As the acoustic time 1 scale in the gas phase is shortcompared to the coarseningtime, the gas pressure (and, consequently, density) remain 0 uniform in space, changing only in time. This character of mass transport in the vapor phase makes the coarsening / t dynamics conserved globally rather than locally. An important characteristics of globally-conserved dynamics is a m interface-controlledkinetics, incontrastto the bulk-diffusion-controlledkinetics typicalfor locally-conservedsystems. Interface-controlledkineticswasinvestigatedinthecontextofgrowthofsmallplatinumparticlessupportedonalumina - d substrates in an oxidizing environment [8]. There are additional examples of cluster growth on surfaces [9], where n it was found possible to single out the interface-controlled kinetics [10]. There is also strong evidence in favor of o globally conserved interface-controlled transport during the coarsening of clusters in granular powders driven by a c low-frequency electric field [11, 12]. : v Part of the theoretical importance of globally conserved phase ordering lies in the fact that it enables an access to i off-critical quenches in the simpler model A dynamics (with global conservation). Thus, it allows one to determine X which characteristics of the system depend on the volume (or area) fraction ε and which do not. r a Dynamic renormalization group arguments show that global conservation should not change the growth law [13]. This early result was confirmed by particle simulations with short range correlations in the initial conditions: for critical (ε = 1/2) [6, 14] and for off-critical (ε < 1/2) [4] quench. Recent phase field simulations of systems with long-range (power-law) correlations in the initial conditions have also shown dynamic scale invariance with the same growthexponentα=1/2[7]. Therefore,α=1/2independentlyofε. Ontheotherhand,theautocorrelationfunction [4, 15] and persistence exponent [15] were found to be ε-dependent. Globally conserveddynamics are relatedto a wide range of multiphase coarseningsystems. Sire andMajumdar [4] showedthatinthelarge-qlimitthedynamicsoftheq-statePottsmodelareequivalenttothedynamicsoftheglobally conserved model with an area fraction ε= 1/q. The large-q limit of the Potts model is of practical importance as it describescorrectlysomeofthe dynamic characteristicsofdrysoapfroths[16]andofthe coarseningofpolycrystalline 2 materials [17]. In the limit of a vanishing volume fraction of the minority phase the late stage of coarsening is describable by the mean field theories of Ostwald ripening. Lifshitz and Slyozov [18] developed such a theory for the bulk-diffusion- controlled (or locally conserved) dynamics. They showed that the size distribution function of the minority domains approaches,atlargetimes,aself-similarform. Correspondingly,theaveragesizeoftheminorityphasedomainsgrows with time like t1/3. Following the seminal work by Lifshitz and Slyozov [18], Wagner developed a similar mean-field theory for the interface-controlled (globally conserved) Ostwald ripening [19]. The Wagner’s theory yields a growth lawt1/2 fortheaveragedomainsize,andadifferent(broader)shapeofthescaleddistributionfunction. Weshallrefer to this scaled distribution function as the Wagner distribution. More recently, it was shown that interface-controlled Ostwald ripening appears in the sharp-interface limit of scalar Ginzburg-Landau equations (and its modifications) with a global conservation law [2, 3, 4, 5, 20]. Although the simple theories of Lifshitz-Slyozov and Wagner were developed more than 40 years ago, they are still very useful in phase ordering theory. For example, Sire and Majumdar [4] employed the Wagner distribution to calculate the equal-time pair correlationfunction. Lee and Rutenberg [15] used the two theories of Ostwald ripening forcalculatingtheautocorrelationexponentandpersistenceexponentforthelocally-andglobally-conservedsystems. Many works were devoted to extensions of the Lifshitz-Slyozov theory to finite volume fractions. Already Lifshitz and Slyozov [18] made an attempt to go beyond their simple model and account for coalescence. Later it became clearthat,inthe locally-conservedsystems,the dominanteffectunaccountedforby thesimple theoryisinter-domain correlations, rather than coalescence. At small area fractions, the relative role of correlations is of order ε1/2 [21], while the relative role of coalescence is of order ε. Therefore, an account of coalescence without a proper account of correlations is an excess of accuracy. The situation is quite different in globally-conserved systems, and this fact has not been recognized until now. Correlationsbetweenneighboringdomainsareexponentially smallinthiscase[3,5]. Therefore,coalescenceisexpected to give the dominant correction to the theory of Wagner [19]. We shall report numerical simulations that show a strongeffectofcoalescenceatmoderateε. Specifically,wefindthat,atlargetimes,thesizedistributionfunctionofthe minority domains has a self-similar form with the “normal” growth exponent 1/2. The scaled distribution function, however,stronglydiffers fromthe Wagnerdistribution. We attribute this difference to coalescenceanddevelopanew theory of Ostwald ripening that takes coalescence into account. Theoutlineoftherestofthepaperisthefollowing. InSectionIIwebrieflyreviewthephasefieldmodelofglobally conserved phase ordering: a scalar Ginzburg-Landau equation with a global conservation law. The sharp-interface asymptotic limit of this equation is introduced and reduced, in 2D, to a simpler model of area-preserving motion by curvature. The criteria for the validity of this reduction are obtained and presented in Appendix. The model of area-preserving motion by curvature is used to obtain dynamic scaling laws for the characteristic coarsening length and for the effective magnetic field. The results of the phase field simulations for the critical quench (ε = 0.50) and off-criticalquench(ε=0.25)arepresentedandanalyzedinSectionIII. InSectionIVanewtheoryofOstwaldripening is developed. The theory leads to a non-linear integro-differential equation for the scaled distribution function. The solution of this equation is in good agreement with the scaled size distribution function found in the phase field simulations. In Section V we summarize our results. II. PHASE FIELD MODEL AND SHARP-INTERFACE LIMIT Globally conserved phase ordering dynamics are describable by a simple phase field model [3, 4, 6]. In this model the free energy functional has the Ginzburg-Landau form: δ F[u]= ( u)2+V(u)+Hu ddr, (1) 2µ ∇ Z (cid:20) (cid:21) and the dynamics follow a simple gradient descent: δF ∂ u= µ =δ 2u+µ u u3 H(t) . (2) t − δu ∇ − − Either no-flux or periodic boundary conditions can be used. I(cid:0)n Eqs. (1) and(cid:1)(2) u(r,t) is the coarse-grained order parameter field, V(u) = (1/4)(1 u2)2 is a symmetric double-well potential, δ is the diffusion coefficient, µ is the − characteristic rate of relaxation of the field u to its stable equilibrium values, and d is the dimension of space. The effective uniform magnetic field H(t) changes in time so as to impose the global conservation law: u(r,t) =L−d u(r,t)ddr=const, (3) h i Z 3 where L is the system size and the integrationis carriedout over the entire system. Integratingboth sides of Eq. (2) over the entire system and using Eq. (3) and the boundary conditions we obtain: H(t)= u u3 =L−d u(r,t) u3(r,t) ddr. (4) h − i − Z (cid:2) (cid:3) Therefore, Eq. (2) takes the form: ∂ u=δ 2u+µ(u u3) µ u u3 , (5) t ∇ − − h − i a globally-constrainedGinzburg-Landau equation (GLE). From now on we shall concentrate on the 2D case. At late stages of the coarsening process, the system consists of domains of ”phase 1” (where u is close to 1) and − ”phase 2” (where u is close to 1) separated by domain walls. The domain walls can be treated as sharp interfaces [5,22],astheir characteristicwidthλ=(δ/µ)1/2 ismuchsmallerthanthe characteristicdomainsizel(t)whichgrows with time. At this stage H(t) is already small, H(t) 1, and slowly varying in time. The phase field in the phases ≪ 1 and 2 is almost uniform and rapidly adjusting to the value of H(t), so u 1 H(t)/2 and u 1 H(t)/2 in ≃ − − ≃ − the phases 1 and 2, respectively. Under these conditions, the so called“sharp-interface theory” holds. The normal velocity of the interface v is given by [5]: n v (s,t)=δκ(s,t)+(δµ)1/2gH(t), (6) n where s is a coordinate along the interface, κ(s,t) is the local curvature, and g = 3/√2. A positive v corresponds n − to the interfacemovingtowardsphase1,while apositiveκ correspondstoaninterfacewhichis convextowardsphase 2. The dynamics of H(t) are described by [5]: 4Λ(t) H˙(t)= δκ(s,t)+(δµ)1/2gH(t) . (7) L2 h i Here κ(s,t) is the interface curvature averagedover the whole interface: 1 κ(s,t)= κ(s,t)ds, (8) Λ(t) I and Λ(t) = ds is the total perimeter of the interface. Equations (6) and (7) provide a general sharp-interface formulation for the GLE with a global conservation law. H Let us denote by A(t) the total area of phase 2: A(t)= d2r. Zu(r,t)>0 Equations (6) and (7) can be used to calculate the rate of change of A(t): A˙(t)= v (s,t)ds=Λ(t) δκ(s,t)+(δµ)1/2gH(t) . (9) n I h i Using Eqs. (7) and (9) we obtain: H˙ =4A˙/L2 which yields the global conservation law: L2H(t) A(t) =const. (10) − 4 The secondterm in Eq. (10)correspondsto the bulk order parameteru being biasedby H(t). One canuse Eq. (10) instead of Eq. (7) in the general sharp-interface formulation of the problem. Insomeimportantcasesthisformulationcanbesimplifiedfurther[5,7,23]. Whenthetwotermsontherighthand side of Eq. (9) approximately balance each other, 1 δ 1/2 H(t) κ(s,t), (11) ≃−g µ (cid:18) (cid:19) the area of each of the two phases remains constant. In this case Eq. (6) takes the form: v (s,t)=δ κ(s,t) κ(s,t) . (12) n − h i 4 Dynamics (12) are known as area-preserving motion by curvature in 2D, and as volume-preserving motion by mean curvature in 3D [3, 5, 24]. Due to the presence of the non-local term κ this model is different from the Allen-Cahn equation [25] v =δκ, which represents the sharp-interface limit for non-conserved(model A) dynamics [1]. n A simple example where the area-preserving dynamics cannot be used is the dynamics of a single circular domain of the minority phase in a “sea” of the majority phase [5, 20]. Another example is the dynamics of a “donut”: a single domain of the minority phase with an inclusion of a majority phase domain [26]. Therefore, the first question we needto addressconcernsthe generalconditions under which the area-preservingdynamics,Eq. (12),representan accurate approximationto the more generalsharp-interfacetheory, Eqs. (6) and (7). These conditions are derivedin Appendix. Now we employ the area-preserving model and do simple dynamic scaling analysis. (In the rest of the paper we are using dimensionless variables and put δ = µ= 1.) For critical quench we have H(t) = 0 and k(s,t)= 0 because of symmetry between the two phases (we neglect finite-size effects). Therefore, the globally conserved dynamics for critical quench are identical to the non-conserved (model A) dynamics. Using the Allen-Cahn equation v (s,t) = n κ(s,t), one arrives at the well-known scaling law l(t) t1/2 [1]. ∼ Turning to the off-critical quench, we notice that, under the scaling assumption, the interface velocity can be estimated as v dl/dt. Each of the two terms on the right hand side of Eq. (12) is of order 1/l(t). Equating and n integrating, we a∼gain obtain l(t) t1/2. Therefore, global conservation does not change the dynamic scaling for any ∼ area fraction. This result was previously obtained by dynamic renormalization group arguments applied to Eq. (5) (withaGaussianwhitenoiseterm)[13],andbyparticlesimulations[4]. Fortheoff-criticaldynamicsofH(t)wehave: H(t) = <κ(s,t)>/g 1/l(t) t−1/2. | | | |∼ ∼ Though the dynamic exponent is independent of the area fraction, other characteristics can depend on it. In the following Section we report numerical simulations that address area-fraction-dependentquantities. III. NUMERICAL SIMULATIONS We performedextensive simulations by directly solving Eq.(5)with initial conditions in the formof ”white noise”. The simulations were done for two different values of the area fraction of the minority phase: ε=0.50,and ε=0.25, correspondingtoacriticalandoff-criticalquench,respectively. Inbothcasestheresultswereaveragedover10different samples. Eq. (5) was discretized and solved on a 1024 1024 domain, with mesh size ∆x = ∆y = 1 and periodic × boundary conditions. The coarsening process was followed up to a time t = 3000. An explicit Euler integration scheme was used to advance the solution in time, and the Laplace operator was discretized by second order central differences. A time step ∆t = 0.1 was required for numerical stability. The accuracy of the numerical scheme was monitored by checking the (approximate) conservation law (10) of the general sharp-interface theory. It was found thatthisconservationlawisobeyedwithanaccuracybetterthan0.02%fort>4andbetterthan0.008%fort>30in the critical quench case. In the off-critical quench the approximate conservationlaw (10) is obeyed with an accuracy better than 1% for t>30. To avoid any misunderstanding here and in the following we notice that, in all cases, the integrated order parameter [see Eq. (3)] is conserved exactly by the numerical scheme. Itisconvenienttointroduceanauxiliarydensityfieldρ(r,t)=(1/2)[u(r,t)+1]. Theminorityphaseisidentifiedas the locuswhere ρ(r,t) 1/2. Typicalsnapshotsofthe coarseningprocessareshownin Fig. 1 forthe criticalquench, ≥ andinFig. 2fortheoff-criticalquench. Forthecriticalquenchthesystemconsistsofinterpenetratingdomainsofthe twophases. Fortheoff-criticalquenchthemorphologyisthatofOstwaldripening[27]: largerdomainsoftheminority phase grow at the expense of smaller ones. As the minority phase area fraction is not very small, binary (and even triple) coalescenceevents areclearly seeninFig. 2. Overall,the coarseningmorphologiesresemble those observedfor locally-conserved system: in numerical solutions of the Cahn-Hilliard equation [28] and in particle simulations [29]. An important difference is an apparent absence of correlations between neighboring domains in Fig. 2. To analyze the coarsening dynamics, the following quantities were sampled and averagedover the 10 initial condi- tions: 1. The area of phase 2. 2. The circularly averagedequal-time pair correlation function: <ρ(r′,t)ρ(r′+r,t)> <ρ(r′,t)>2 C(r,t)= − . (13) <ρ2(r′,t)> <ρ(r′,t)>2 − 3. The characteristic coarsening length scale l(t), determined from the condition C(l,t)=1/2. 4. The effective magnetic field H(t) computed from Eq. (4). 5 5. The size distribution function of the minority phase domains (for the off-critical quench). For critical quench we found that the area of the minority phase is constant with an accuracy better than 0.03% at all times. The situation is quite different for the off-critical quench. Here there is a systematic trend in the area fraction of the minority phase. Still, with time this quantity approaches a constant value. Deviations from this constantvaluebecomelessthan3%fort>100. Thisapproximateareaconservationplaysacrucialroleinthetheory of Ostwald ripening (see Sec. IV). Figure 3 shows, on a single graph, the scaling forms of the correlation function C(x), where x = r/l(t), for the critical and off-critical quench. The l(t)-dependence is presented in Fig. 4. A comparison of the scaling forms C(x) with those obtained in particle simulations of globally conserved [6] and non-conserved [30] dynamics for critical quench is also shown. The three curves for the critical quench almost coincide. For off-critical quench, the C(x) curve is slightly different from the curvesfor the criticalquench. A similar weak dependence of the scaledcorrelation function on the area fraction of the minority phase was observed in locally conserved systems [28, 29, 31]. Figure 4 shows corrected power-law fits l(t) = l +btα which yield α = 0.50, l = 0.5, and b = 1.2 for the critical 0 0 quench, and α = 0.51, l = 1.3, and b = 0.9 for the off-critical quench. A pure t1/2 power-law line serves as a 0 reference for the expected late-time dynamic behavior. Therefore, l(t) obeys the expected t1/2 dynamical scaling law,in agreementwith the predictions ofthe dynamic renormalizationgroupanalysis[13]and area-preservingsharp- interface theory. The difference in the values of the amplitudes b againindicates a dependence on the minority phase area fraction. The time history of 1/H(t) for the off-critical quench is presented in Fig. 5. The data is fitted by a corrected power-law: 1/H(t) = a+| ctα|with a = 7.4, c = 2.3, and α = 0.51. Also shown is a 2.3t1/2 power law, serving as a reference to th|e exp|ected late-time dynamics. We conclude that H(t) t−1/2, as predicted by the sharp-interface | | ∼ theory. The significanceofthe value ofthe amplitude c willbe discussedinSectionIV in the contextofour theoryof Ostwald ripening with coalescence. One can distinguish in Fig. 5 small ”fluctuations” of 1/H(t) around a smooth | | trend. (This is incontrastto thel(t)-dependence where no fluctuationsareobserved.) To interpretthese fluctuations weuseEqs. (11)and(50)toobtain: 1/H(t) Λ(t)/N (t). Λ(t)isacontinuousfunctionoft,whereasN (t)behaves 2 2 | |∼ discontinuouslyatthetimemomentswhendomainsdisappearduetoshrinkingandmergingevents. ThusH(t)serves as a “domain counter”. For critical quench, H(t) exhibits very small irregular fluctuations around zero. The typical values of H(t) in this case are of the order of 10−5, and we interpret these fluctuations as finite-size effects. Aswehaveshown,thescaledcorrelationfunctiononlyweaklydependsontheareafraction. Amuchmoresensitive diagnostics of the off-critical quench dynamics is provided by the size distribution function of the minority phase domains. We found that, at late times, this function exhibits dynamic scaling. Figure 6 shows the scaled form Φ num of the distribution function obtained in the simulations with the GLE. The scaled variable on the horizontal axis of Fig. 6 is ξ = R/t1/2, where the effective radius of each domain is defined as R = (A /π)1/2 and A is the domain d d area. Function Φ was obtained, at each moment of time, by multiplying the values of the distribution function, num found in the simulations, by t3/2. The dynamic exponents 1/2 and 3/2 are the same as in the classical theory of Wagner [19]. Hereisamoredetailedaccountofourcalculationofthescaleddistributionfunction. Wechoseforsampling13time moments in the interval 120 < t < 2900. The domain statistics is obviously better at earlier times of this interval, and it deteriorates at later times, as many domains shrink and disappear. On the other hand, the dynamic exponent 1/2 shows up, with a good accuracy, only at relatively late times (see Figs. 4 and 5). Therefore, we had to include the relatively late times in our sampling, which led to relatively big error bars in Fig. 6. The area fraction ε = 0.25, used in our simulations, is moderately large. Therefore, one could expect significant deviations of the scaled distribution function, found numerically, from the Wagner distribution [19] corresponding to the same area fraction (that is, having the same second moment). The Wagner distribution has the following form: ξ 2√2 Φ (ξ)=Cε exp (14) W (ξ √2)4 −√2 ξ! − − for ξ <√2, and Φ (ξ)=0 for ξ √2. The normalization constant W ≥ 1 C = 16.961, π[(2e2)−1+Ei( 2)] ≃ − where Ei(...) is the exponential integral function [32]. Thetwodistributions,Φ (ξ)andΦ (ξ),areshowninthe sameFig. 6. Onecanseethatthedifferencebetween num W them is enormous (in order to show the Wagner distribution on the same graphwith Φ , we had to multiply it by num a factor of 0.5). Therefore, at moderate area fractions, the Wagner’s theory is inapplicable. 6 It is instructive to compare the zero moments M of the two distributions. The zero moment is the amplitude of 0 the scaling law for the number density of domains at large times: n(t) = M t−1. We obtained M = 4.72 10−2 for 0 0 Φ and M 1.43 10−1 for Φ . Therefore, for ε = 0.25 the Wagner distribution overestimates the n·umber of num 0 W ≃ · domains at late times by a factor of 3. An additional difference is the pronounced tail in Φ which extends much num further thanthe edgeofthe compactsupportofthe Wagnerdistribution. Coalescenceprovidesanaturalexplanation tothesetwofacts: coalescenceeventsreducethetotalnumberofdomainsandproducedomainsofprogressivelylarger size. We shall see in the next Section that an account of coalescence leads to a good quantitative agreementbetween theory and simulations. IV. THEORY OF OSTWALD RIPENING WITH COALESCENCE In this Section we present a new theory of the globally-conserved (interface-controlled) Ostwald ripening that accounts for coalescence. One of our assumptions is that each domain can be represented by an equivalent circular domain, or droplet, the area of which is equal to the area of the domain. We shall denote by f(R,t) the distribution ∞ function of the droplets with respect to their radii. f(R,t) is normalized by the condition f(R,t)dR = n(t), 0 where n(t) is the number density of the droplets. We start with a brief review of the “classical” theory that neglects R coalescence and goes back to Wagner [19]. Then we derive a kinetic equation that accounts of coalescence. We shall focus on the long-time, self-similar asymptotic solutions to that kinetic equation, find the solution by an iteration procedure and compare it with the result of the phase-field simulations. A. Ostwald ripening without coalescence: a brief review At a late stage of coarsening H(t) 1, so there is no nucleation of new domains. Then, neglecting coalescence, | | ≪ one can write a simple continuity equation in R-space for the size distribution function of domains, or droplets: ∂ f +∂ (R˙f)=0 . (15) t R Whencriterion(51)is satisfied,the dynamicsaredescribablebythe area-preservingmotionbycurvature(12)(where we put δ =1). This leads immediately to 1 1 R˙ = , (16) R (t) − R c where the time-dependent critical radius R (t) = √2/(3 H(t)) is determined, at a late stage of coarsening, by the c | | conservation of the total area of the minority phase: ∞ π R2fdR=ε=const. (17) Z0 Equations(15)-(17)representtheclassicalmodelofinterface-controlledOstwaldripening. Thismodelwasformulated by Wagner [19] by analogywith the theory of Lifshitz andSlyozov[18] developed for the locally conserved(diffusion- controlled) dynamics. Using Eqs. (15)-(17), one obtains ∞ RfdR Rc(t)= 0∞ = R(t) , (18) fdR h i R 0 where R(t) is the time-dependent average radius ofRthe droplets. h i Droplets with R > R (t) grow at the expense of droplets with R < R (t) which shrink. The late-time asymptotic c c behaviordescribedbyEqs. (15)-(17)isthefollowing[19]. Thecriticalradiusgrowswithtime(thiscorrespondstothe decreasewithtimeoftheeffectivemagneticfieldwhichplaystheroleofsupersaturation). Asaresult,adropletwhich was growing at an early time begins to shrink at a later time. Since all the quantities are position-independent, this model represents a mean-field theory. It should be noticed that the mean-field approximationis much more accurate for the globally-conserved(interface-controlled)Ostwaldripening thanfor the locally-conserved(diffusion-controlled) Ostwaldripening[5]. First,intheglobally-conservedcase,the”meanfield”H(t)istheactualfieldinthesystem. This is in contrast to the diffusion-controlled Ostwald ripening [18], where a mean field description of the supersaturation isanapproximationvalidonlywhenthe typicaldistancesbetweenthedropletsareverylargecomparedtothe typical droplet radius. The second difference concerns the role of correlations. In the locally-conserved case, correlations 7 between droplets result from the Laplacian screening effect, and their relative contribution to the size distribution function is of order ε1/2 (see, e.g. Ref. [21]). The effect of coalescence scales like ε (see below) so, at small ε, correlation effects should be much less significant. By contrast, in the interface-controlled case direct correlations between droplets are exponentially small, and significant correlations can be caused only by coalescence events. Therefore, in the interface-controlled case, it is legitimate to account for coalescence while neglecting correlations. Wagner [19] obtained a self-similar solution to Eqs. (15)-(17) (the Wagner distribution) that corresponds to a long-time asymptotics of the initial-value problem. The similarity Ansatz is 1 R t1/2 f(R,t)= Ψ , R (t)= , (19) t3/2 β t1/2 c β (cid:18) (cid:19) where β is a constant number. The scaled distribution Ψ (ξ) obeys an ordinary differential equation: β ξ 1 3 1 ′ +β Ψ (ξ)+ + Ψ (ξ)=0. (20) −2 − ξ β −2 ξ2 β (cid:18) (cid:19) (cid:18) (cid:19) The total area conservation (17) leads to normalization condition ∞ π ξ2Ψ (ξ)dξ =ε=const. (21) β Z0 Formally solving Eq. (20), one actually obtains a family of solutions parameterized by β. For √2 β 2 2/3 ≤ ≤ these solutions have compact support: they are positive on an interval 0<ξ <ξ (β), and zero elsewhere. Similar max p solutions in 3D were investigated in Refs. [5, 20, 33]. We call these solutions localized. For 0<β <√2 the solutions of Eq. (20) are extended: they have an infinite tail. These solutions can be written as const Ψ (ξ), where 0β · ξ Ψ (ξ)= 0β (ξ2 2βξ+2)2× − 2β ξ β exp arctan − . (22) − 2 β2 2 β2! − − Extended solutions fall off like ξ−3 as ξ .pAs a result, thpe integral in Eq. (21) diverges logarithmically, so → ∞ the extended solutions are non-normalizable. Still, as we shall see, they play a crucial role in the theory of Ostwald ripening with coalescence. Which of the similarity solutions is selected by the dynamics (that is, represents a long-time asymptotics of the initialvalueproblem)? Itturnsoutthatselectionis“weak”,thatis,determinedbytheinitialconditions. TheWagner distribution is selected for (normalizable) extended initial distributions. On the contrary, if the initial distribution f(R,t = 0) has compact support, one of the localized distributions is selected. The selection is determined by the asymptotics of f(R,t=0) near the upper edge of its support [5, 20, 33]. However, this weak selection rule was obtained in the framework of the classical formulation of the problem, Eqs. (15)-(17). One can expect that strong selection (independent of the initial conditions) can be obtained if one goes beyondtheclassicalformulation. Indeed,itwasshowninRef. [34](seealso[35])thatanaccountoffluctuationsleads to strong selection. Fluctuations produce a tail in the time-dependent distribution function and drive the solution towards the Wagner distribution. We shall see in the following that an account of coalescence also leads to strong selection, even in the absence of fluctuations. B. Kinetic equation with coalescence We shall now take into account the processes of binary coalescence. Coalescence events occur when two droplets contact each other. Within the framework of the GLE, the positions of the droplet centers remain fixed. Therefore, for coalescence to happen, at least one of the droplets must be expanding. Consider a droplet of radius R < 1 R < R + ∆R . The number density of such droplets is f(R ,t)∆R . Now consider another droplet of radius 1 1 1 1 R <R<R +∆R inthe vicinityofthefirstdroplet. IfR˙ +R˙ >0then, duringthetime interval∆t, thedistance 2 2 2 1 2 between the boundaries of these droplets will decrease by (R˙ +R˙ )∆t. If the distance r between the centers of the 1 2 droplets obeys the double inequality (R +R ) r (R +R )+(R˙ +R˙ )∆t, (23) 1 2 1 2 1 2 ≤ ≤ 8 (which assumes that the condition R˙ +R˙ > 0 is fulfilled), then these two droplets will collide during the time 1 2 interval ∆t. Therefore, for the two droplets to collide, the center of the second droplet should be located within a circular ring, concentric with the first droplet, with radius R +R and width (R˙ +R˙ )∆t. The area of this ring is 1 2 1 2 equal to 2π(R +R )(R˙ +R˙ )∆t . (24) 1 2 1 2 Hence, the average number of such second droplets is equal to 2M(R ,R )f(R ,t)∆R ∆t , 1 2 2 2 where M(R ,R )=π(R +R )(R˙ +R˙ )θ(R˙ +R˙ ) . 1 2 1 2 1 2 1 2 The total number of the collision events per unit area is equal to [2M(R ,R )f(R ,t)∆R ∆t]f(R ,t)∆R . (25) 1 2 2 2 1 1 Each collision leads to coalescence: disappearance of a droplet of radius R and a droplet of radius R , and creation 1 2 of a new droplet. Now we make two assumptions that will enable us to construct a closed theory. First, we assume that the area of a new droplet, formed by a binary coalescence event, is equal to the sum of the areas of the two merging droplets. Second, we assume that new droplet instantaneously becomes circular [36], so its radius is (R2+R2)1/2. The kinetic 1 2 equationfor the size distribution function includes the rates ofgains andlossesof droplets by coalescence. This leads to the following equation: ∞ ∞ 1 ∂ f +∂ (R˙f)= 2M(R ,R ) t R 1 2 −2 { × Z Z 0 0 δ(R R )+δ(R R ) δ R R2+R2 − 1 − 2 − − 1 2 × (cid:20) (cid:18) q (cid:19)(cid:21) f(R ,t)f(R ,t)dR dR , (26) 1 2 1 2 } where δ(...) is the Dirac’s delta-function and the factor 1/2 is introduced in order to avoidcounting each coalescence eventtwice. Performingintegrationwithδ(R R )andδ(R R )andtakingintoaccountthesymmetryofM(R ,R ) 1 2 1 2 − − under a transposition of its arguments: M(R ,R )=M(R ,R ), 1 2 2 1 we obtain: ∂ f +∂ (R˙f)= t R ∞ = 2f(R,t) M(R,R )f(R ,t)dR + 1 1 1 − Z 0 ∞ ∞ + M(R ,R )δ R R2+R2 1 2 − 1 2 × Z0 Z0 (cid:18) q (cid:19) f(R ,t)f(R ,t)dR dR . (27) 1 2 1 2 Integration of the right hand side of Eq. (27) over R2dR yields zero, so the new kinetic equation preserves the conservation law (17) as it should. In addition, the simple relation (18) continues to hold. Integrating the right 9 hand side of Eq. (27) over dR, and over RdR, respectively, one can show that the coalescence term reduces the number density of the droplets and the total interface length. Moreover, the new equation preserves the dynamic scaling. Indeed,ifsomef(R,t)andR (t)giveasolutiontoEqs.(27),(17)and(18),thenf′(R,t)=η3f(ηR,η2t)and c R′(t) = η−1R (η2t) give another solution to the same equations. This invariance under a stretching transformation c c implies the existence of a self-similar solution that will be considered in the next subsection. WhilederivingEq. (27),weneglectedeffectsofinteractionsofthreedroplets. Bythiswerefertocaseswherethere are three closely lying droplets. In these cases triple coalescence events may occur. In addition, an excluded area in the ring (23) appears. The effects of this excluded area, and of the triple coalescence events were not taken into account in our theory. These effects are expected to be of order ε3, while the effects of binary coalescence events are of order ε2. Therefore, Eq. (27) is expected to be valid for small area fractions ε. We shall see, however, that a very good accuracy is obtained even for the moderate value of ε = 0.25 used in our simulations, when triple coalescence events do occur (see Fig. 2). Another limitation of our theory concerns the large-R tail of f(R,t). The tail shape is affected by higher-order coalescence events unaccounted for in our theory. This limitation is not very important in practice. The main contributiontothecriticalradiusR comes,fornormalizabledistributions,bythe“body”ofthedistributionfunction, c rather than by the tail. Weconcludethissubsectionbyabriefdiscussionofadifferenttypeofcoalescence: Browniancoalescence. Following thepioneeringworkofSmoluchowski[37],BinderandStauffer[38]suggestedamean-fieldscenarioofphaseseparation in alloys in which clusters of the minority phase are regarded as Brownian particles: they perform random walk in space. Whentwoclusterscollide,they mergeinto alargersinglecluster. The correspondingkinetic equationincludes an integral term whose general structure resembles that of the integral term in Eq. (27), but with a different kernel M(R ,R ). If the cluster diffusivity is a power-law function of the cluster size, one arrives at a self-similar solution 1 2 for the size distribution function of the droplets. An important further development was the work of Siggia [39] who considered hydrodynamic interactions between randomly moving and coalescing droplets in phase separating binary fluids. Following the work of Siggia, the Brownian coalescence in binary fluids has been extensively studied theoretically and experimentally. Among important issues here is a crossover from Ostwald ripening (the Lifshitz- Slyozov-Wagnermechanism)toBrowniancoalescence[1,39,40],plethoraofhydrodynamicinteractionsintheprocess of coalescence [39, 41, 42], scaling violations [43] etc. In parallel, Brownian coalescence has been investigated in the context of coarsening of clusters of atoms or vacancies diffusing on surfaces, following particle deposition [44]. It is clear that Brownian coalescence is different in its nature from the coalescence process considered in this work. In contrast to Brownian coalescence, droplets in our system do not move: they coalesce only because they grow. C. Self-similar solution with coalescence Equations (27) and (17) admit the same similarity Ansatz as Eqs. (15) and (17): 1 R f(R,t)= Φ (28) t3/2 t1/2 (cid:18) (cid:19) and R (t)=β−1t1/2, (29) c whereβ isagainanunknownyetconstantnumber. ThescaleddistributionfunctionΦ(ξ)obeysthe followingintegro- differential equation: ξ 1 3 1 ′ +β Φ(ξ)+ + Φ(ξ)= −2 − ξ −2 ξ2 (cid:18) (cid:19) (cid:18) (cid:19) ∞ 2Φ(ξ) w(ξ,ξ )Φ(ξ )dξ + 1 1 1 − Z 0 ∞ ∞ w(ξ ,ξ )δ ξ ξ2+ξ2 Φ(ξ )Φ(ξ )dξ dξ (30) 1 2 − 1 2 1 2 1 2 Z0 Z0 (cid:18) q (cid:19) 10 subject to normalization condition ∞ π ξ2Φ(ξ)dξ =ε. (31) Z0 In Eq. (30) we denoted 1 1 w(ξ ,ξ )=π(ξ +ξ ) 2β 1 2 1 2 − ξ − ξ × (cid:18) 1 2(cid:19) 1 1 θ 2β , (32) − ξ − ξ (cid:18) 1 2(cid:19) where θ(...) is the theta-function. It is convenient to rewrite Eq. (30) in a symbolic form: Φ= [Φ], (33) β β L N where ξ 1 3 1 ′ Φ(ξ)= +β Φ(ξ)+ + Φ(ξ) , (34) Lβ −2 − ξ −2 ξ2 (cid:18) (cid:19) (cid:18) (cid:19) and ∞ [Φ](ξ)= 2Φ(ξ) w(ξ,ξ )Φ(ξ )dξ + β 1 1 1 N − Z 0 ∞ ∞ w(ξ ,ξ )δ ξ ξ2+ξ2 Φ(ξ )Φ(ξ )dξ dξ . (35) 1 2 − 1 2 1 2 1 2 Z0 Z0 (cid:18) q (cid:19) One important property of Eq. (30) can be noticed immediately: the coalescence term vanishes identically at 0 ξ < 1/(2β). As a result, the scaled distribution function at 0 ξ < 1/(2β) should coincide (up to a ξ- ≤ ≤ independent multiplier) with one of the solutions of the classical Wagner’s problem. A simple argument shows that parameter β, parameterizing this solution, should be less than √2. Indeed, inverting the linear operator , we β L rewrite Eq. (33) as an integral (rather than integro-differential) equation: ∞ ′ ′ [Φ(ξ )] dξ β Φ(ξ)=Ψ (ξ) N +C , (36) β  ξ′ β+ 1 Ψ (ξ′) 1 Zξ 2 − ξ′ β    (cid:16) (cid:17)  where functions Ψ (ξ) were introduced in subsection A and C is a constant. Unless β < √2, the integral over dξ′ β 1 diverges. Therefore, Ψ (ξ) should be one of the extended solutions Ψ (ξ), given by Eq. (22). In addition, since the β 0β secondterminthesquarebracketsofEq. (36)wouldleadtodivergenceoftheintegralappearinginthenormalization condition (31), we must choose C =0. Hence, Eq. (36) reads: 1 ∞ ′ ′ [Φ(ξ )] dξ β Φ(ξ)=Ψ (ξ) N . (37) 0β Zξ ξ2′ −β+ ξ1′ Ψ0β(ξ′) (cid:16) (cid:17) Integral equation (37) and normalization condition (31) make a complete set. For a given ε, the scaled distribution function Φ = Φ (ξ) and parameter β = β(ε) are uniquely determined. Therefore, an account of coalescence does β provide strong selection to the problem of Ostwald ripening.