Biscale Chaos in Propagating Fronts Anatoly Malevanets, Agust´ı Careta∗ and Raymond Kapral Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, Canada M5S 1A1 6 9 9 1 Abstract n a J 0 Thepropagatingchemicalfrontsfoundincubicautocatalytic reaction- 1 1 diffusion processes are studied. Simulations of the reaction-diffusion v 4 equation near to and far from the onset of the front instability are 0 0 performed and the structure and dynamics of chemical fronts are stud- 1 0 6 ied. Qualitatively different front dynamics are observed in these two 9 / n regimes. Close to onset the front dynamics can be characterized by a y d single length scale and described by the Kuramoto-Sivashinsky equa- - o a tion. Far from onset the front dynamics exhibits two characteristic h c lengths and cannot be modeled by this amplitude equation. An am- : v i X plitude equation is proposed for this biscale chaos. The reduction of r a the cubic autocatalysis reaction-diffusion equation to the Kuramoto- Sivashinsky equation is explicitly carried out. The critical diffusion ratio δ , where the planar front loses its stability to transverse pertur- c bations, is determined and found to be δ = 2.300. c Typeset using REVTEX 1 I. INTRODUCTION Propagatingfrontsseparatingregionswithdifferentcharacteristicsoccurinmany physical contexts. Often such fronts have a complicated structure of their own as in the fractal forms that arise in some diffusion-limited aggregration processes or viscous fingering [1]. In this article we consider chemical fronts separating regions of distinctly different chemical concentrations and study the nature of the spatio- temporal dynamics that these fronts exhibit. In these systems, under appropriate conditions, aplanarfrontmaynotbestableandcomplex, evenchaotic, frontdynam- ics can arise. For example, this is the situation often encountered in propagating flame fronts which have been studied extensively [2]. The complex dynamics of fronts can also underlie and determine the character of spatio-temporal chaos and dynamics in a variety of other chemical contexts [3]. We study a specific system here, the fronts in the cubic autocatalysis reaction A + 2B 3B, but several parts of our analysis are general and should apply in → other situations. These cubic autocatalysis fronts have many features in common with propagating flames [4]. An experimental realization of such chemical front in- stabilities occurs in the iodate-arsenous acid reaction carried out in a gelled medium to suppress fluid flow [5]. The experiments show patterned fronts, much like those in quadratic and cubic mixed-order models. For the cubic autocatalysis model, when the diffusion coefficient of the fuel A is sufficiently larger than that of the autocatalyst B, the planar front is unstable to transverse perturbations. We have investigated the nature of the resulting front dynamics as a function of the diffusion ratio D /D in large systems. Several A B results have emerged from our study of this system. For small enough diffusion ratios within the unstable regime the front exhibits complicated dynamics with statistically stationary properties and is characterized by a single length scale. Far 2 beyond the instability point a new chaotic regime is encountered, termed biscale chaos, in which there are two characteristic length scales. We have carried out a detailed analysis of the front dynamics in terms of am- plitude equations. Just beyond the instability threshold the front dynamics is de- scribed in terms of the Kuramoto-Sivashinsky equation [6,7]. We document this relation through a statistical characterization of the interface. The regime far from the instability point where biscale chaos is observed cannot be described in terms of the Kuramoto-Sivashinsky equation and we construct an amplitude equation in which the nonlinear mode coupling is generalized. This amplitude equation is able to capture the principal qualitative features of the biscale chaos. This part of the analysis may generalize to other situations. We also give a detailed reduction of the cubic autocatalysis reaction-diffusion equation to the Kuramoto-Sivashinsky equation. The coefficients that enter in this equation are estimated analytically using solutions for the right and left eigenvectors oftheeigenvalue problem thatenters theanalysis, andanumerical scheme isdevised that allows the computation of these coefficients for any values of the diffusion coefficients. As a by-product of this analysis we may easily determine the critical diffusion ratio where the planar front becomes unstable. II. CUBIC AUTOCATALYSIS FRONTS Consider the autocatalytic reaction A+ 2B 3B described by the reaction- −→ diffusion equation, ∂α = αβ2 +D ∆α , A ∂t − ∂β = αβ2 +D ∆β , (1) B ∂t 3 where α(t,r) and β(t,r) are the (scaled) concentrations of the A and B species, respectively, with diffusion coefficients D and D . In one dimension, with suitably A B defined initial conditions [8], this system supports isothermal chemical fronts [9]. In two space dimensions the planar front is unstable to transverse perturbations when D is sufficiently larger than D . The origin of such instabilities and the A B nature of the resulting front dynamics have been the subject of earlier studies [4]. Horvath et al. [10] performed simulations of the cubic autocatalysis model (1) and investigated the bifurcation structure as a function of the system length. For small system sizes, when the front possesses a single minimum, the onset of chaotic be- haviour was observed to occur through a period doubling cascade in the temporal dynamics of the minimum. For large system sizes the resulting front dynamics is best analysed in statistical terms. We now briefly describe and quantitatively char- acterize the nature of the front dynamics for large system lengths as a function of the diffusion coefficient ratio, δ = D /D . A B We consider a two-dimensional system which is infinite in the x direction and has length L along y, where periodic boundary conditions are applied. The initial conditions are taken to be α(0,x,y) = 1, β(0,x,y) = 0 x w/2, y and { | ≤ − ∀ } α(0,x,y) = 0, β(0,x,y) = 1 x w/2, y . The region x,y w/2 < x < { | ≥ ∀ } { | − w/2,0 y L was divided into segments oflength ℓ alongy and each segment was ≤ ≤ } assigned the values (α(0,x,y),β(0,x,y))= (1,0) or (0,1) with probability p = 1/2. It is convenient to represent the front dynamics in terms of the isoconcentration profiles, h (t,y) and h (t,y) defined by α(t,h (t,y),y) = α and β(t,h (t,y),y) = A B A r B β , where α and β are reference concentrations. Henceforth we focus on h (t,y) r r r B measured relative to its average value at time t and denote this quantity by h(t,y). We note that for these initial conditions this average value moves with the minimum front speed, which is selected from the continuum of allowed front speed values. 4 Figure 1(a) shows a space-time plot of the minima of h(t,y) for δ = 5. In the integration [11] of (1) we have scaled the length variables so that D = 1 and A D = δ−1. After an initial transient period during which the interface develops, a B moving front with statistically stationary properties is formed. In this stationary regime the spatio-temporal behaviour can be described in terms of the dynamics of theminimawhichplaytheroleof“particles”inthesystem. Theparticlescollideand coalesce and new particles are born so that the average density of particles per unit length remains constant. The front dynamics may be characterized quantitatively by the space-time correlation function L C(t,y) = L−1 dy′[h(t,y′) h(t,y +y′)]2 , (2) h − i Z0 where the angle brackets signify an average over realizations of the front evolution starting from the random initial conditions given above. We also consider the time average of C(t,y) in the stationary regime, t0+T C¯(y) = T−1 dt′ C(t′,y) , (3) Zt0 where t is a time longer than the transient period. The correlation function (2) 0 is sketched in Fig. 2 for several values of the time t and shows the development of correlations during the transient period and the approach to the stationary regime. The function C(t,y) exhibits short range correlations with a tendency to reach a nearly parabolic form, seen in the C¯(y) curve, which is characteristic of diffusive transverse front dynamics and will be discussed further below. The persistent short range correlations are due to the existence of a characteristic length ℓ∗ arising from the front instability. This type of behaviour persists up to about δ = 6 when a new type of chaotic frontdynamics, distinguishedbytheappearanceofasecondlengthscale, isobserved. 5 The second length scale is easily discerned in both the space-time plot of the minima of h(t,y) shown in Fig. 1(b) for δ = 8 and the power spectrum, t0+T E(k) = T−1 dt h(t,k) 2 . (4) h| | i Zt0 The power spectrum is presented in Fig. 3 for both δ = 5 and δ = 8. For δ = 5 one sees a single peak corresponding to the wavenumber of the most unstable mode. However, for δ = 8 the minimum in the power spectrum has filled in, indicating the growth of modes with smaller wavenumbers and thus the appearance of structure on longer length scales. We now turn to an analysis of these results in terms of amplitude equations and provide a foundation for the phenomenological picture of the front dynamics. III. ANALYSIS OF FRONT DYNAMICS A description of the origin of the front dynamics is most conveniently given in terms of amplitude equations for the front profile derived from the reaction-diffusion equation (1). We analyse the dynamics close to the instability point in terms of the Kuramoto-Sivashinsky equation and show how this equation must be generalized in order to explain the regime far beyond the instability where a second characteristic length appears. A. Dynamics of small perturbations We present an adaptation of Kuramoto’s derivation [12] for the case of small amplitude perturbations which yields the amplitude equation that will be used in the subsequent analyses. We shall confine our attention to the question of how the presence of a second dimension affects the dynamics of the general reaction-diffusion system, 6 ∂z = F(z)+D∆z . (5) ∂t Here z is a vector of concentration fields and F(z) is a vector-valued function de- scribing chemical reactions. We assume that in one dimension (5) possesses a stable solution with a propagating front profile z(t,x,y) = z (x ct), where c is the veloc- 0 − ity of the front. Following the treatment by Kuramoto [12] we seek the dynamics of a perturbed solution in the form, z(t,x,y) = z (ξ +φ (t,y))+ φ (t,y)u (ξ) , (6) 0 0 i i i>0 X where ξ = x ct. As is customary, we use an abstract notation ξ u = u (ξ). The i i − h | i vectors u are the solutions of the eigenvalue problem, i | i ∂2 ∂ "DFˆ(z0)+D∂ξˆ2 +c∂ξˆ#|uii ≡ Lˆ|uii = λi|uii , (7) where F is the Jacobian of the vector-function F and ˆ is the operator in square D L brackets. Here a hat signifies an abstract operator while the corresponding quantity without a hat is its ξ representation. We expand F(z) in a Taylor series near z 0 and neglect terms of second order and higher. Since the main features of the front dynamics are embodied in F, the truncation of the Taylor series does not discard D any important information. Due to the translational symmetry of equation (5) we have a straightforward solution of (7), namely u = ∂ z /∂ξˆ corresponding to a 0 0 | i | i zero eigenvalue. Taking this solution into account we write ∂φ ∂2 ∂ i u = c u +Fˆ(z )+D z +( φ )2D u + ∂t | ii | 0i 0 ∂ξˆ2 0 ∇ 0 ∂ξˆ| 0i (8) ∂2 ∂ +∆φ D u + φ F u +D u +c u . i | ii i"D | ii ∂ξˆ2| ii ∂ξˆ| ii# i>0 X ∂2 Using the identity c u = F(z ) + D z and definition (7) we arrive at the − | 0i 0 ∂ξˆ2| 0i expression 7 ∂φ ∂ i u = λ φ u +∆φ D u +( φ )2D u . (9) ∂t | ii i i| ii i | ii ∇ 0 ∂ξˆ| 0i Here and below we use the Einstein summation convention. Alternatively, after multiplication by u , (9) may be represented by a set of coupled equations i h | ∂φ ∂ i = λ φ + u D u ∆φ + u D u ( φ )2 . (10) ∂t i i h i| | ji j h i| ∂ξˆ| 0i ∇ 0 This set of equations can be formally solved for modes φ (i = 0) by treating i 6 φ as an independent function and applying Duhamel’s principle to the resulting 0 system of inhomogeneous linear equations. We find t φ(t) = eWtφ(0)+ dseW(t−s) a∆φ (s)+b( φ (s))2 . (11) 0 0 ∇ Z0 h i c c In (11) we use the following notation: the matrix operator W has elements W = δ λ + u D u ∆, the vectors a and b have elements ca = u D u , ij ij i i j i i 0 h | | i h | | i ∂ bc= u D u , respectively, and φ is a vector with elements φ (i > 0). i h i| ∂ξˆ| 0i i Assuming that exp(Wt) decays rapidly compared to φ we may perform the 0 integration in (11) and scubstitute the result into the equation for φ to obtain 0 ∂φ 0 = a c W−1 a ∆ ∆φ + b c W−1 b ∆ ( φ )2 , (12) 0 i ij j 0 0 i ij j 0 ∂t  − { }   − { }  ∇ i,j>0 i,j>0 X X  c   c  where c = u D u . This is the generalized amplitude equation which will form i 0 i h | | i the basis of the analysis of the biscale front chaos given below. If the dynamics is described by small k modes one can take W−1 to bea diagonal matrix δ /λ and by omitting second order terms in the coecfficient of ( φ)2 one ij i ∇ obtains [13] ∂φ ∂ u D u u D u 0 = u D u ∆φ + u D u ( φ )2 h 0| | iih i| | 0i∆2φ . (13) ∂t h 0| | 0i 0 h 0| ∂ξˆ| 0i ∇ 0 − λ 0 i>0 i X The coefficient of the ( φ )2 term is obtained from the following identity: 0 ∇ 8 x 0 u ˆξ ξ u u ξ ξ ˆu dξ 0 0 0 0 ≡ h |L| ih | i−h | ih |L| i −Z∞ h i ∂ ∂ = 2 u ξ ξ D u + c D u ξ ξ u . (14) " h 0| ih | ∂ξˆ| 0i − ∂ξ!h 0| ih | 0i#(cid:12) (cid:12)x (cid:12) (cid:12) The above formula vanishes for all x so that after integrating(cid:12)over the entire do- ∂ c main we obtain u D u = and observe that the expression has a universal h 0| ∂ξˆ| 0i −2 character. This gives the Kuramoto-Sivashinsky equation, ∂φ c 0 = ν∆φ ( φ )2 κ∆2φ , (15) 0 0 0 ∂t − 2 ∇ − with ν = u D u , (16) 0 0 h | | i and u D u u D u 0 i i 0 κ = h | | ih | | i . (17) λ i>0 i X A stability analysis of (10) shows that the planar front is unstable with respect to long-scale, small-amplitude perturbations when ν < 0. Inthenextsubsection weanalysethecubicautocatalysisfrontdynamicsinterms ofboththegeneralizedamplitudeequation(12)aswellastheKuramoto-Sivashinsky equation (15). B. Amplitude equation description of front dynamics 1. Dynamics near instability The Kuramoto-Sivashinsky equation provides insight into the nature of the cubic autocatalysis fronts close to the instability point, and we briefly comment on this connection. Simulations of the Kuramoto-Sivashinsky equation produce height cor- relation functions and space-time plots (cf. Fig. 4) of the front dynamics similar to 9 those of the cubic autocatalysis model. In terms of the Kuramoto-Sivahinsky equa- tion, one sees that the front dynamics for δ > δ is determined by the instability at c long wavelengths, due to the fact that ν < 0 in the ∆φ term, and the dissipation 0 at short wavelengths controlled by the ∆2φ term with positive κ. The dispersion 0 relation of the Kuramoto-Sivashinsky equation is ω(k) = νk2 κk4 with maximum − − at k = ( ν/2κ)1/2. The long wavelength structure drives the dynamics of the ex- m − trema (“particles”) whose collision dynamics provides the dissipation of the energy. The characteristic distance between extrema can be related to the wavevector by ℓ∗ = 2π/k . Using a zeroth order approximation to the ν and κ coefficients given in m the next section, for δ = 5 we obtain ℓ∗ = 30 which is is comparable to the value of ℓ∗ = 39 found in the numerical simulations. Using the methods given in Sec. IV it is possible to obtain ν and κ accurately and improve this estimate. For the present purposes thisroughestimate plustheappearanceofthespace-timeplot [14]inFig.4 suffice to confirm the general character of the Kuramoto-Sivashinsky mechanism for the chaotic cubic autocatalysis front dyanmics for large system lengths. It has been argued [15] that the long-scale dynamics of the one-dimesional Kuramoto-Sivashinsky equation can be described by the Kardar-Parisi-Zhang equa- tion [16], ∂φ (t,y) ∂2φ (t,y) c 0 = 0 + ( φ (t,y))2 +ξ(t,y) , (18) ∂t D ∂y2 2 ∇ 0 where c is the front speed, is a diffusion coefficient and ξ(t) is a Gaussian white D noise source with correlation function ξ(t,y)ξ(t′,y′) = 2Γδ(t t′)δ(y y′). Al- h i − − thoughthereisdebateaboutsomeaspectsofthereductioninhigher dimensions[17], numerical simulations [18] of the scaling properties of the Kuramoto-Sivashinsky equation have confirmed this reduction in one dimension. The diffusive character of the cubic autocatalysis front dynamics discussed in Sec. II is consistent with such a description and provides additional confirmation of the applicability of the 10