Effects of intrinsic noise on a cubic autocatalytic reaction diffusion system
EEffects of intrinsic noise on a cubic autocatalytic reaction diffusion system
Fred Cooper,
1, 2, ∗ Gourab Ghoshal, † and Juan P´erez-Mercader
1, 2, ‡ Department of Earth and Planetary Sciences, Harvard University, Cambridge, MA 02138 The Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA (Dated: September 9, 2018)Starting from our recent chemical master equation derivation of the model of an autocatalyticreaction-diffusion chemical system with reactions U + 2 V λ → V ; and V µ → P , U ν → Q , we determinethe effects of intrinsic noise on the momentum-space behavior of its kinetic parameters and chemicalconcentrations. We demonstrate that the intrinsic noise induces n → n molecular interaction pro-cesses with n ≥
4, where n is the number of molecules participating of type U or V . The momentumdependences of the reaction rates are driven by the fact that the autocatalytic reaction (inelasticscattering) is renormalized through the existence of an arbitrary number of intermediate elasticscatterings, which can also be interpreted as the creation and subsequent decay of a three bodycomposite state σ = φ u φ v , where φ i corresponds to the fields representing the densities of U and V .Finally, we discuss the difference between representing σ as a composite or an elementary particle(molecule) with its own kinetic parameters. In one dimension we find that while they show markedlydifferent behavior in the short spatio-temporal scale, high momentum (UV) limit, they are formallyequivalent in the large spatio-temporal scale, low momentum (IR) regime. On the other hand intwo dimensions and greater, due to the effects of fluctuations, there is no way to experimentallydistinguish between a fundamental and composite σ . Thus in this regime σ behave as an entity untoitself suggesting that it can be effectively treated as an independent chemical species. PACS numbers: 82.40.Ck, 11.10.–z, 05.45.–a, 05.65.+b
I. INTRODUCTION
Reaction diffusion (RD) systems are a versatile class ofmodels capable of encoding a variety of phenomena ob-served in Nature in areas encompassing physics, biology,ecology, chemistry and many other fields [1–4]. Their ap-plication to chemical systems are of particular interest,as they include biologically relevant phenomena such aspattern formation and self-replication, and therefore canbe used as proxies for high-level biological systems [5, 6].While RD systems have been mostly studied from a de-terministic standpoint, any faithful application to biolog-ical systems must take into account the effects of noise,since an important facet of such systems is its exchangeof matter and energy with the environment; a processwhich clearly brings in some amount of stochasticity.To reflect this, there have been recent efforts to studystochastic chemical reaction diffusion systems [7–9]. Inparticular, when such systems are coupled with external noise it is known that there are renormalization effectsdue to the fluctuations represented by the noise (by ex-ternal we mean fluctuations that are not inherent to thechemistry itself). These affect for example, the strengthof the chemical interactions that, in turn, induce newinteractions not originally present in the “macro-level”chemistry [8]. However, independently of the above,there is also some form of intrinsic noise in the chemical ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] system whose effects are less understood. Qualitatively,one can interpret this intrinsic noise as a manifestationof the underlying mechanisms that lead to the observedbehavior of reaction diffusion systems at the level of theirchemical kinetics. Because of this, it is important to un-derstand the precise nature and effects of intrinsic noise,as it might give intuition and provide hints for under-standing the internal structure of the system [10].In light of this, in this paper we seek to determinewhether the inherent stochasticity in the nature of thechemical reactions themselves leads to effects similar tothose induced by the external or extrinsic noise. Ofcourse this particular stochasticity is restricted by certainassumptions when we attempt to model these reactions interms of kinetics. The most basic is that the moleculesare random walkers in a d -dimensional space and thatcollisions between them occur as a function of the proba-bility of encounters between these random walkers. Mostcollisions are elastic and do not result in a chemical re-action, whereas comparatively few are inelastic and leadto the actual chemistry that we are interested in. Therelative scarcity of the latter with respect to the formerimplies that the chemically interesting inelastic collisionsare effectively statistically independent and therefore the chemical reactions (occurring at large scales) are Marko-vian in nature. (Note that the assumption of the chemi-cal reactions as a Markovian process is valid only up to aresolution limit which corresponds to the mean free path of the molecules involved in the reaction.)In a previous paper [9] we considered as a test casea generic spatially extended set of macroscopic chemical a r X i v : . [ c ond - m a t . s t a t - m ec h ] J un reactions, U + 2 V λ → V, V µ → P, U ν → Q, f → U. (1)There is a cubic autocatalytic step for V at rate λ , anddecay reactions at rates µ, ν that transform V and U into inert products P and Q . Finally, U is fed into thesystem at a rate f and both U and V are allowed todiffuse with diffusion constants D u and D v respectively.We determined the form of the intrinsic noise associatedwith this system of reactions through a procedure whichtook us from the Master equation describing its chemistry(and that captures its Markovian nature) to an effectivenon–equilibrium field theory action (Cf. Appendix. A).This enabled us to derive a set of Langevin equationsthat incorporated the effects of the intrinsic noise. Thestructure of the noise was described through unique cor-relation functions. These required the existence of a col-lective mode σ = λ φ u φ v where the φ i are the fieldsencoding the chemical concentrations of the species.In this sequel paper we focus on the effects of this noise.Specifically we seek to determine whether the physi-cal parameters of the model inherit a scale-dependenceas a function of the noise, and if so, whether this in-duces new interactions (relevant and irrelevant) apartfrom those initially present in the macroscopic chemistryspecified by Eq. (1). We answer this question positivelyand consequently investigate the spatiotemporal scalesat which these induced interactions manifest themselvesalong with the momentum space behavior of the parame-ters and chemical concentrations in the limiting regimes.To denote the limits we employ the field theory “jargon”whereby the large spatiotemporal scales correspondingto low wave-number and frequency are referred to collec-tively as the infrared (IR) regime, whereas the small spa-tiotemporal scales corresponding to high wave-numberand frequency is the ultraviolet (UV) regime.Since our goal is to determine the scaling behavior be-tween the two limits, a natural way to proceed is to usea renormalization group approach [11, 12] on the many-body description of (1) that we studied in [9]. Inter-estingly, we find that the irreversibility of the reactionsleads to a time-directionality associated with the Feyn-man diagrams describing the interactions. This in turn severely restricts the possible topologies of the graphsand allows us to carry out a systematic exact calculationto determine the effect of the noise.We find that the strength of the coupling λ that reg-ulates the autocatalytic part of the chemical reaction isrenormalized due to the Markovian nature of the pro-cess. Specifically the chemically relevant inelastic colli-sions ( U +2 V → V ) proceed through an arbitrary inter-mediate number of elastic collisions ( U + 2 V → U + 2 V )modifying the coupling strength as we change scales.This scale dependence manifests itself only up to a criticaldimension d c = 1, above which (in the absence of cutoffs)the coupling constant is formally zero in the IR. Thus be-yond d c = 1 the system must be considered an effectivefield theory , whose parameters must be determined by low momentum experiments. The number of parame-ters of this effective theory are most easily described interms of recasting it via composite fields for the concen-trations φ u φ v , φ u φ v and φ v . In particular the reaction U + 2 V → V is then interpreted as proceeding throughthe formation of an intermediate state σ = λ φ u φ v thatwas also essential in determining the Langevin equationsfor (1) that were derived in [9].When probing the system at short spatiotemporal(UV) scales it turns out that there are two separatemanifestations of σ (i) a composite bound state or (ii)an elementary particle (molecule) with its own bare ki-netic terms each leading to different equations of motion.However, starting from d = 2 both versions of σ leads tothe same equations of motion. The physical implicationof this is that one cannot experimentally resolve σ into itsconstituents and that it behaves as an elementary parti-cle. In order to see its composite nature one would needto probe at spatio-temporal scales shorter than that as-sociated with the mean free path of the chemicals. How-ever in this limit the Markovian assumption is violatedand new physics is required to describe the chemistry.The fluctuations also lead to relevant new noise-induced interactions of the form n → n for n ≥ n is the number of molecules entering or leaving the re-action zone. Although all these higher order processesare naively divergent, to regulate them, one might thinkthat an infinite number of “counter terms” need to beadded. In two dimensions, however, it turns out that byagain introducing composite concentration fields for thedi-molecules φ u φ v and φ v , the divergences in these in-duced interactions can be regulated by introducing justtwo new effective field theory parameters: the decay rates(masses) M , for φ u φ v and φ v . Finally in three dimen-sions no new parameters are needed, however the equa-tions of motion are modified: in order to describe the in-frared physical chemistry at small momentum one needsto introduce higher order kinetic terms into the actionfor σ . II. MOMENTUM-DEPENDENT REACTIONRATE
In order to uncover the momentum scale behavior ofthe constituents of the model, we will resort to the manybody description of Eq. (1). In this description the reac-tion U +2 V → V is interpreted in terms of a three bodyinelastic collision where two particles of V and one of U are destroyed at an interaction vertex to create three of V . When one writes down the master equation describ-ing this reaction (see Eq. (A1)), one finds that in orderto conserve probabilities, it is also necessary to includethe elastic collision U + 2 V → U + 2 V . A graphicalrepresentation of this is shown in in Fig. 1, where the in-elastic scattering is proportional to − λ , while the elasticscattering is proportional to λ . Following this, the Doi-Peliti operator technique [13, 14] is used to write down u v v ˜ v ˜ v ˜ u u v v ˜ v ˜ v ˜ v (a) (b) FIG. 1: Many body interpretation of the reaction U + 2 V → V . (a) Inelastic scattering where one molecule of U and twomolecules of V are destroyed at an interaction vertex (shaded)to create three molecules of V . (b) The elastic scatteringwhere U + 2 V → U + 2 V which must be present in orderfor particle number and probability conservation. Time flowsfrom right to left. an equivalent non-equilibrium field theory action [15–17]thus, S = (cid:90) dx (cid:90) τ dt (cid:2) φ (cid:63)v ∂ t φ v + D v ∇ φ (cid:63)v ∇ φ v + φ (cid:63)u ∂ t φ u + D u ∇ φ (cid:63)u ∇ φ u + µ ( φ (cid:63)v − φ v + ν ( φ (cid:63)u − φ u − f ( φ (cid:63)u − − λ ( φ (cid:63)v − φ (cid:63)u ) φ (cid:63) v φ v φ u (cid:3) , (2)where the φ i ’s are fields representing the concentrationsof the chemicals. (An outline of this method is shownin Appendix A, see [9] for more extensive details of thederivation.)Since the action must be dimensionless, if we introducea momentum scale [ q ] = κ along with the diffusive tem-poral scaling [ t ] = κ − , it is straightforward to see thatwe have associated scaling dimensions:[ φ v,u ] = κ d , [ φ (cid:63)v,u ] = κ , [ λ ] = κ (cid:15)/d c , (3)where the starred fields are chosen to be dimensionlessby convention and (cid:15) = d c − d , where d c is the criticaldimension (i.e. the dimension below which λ has non-vanishing finite value in the high-momentum UV limit).Generalizing the cubic interaction φ v φ u in (2) to the form φ kv φ lu , the general form of the critical dimension is seento be d c ( k, l ) = 2 / ( k + l −
1) which in this case (where k = 2 , l = 1) implies that d c = 1. Consequently the di-mensionless version of the coupling constant λ is there-fore the combination λ κ − − d ) .It is worth noting that the action shown in (2) cor-responds to a “symmetric” phase in the sense that itpossesses a U (1) symmetry due to particle number con-servation (essentially equivalent to a form of the classicLavoisier’s principle). The graphs that contribute to any one particle irreducible vertex process or any intermedi-ate state, must be constructed from the (directed) basicinteractions shown in Fig. 1. The U (1) symmetry of thebasic 3 → U + 2 V → V (as well as the correspondingelastic reaction U + 2 V → V ), along with the unidirec-tionality of the reactions (being non-reversible) severelyrestrict the set of graphs that can be constructed. In fact looking at the topological structure of the bare vertices inFig. 1, it becomes clear that there is no combination thatcan generate any diagram contributing to corrections tothe propagator ( φ i ↔ φ i ), since this requires (graphi-cally) one incoming and one outgoing line. This impliesin turn that there is no momentum dependence of thediffusion constants D u , D v or the decay terms µ, ν .On the other hand, as shown in Fig. 2, it is clearly pos-sible to rewrite the bare three body interaction in termsof a diagrammatic expansion corresponding to a pertur-bation series in λ . Each component in the expansion cor-responds to an elastic re-scattering represented througha loop I ( p, t ) (the subscript 2 reflects the fact that theelastic re-scattering consists of two loops). Mathemati-cally this is expressed through the vertex function Γ (3 , (where the 3 refers to the number of incoming and out-going lines in Fig. 1) whose formal expression isΓ (3 , ( p, t ) = λ δ ( t − t ) − λ I ( p, t − t ) + λ × (cid:82) t t dt (cid:48) I ( p, t − t (cid:48) ) I ( p, t (cid:48) − t ) − . . . , (4)where the subscript i refers to inelastic. (The vertex func-tion for elastic scattering Γ e is just -Γ i .)Assuming there is an external momentum p flowinginto the graph from the three external legs, the explicitexpression for the loop I ( p, t ) is, I ( p, t ) = 3! (cid:90) (cid:89) i =1 (cid:18) d d p i (2 π ) d (cid:19) e − [ D v ( p + p )+ D u p +2 µ + ν ] t × (2 π ) d δ ( p − p − p − p ) . (5)Through appropriate linear transformations and (legal)shifts of the momentum variables, it is easily shown that I ( p, t ) = B ( d ) t − d exp (cid:104) − (cid:16) α + ˜ D p (cid:17) t (cid:105) , (6)where to avoid clutter we use the abbreviations B ( d ) = πD ) d , D = (cid:112) D v (2 D u + D v ) , ˜ D = D u D v D and α = 2 µ + ν . Taking the Laplace transform of (4) rendersit into a geometric sum via the convolution theorem andyields, Γ (3 , (˜ s ) = λ λ I (˜ s ) ,I (˜ s ) = B ( d )Γ(1 − d ) (˜ s + α ) − (1 − d ) , (7)where Γ( x ) is the Euler-Gamma function and ˜ s = s +˜ D p . From (7) it is clear that λ is the reaction rate inthe absence of fluctuations (since in that case I (˜ s ) = 0),whereas the full expression Γ (3 , reflects its function onthe energy and momentum scale at which it is measured.It is therefore interpreted as the dimension-full running reaction rate.The fluctuation term I (˜ s ) has a pole in one dimensionand in every positive integer dimension. For physicalreasons we restrict ourselves to the range 1 ≤ d ≤ u v v ˜ v ˜ v ˜ v u v v ˜ v ˜ v ˜ v u v v ˜ v ˜ v ˜ v ++ + . . . = u v v ˜ v ˜ v ˜ v FIG. 2: The diagrammatic expansion for the vertex function Γ (3 , for inelastic scattering. The same set of primitive divergencescontribute to that for elastic scattering, which is the same except for the substitution of a single outgoing leg from φ v to φ u .Consequently both processes are renormalized through a single renormalized coupling constant. therefore we expose the relevant divergences through theidentity Γ( n + 1) = n Γ( n ), to rewriteΓ(1 − d ) = Γ(4 − d )(1 − d )(2 − d )(3 − d ) . (8)In order to carry out the correspondingrenormalization—as in standard quantum field theory—it is convenient to introduce a dimensionless counterpartto λ . As mentioned before the dimensionless bare reaction rate in the presence of a momentum scale κ is λ κ − − d ) . However, in order to simplify the algebrawe will instead find it convenient to work with a slightlymodified form thus g ( κ ) = λ κ − − d ) B ( d )Γ(4 − d ) (9)in terms of which the dimensionless running reaction rateis modified to g R (˜ s ) = Γ (3 , ( p, s ) B ( d )Γ(4 − d ) κ − − d ) = g ( κ )1 + g ( κ )(1 − d )(2 − d )(3 − d ) (cid:0) ˜ s + α κ (cid:1) − (1 − d ) . (10)By itself, this does not seem particularly useful, since weneed to translate this into a physically measurable quan-tity. To do so we begin by defining a renormalized cou-pling constant g R ( κ ) at the convenient renormalizationpoint ˜ s = κ − α .Physically, this implies that we measure (experimen-tally) the coupling at ˜ s and use it to determine whatvalue of the bare dimensionless coupling g correspondsto the physically measured running g R . Of course, oncethis measurement is made, the value of the coupling atany other momentum scale is determined.For this choice of measurement scale the connectionbetween g R and g simplifies to g R (˜ s ) = g ( κ )1 + g ( κ )(1 − d )(2 − d )(3 − d ) . (11) with g ( κ ) given by Eq. (9). Finally, combiningEqns. (10) and (11), the expression for the (running) cou-pling constant for arbitrary momentum and energy scaleis g R (˜ s ) = g R (˜ s )1 + g R (˜ s )(1 − d )(2 − d )(3 − d ) (cid:16)(cid:0) ˜ s + α κ (cid:1) − (1 − d ) − (cid:17) , (12)and of course by inspection, it is apparent that at therenormalization point ˜ s = κ − α one obtains the equiv-alence g R (˜ s ) = g R (˜ s ) as should be the case by defini-tion.We see that for d < s → ∞ we get g R → g . At the critical dimension d c = 1 we can expandthe denominator in a power-series to get g R (˜ s ) = g R (˜ s )1 − g R (˜ s )2 ln (cid:0) ˜ s + α κ (cid:1) , (13)which suggests that there is a Landau pole at g R (˜ s ) =2 / ln (cid:0) ˜ s + α κ (cid:1) ). A similar situation is found in quantumelectrodynamics where such a singularity is what leadsto the divergence of the bare charge e in the UV limit.Correspondingly, this is interpreted as the harbinger of new phenomena or degrees of freedom in the UV or,equivalently, at the shorter length scales, and which ina chemical context hints at the presence of short lived orintermediate substances in the mechanism of the originalchemical reaction.A naive use of this dimensionally regulated answeryields as d → g R (˜ s ) = (2 − d ) g R (˜ s )(2 − d ) − g R (˜ s ) (cid:0) ˜ s + α κ − (cid:1) , (14)implying that the coupling goes to zero in the contin-uum limit. This can also be seen through the associated -0.4-0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 ( g R ) g R g ?? g ? -0.4-0.2 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 g ? g ?? < d < d < d > < d < g R FIG. 3: The Beta function (15) as a function of the dimensionless coupling g R (˜ s ). Note that this refers to IR and not UVflow. There are two fixed points in the regions d < < d <
3, whereas there is a single fixed point when 1 < d < d > β function for g R (˜ s ) which has the particularly simpleform, β ( g R (˜ s )) = κ dg R (˜ s ) dκ = 2 g R (˜ s ) (cid:18) d − g R (˜ s )(2 − d )(3 − d ) (cid:19) . (15)This has two fixed points: a trivial one at g (cid:63)R = 0 anda non-trivial one at g (cid:63)(cid:63)R = (1 − d )(2 − d )(3 − d ). Notethat the non-trivial fixed point exists only when d < < d <
3. In the former case g (cid:63)(cid:63)R is IR stable while g (cid:63)R isUV stable, with the situation being reversed in the lattercase. On the other hand when 1 < d < d > g (cid:63)R which is IR stable but UVdivergent (see Fig. 3).Taken naively this seems to indicate that there is nomomentum dependence of the coupling once we are be-yond one dimension. In reality, however, there is a char-acteristic length scale in the system. Indeed, as discussedin the introduction, beneath the development of the mas-ter equation and its corresponding action (2) lurks theassumption of a lattice on which the chemicals hop be-tween sites. In addition, the size of the lattice must belarger than the mean free path of each chemical speciesin order to preserve the assumed Markovian nature ofthe collisions. In other words the term ( d −
2) in (14) re-ally indicates a cut-off of the form 1 / ln(Λ ), where Λ isthe maximum limit of the lattice momentum. Thereforein two dimensions, the action (2) represents an effective field theory with an UV cutoff.In three dimensions we get the dimensionally regulatedanswer: g R (˜ s ) = (3 − d ) g R (˜ s )(3 − d ) + g R (˜ s ) (cid:16)(cid:0) ˜ s + α κ (cid:1) − (cid:17) , (16) which suffers from the same pathologies as in the two-dimensional case.In order to obtain the correct effective theory, we makeuse of auxiliary fields that allows us to interpret the sumof loops in the vertex function Γ (3 , (i.e. the sum of elas-tic scatterings) as going through a single composite state σ . We can imagine the field σ as representing the densityof a “cloud” of chemicals involved in the elastic scatter-ing. This field then requires a “mass” (decay rate) andwave-function renormalization to render the system fi-nite. The effective theory can then be described in termsof some low energy parameters such as long distance re-action rates, as well as the low momentum decay ratesof composite fields. As is well known, this can be madeexplicit by making use of the well known technique of theHubbard-Stratonovich [18, 19] transformation, which wedescribe and apply next to this problem. III. COMPOSITE FIELD OPERATORS
We now introduce the composite fields σ and σ (cid:63) intothe previous field theory via a Hubbard-Stratonovichtransformation and discuss the differences between an“elementary” σ field and a composite one. We will findthat at large times and distances, these two theories areindistinguishable. However at short scales (high momen-tum), these theories differ for d = 1.Our starting point, once again, is the unshifted ac-tion (2) S = (cid:90) dx (cid:90) τ dt (cid:2) φ (cid:63)v ∂ t φ v + D v ∇ φ (cid:63)v ∇ φ v + φ (cid:63)u ∂ t φ u + D u ∇ φ (cid:63)u ∇ φ u + µ ( φ (cid:63)v − φ v + ν ( φ (cid:63)u − φ u − f ( φ (cid:63)u − − λ ( φ (cid:63)v − φ (cid:63)u ) φ (cid:63) v φ v φ u (cid:3) . (17) u v v ˜ ˜ v ˜ v ˜ u (b) u v v ˜ ˜ v ˜ v (a) ˜ v FIG. 4: The vertex function Γ (3 , shown in Fig. 2 reinterpreted as the production and decay of the composite field σ . Thecomposite field σ is formed out of the combination λφ u φ v and leads to the two decay channels shown in (a) and (b). To carry out the Hubbard-Stratonovich transformation,we construct a second action which defines the compositefields σ and σ (cid:63) thus, S (cid:48) = (cid:90) dx (cid:90) τ dtλ (cid:2)(cid:0) ( σ (cid:63) − − ( φ (cid:63)v − φ (cid:63)u ) φ (cid:63) v (cid:1) × ( σ − φ v φ u ) (cid:3) , (18)and add it to S obtaining, S comp = S + S (cid:48) = (cid:90) dx (cid:90) τ dt (cid:2) φ (cid:63)v ∂ t φ v + D v ∇ φ (cid:63)v ∇ φ v + φ (cid:63)u ∂ t φ u , + D u ∇ φ (cid:63)u ∇ φ u + µ ( φ (cid:63)v − φ v + ν ( φ (cid:63)u − φ u − f ( φ (cid:63)u − − λ σ ( φ (cid:63)v − φ (cid:63)u ) φ (cid:63) v − λ ( σ (cid:63) − φ v φ u + λ ( σ (cid:63) − σ (cid:3) . (19)Eliminating the constraint equation defining σ leads backto the original equation of motion for the fields φ i . From(19), it is apparent that while the composite field σ isformed through a combination of a single φ u and two φ v ’s, it has two potential decay channels (i) it convertsinto either 3 φ v ’s or (ii) back again to the original con-stituents, a single φ u and two φ v ’s (Cf. Fig. 4).It is instructive to perform a comparison of the above,with the case where we consider σ as an elementary scalarparticle. In this version, instead of the local interaction λ ( σ (cid:63) − σ found in S comp there is a fundamental (orrather bare) kinetic energy term for the σ field, as wellas an unrenormalized decay rate M . To keep the dimen-sionality of the elementary σ the same as in the compos-ite case, the “free” part of the σ field action needs to bedivided by κ d . Consequently, the action is now S elem = (cid:90) dx (cid:90) τ dt (cid:2) φ (cid:63)v ∂ t φ v + D v ∇ φ (cid:63)v ∇ φ v + φ (cid:63)u ∂ t φ u + D u ∇ φ (cid:63)u ∇ φ u + ( σ (cid:63) ∂ t σ + D σ ∇ σ (cid:63) ∇ σ ) κ − d + µ ( φ (cid:63)v − φ v + ν ( φ (cid:63)u − φ u − f ( φ (cid:63)u − M κ − d ( σ (cid:63) − σ − λ σ ( φ (cid:63)v − φ (cid:63)u ) φ (cid:63) v − λ ( σ (cid:63) − φ v φ u (cid:3) . (20) Note that in the models defined by either (19) or (20),the structure of the Schwinger-Dyson equations are quitesimilar.Having introduced σ , the process of elastic scatteringis now interpreted as proceeding through the exchange ofa composite field. This involves the composite propaga-tor, whose inverse is given by the dimensionally regulatedexpression G − (˜ s ) = g ( κ ) + g ( κ )(1 − d )(2 − d )(3 − d ) × (cid:18) ˜ s + α κ (cid:19) − (1 − d ) . (21)In the case where σ is an elementary particle, this isinstead G − (˜ s ) = ˜ s + M κ + g ( κ )(1 − d )(2 − d )(3 − d ) × (cid:18) ˜ s + α κ (cid:19) − (1 − d ) , (22)where for the sake of simplicity we have chosen the samediffusion constant for the free part as in the induced one,i.e. D σ = ˜ D .To complete the renormalization procedure in termsof the composite field σ we must now allow for “wavefunction” renormalization for the σ field as well as “mass”(decay rate) renormalization. Since the vertex functionΓ (3 , depends on the combination ˜ s = s + ˜ D p , theinverse propagator (21) can therefore be expanded in apower series in ˜ s thus G − = g ( κ ) + g ( κ ) F (˜ s ) ≡ g ( κ ) + g ( κ )(1 − d )(2 − d )(3 − d ) (cid:18) α + ˜ sκ (cid:19) − (1 − d ) = g ( κ ) + g ( κ )(1 − d )(2 − d )(3 − d ) (cid:16) α κ (cid:17) − (1 − d ) × (cid:20) d − sα + ( d − d − s α + ( d − d − d − s α + . . . (cid:21) , (23)where F (˜ s ) can be thought of as a self-energy function.The expansion for the elementary σ is identical, with theexception that the first term in (23) i.e g ( κ ), is replacedby ( M + ˜ s ) /κ . Otherwise the renormalization proce-dure is identical in both cases. Let us first consider thesituation for d < A. One and two dimensions
For both the composite and elementary σ ’s we canrewrite G − in the form G − = Z − (cid:20) M σ + ˜ sκ + g ( κ ) Z F sub (˜ s ) (cid:21) , (24)for d < Z and M σ are obtained fromthe first two terms in the expansion of Eq. (23). (We willhenceforth drop the subscript for G and unless mentionedotherwise it refers to both the elementary and compositeversions of the model.) The subscript sub refers to thesubtraction of the first two terms in the power series (23), M σ is the renormalized decay rate for the σ field and Z is the wave function renormalization. Note that at lowmomentum, apart from the rescaling by Z − , the inversepropagator resembles the free field one.When σ is composite this leads to the identities Z − M σ κ = g ( κ ) + g ( κ )(1 − d )(2 − d )(3 − d ) (cid:16) α κ (cid:17) − (1 − d ) ,Z − ˜ sκ = − g ( κ )(2 − d )(3 − d ) (cid:16) α κ (cid:17) − (2 − d ) ˜ sκ . (25)whereas when it is considered elementary we have, Z − M σ κ = M κ + g ( κ )(1 − d )(2 − d )(3 − d ) (cid:16) α κ (cid:17) − (1 − d ) ,Z − ˜ sκ = (cid:18) − g ( κ )(2 − d )(3 − d ) (cid:16) α κ (cid:17) − (2 − d ) (cid:19) ˜ sκ . (26)Next we introduce the renormalized vertex g σ for σ → φ u φ v making use of the identities g σ = g ( κ ) Z ; G R = Z − G, (27) where G R is the renormalized propagator for σ . Conse-quently the combination g ( κ ) G = g σ G R (28)is invariant under renormalization. Combining theseleads us to a finite expression for the running reactionrate thus g R (˜ s ) = g σ G R = g σM σ +˜ sκ + g σ F sub (˜ s ) . (29)In two dimensions F sub is just zero and therefore thisreduces to g R (˜ s ) = g σ κ M σ + ˜ s . (30)In particular using (28) and the definition of G R we havethat g R (˜ s = 0) = g σ κ M σ = g ( κ ) κ Z − M σ . (31)Finally, employing Eqns. (25) or (26) allows us to relate g R (0) to g ( κ ).We thus conclude that in d = 2, the renormalized cou-pling g R (˜ s ) is proportional to d − Z . This suggests that whenwe use actual cutoffs to regulate the integrals, the renor-malized coupling is related to the inverse of the phys-ical cutoff of the problem. This confirms our previousintuition that the system can only be described by an effective field theory with a momentum space (or spa-tial) cutoff. It is also apparent that there is no differencebetween a fundamental or composite σ at this level ineither the IR or UV limits. The dependence on the pa-rameter M σ can be eliminated by evaluating the runningcoupling constant (momentum dependent reaction rate)at a particular reference point ˜ s . In terms of this ref-erence point ˜ s = κ − α (as chosen for Eq. (11)) onefinds the reaction rate is given by g R (˜ s ) = g R (˜ s )1 + g R (˜ s ) (˜ s − ˜ s ) κ (32)thus explicitly showing that M σ can be replaced by thephysical quantity g R (˜ s ).Moving onto one dimension which is the critical di-mension, from Eqns. (25) and (26) we find that the wavefunction renormalization is finite. Thus only mass renor-malization is needed which can again be translated intodefining g R (˜ s ) at a particular reference value ˜ s . In thiscase the second term in (23) is g ( κ ) F (˜ s ) = g ( κ )2(1 − d ) − g ( κ )2 ln ˜ s + α κ , (33)which leads to G − = Z − (cid:20) M σ + ˜ sκ − g σ (cid:18) ln (cid:18) ˜ s + α α (cid:19) − ˜ sα (cid:19)(cid:21) . (34)Once again making use of (27) we find that g σ = − α /κ and therefore the terms linear in ˜ s cancel. Thusthe renormalized propagator for the composite σ is G − R = M σ κ − g σ (cid:18) sα (cid:19) , (35)Once again the parameter M σ can be eliminated by defin-ing g R (˜ s ) at some scale ˜ s such that g R (˜ s ) = g R (˜ s )1 − g R (˜ s )2 ln (cid:16) ˜ s + α ˜ s + α (cid:17) (36)implying that it goes to zero as (ln ˜ s ) − in the UV limit.Note that this is equivalent to our previous result Eq.(13).Through a similar sequence of manipulations it can beshown that for the elementary σ we have the relation g σ = Z g ( κ ) = g ( κ )1 − g ( κ ) κ α (37)and therefore this is equivalent to the composite σ onlywhen g → ∞ in which case Z = 0 (the standard “com-positeness” condition in field theory). The correspondinganalog to Eq. (36) is g R (˜ s ) = g R (˜ s )1 − g R (˜ s )2 ln (cid:16) ˜ s + α ˜ s + α (cid:17) + g R (˜ s ) (˜ s − ˜ s ) f ( α , g σ , κ )(38)where f ( α , g σ , κ ) = α + k g σ α k g σ . Here, the running cou-pling goes to zero linearly with ˜ s as opposed to logarith-mically. Thus unlike in the two dimensional case there isa difference between the elementary and composite man-ifestations of σ . B. Higher dimensions ( d ≥ The renormalization process for d ≥ σ . This is due to the fact that the self-energy function F (˜ s ) in (23) has an increasing number ofdivergent terms.In one and two dimensions, as discussed previously,the first two terms in the series diverge and are regulatedvia decay rate and wave function renormalization respec-tively. Starting from d = 3 one begins to induce termsthat are not originally present in Eqns. (19) and (20).Consequently one needs to subtract three terms from F (˜ s ) to obtain a finite contribution. (And four in d = 4and so on.) The resulting renormalized propagator canbe written as˜ g R G R = ˜ g RM σ +˜ sκ +
12 ˜ s κ + ˜ g r F sub (˜ s ) , (39)where the term
12 ˜ s κ corresponds to the addition of a newinduced term that needs to be inserted in (19) and (20).This term has the form S ind = (cid:90) dt (cid:90) dx (cid:104) σ (cid:63) ( x, t ) (cid:0) ∂ t − D σ ∇ (cid:1) σ ( x, t ) (cid:105) κ − d +1) . (40) IV. FLUCTUATION INDUCED PROCESSES
The interactions shown in Fig. (1) are such that theycan induce processes (via fluctuations) that are not ex-plicitly present in the original action described by Eq.(2) . (For examples of this phenomenon in out of equilib-rium RD systems, but generated by extrinsic noise, cf.Ref. [8].) While this is not problematic if the inducedterms are f inite , a serious problem appears if the termsbring with them a divergence. If each order in perturba-tion theory leads to new interactions which are divergent,then usually an infinite number of parameters are neededto define the theory. In such a situation physical predic-tions cannot be made as they depend on an arbitrarynumber of parameters. In that case the theory is callednon-renormalizable. Although we will find that indeedthere are an infinite number of induced processes thatappear superficially divergent, by an appropriate intro-duction of composite field operators, we will find thatonly three long distance (infrared) parameters will beneeded to describe all the induced processes.The induced processes in this case can be constructedfrom diagrams that represent n → n scattering for n ≥ n = 4 , U + 2 V → U + 3 V in (a) and 3 U + 2 V → U + 3 V in (b). Although at firstsight it appears these diagrams are naively divergent, ifone thinks of them as proceeding through “tree level”diagrams in terms of the production of di-molecule andtri-molecule states, then the only renormalizations neces-sary are those of the tri-molecule σ propagator as well asdecay constant renormalization for the di-molecule corre-lation functions. In terms of the renormalized parameters u v v ˜ v ˜ v ˜ v u u u u u v v ˜ v ˜ v ˜ v u u (a) (b) u v v ˜ v ˜ v ˜ v (c) v u v v ˜ v ˜ v ˜ v u u (d) v v v FIG. 5: Fluctuation induced n → n scattering for (a) n = 4; 2 U + 2 V → U + 3 V and (b) n = 5; 3 U + 2 V → U + 3 V .The shaded circle represents the renormalized original 3 → φ v ’s. which can be thought of as a composite field correlation function for ψ = φ v . We can haverelated scattering processes through the composite field correlation function for ψ = φ u φ v as shown in (c) for the reaction U + 3 V → V and indeed a combination of the two composite field correlation functions as shown in (d) for the reaction2 U + 3 V → U + 4 V . of the σ propagator and di-molecule propagators, theseinduced processes are then rendered finite.Going back to Fig. 5a,b, we see the formation of a sin-gle loop which now consists of two φ v ’s (or, equivalently,a correlation function for the composite field ψ = φ v ).If we instead considered the process U + 3 V → V (c)this would have been facilitated instead by the correla-tion function for the composite field ψ = φ u φ v . In termsof the composite fields, these processes are “tree graphs”in the Green’s functions for σ , ψ and ψ . These internalprocesses are reversible in that φ u + φ v → ψ → φ u + φ v and 2 φ v → ψ → φ v etc. Starting at n = 5 we canhave both the ψ and ψ correlation functions occurringin the composite field “tree” diagram as shown in (d) forthe process 2 U + 3 V → U + 4 V .To expose and study the divergence structure of the ψ , correlation functions, we only need to calculate asingle loop (unlike for the case of the elementary σ cor-relation function which is a geometric sum of two loopgraphs). We will assume that an external momentum p flows into the right-most vertices in Fig. 5 and to simplifythe discussion we will set the energy and momentum ofthe “new” incoming and outgoing particles (beyond thebasic U + 2 V → V reaction) to zero. We will then onlyhave to consider the internal one loop graphs at somearbitrary momentum p . The expression for the one loopintegral I ( p, t ) is given by, I ( p, t ) = 2! (cid:90) (cid:89) i =1 (cid:18) d d p i (2 π ) d (cid:19) e − [ D v ( p + p )+2 µ ] t × (2 π ) d δ ( p − p + p ) . (41)This can be calculated exactly and after Laplace trans- forming to s –space one gets, I ( p, s ) = B D − d/d c Γ(1 + (cid:15)/d c ) (cid:15)/d c (cid:18) s + α + D p (cid:19) − (cid:15)/d c , (42)where now B = 2!(4 π ) d/d c , D = 2 D v , α = 2 µ. (43)Note that for the single loop the critical dimension is now d c = 2, and therefore (cid:15) = 2 − d . At d = 1, (cid:15) = 1 and weget the finite result: I ( p, t ) = B D − / (cid:18) s + α + D p (cid:19) − / , (44)so that this process vanishes as p − for large momenta.To evaluate I ( p, s ) in the critical dimension d c = 2,we exponentiate the term in s above and expand theexponential thus, (cid:18) s + α + D p (cid:19) − (cid:15)/ = e − ( (cid:15)/
2) ln( s (cid:48) ) = 1 − (cid:15) s (cid:48) + O( (cid:15) ) . (45)Inserting the result of this expansion back into (42), weimmediately see that the only term that diverges as (cid:15) → I corresponds to the correlation function for the composite0field ψ . Next identifying the zero energy-momentumpiece of the dimensionless propagator as α /M we get I (0) ≡ α M = C (cid:18) (cid:15) −
12 ln α (cid:19) . (46)Consequently the renormalized version of the correlationfunction is now ˜ I ( p, s ) = α M − C s (cid:48) α (47)where C = 2 B D − .Thus in two dimensions, all n → n processes are ren-dered finite by the introduction of only two new parame-ters M , M corresponding to the decay rates for the twocomposite fields ψ and ψ . In terms of these two param-eters (and the renormalized coupling g R ) all the inducedcouplings can be determined.Finally, the full expression for the Laplace transform ofthe 4 → g R (˜ s ) I ( p, s ),where at d c = 2 we need to replace I by the regulated˜ I . It is not hard to see that for processes like Fig. 5band extensions to the form nU + 2 V → ( n − U + 3 V this generalizes to g R (˜ s ) n ˜ I ( p, s ) n − . (48) V. CONCLUSIONS
In this paper we presented a detailed analysis of theeffects of intrinsic noise on a spatially extended reaction-diffusion chemical model (1). We found that remarkably,the short distance behavior of the system could be de-termined analytically by studying the model in the U (1)symmetric phase of the field theory [9] corresponding tothe chemical reactions described by Eq. (1). The sym-metric phase of the field theory reflects particle numberconservation and consequently the only allowed graphs inits Feynman diagram representation are for n → n , with n ≥
3. The fluctuations due to the intrinsic noise leads totwo types of potentially divergent graphs in the theory.The first divergent graph, is a 2-loop graph which firstdiverges in the critical dimension d c = 1. This graphwe relate to the renormalization of the 3 → d = 2 in the induced n → n processes. These we regulateusing the standard technique of dimensional regulariza-tion. We then investigated what (if any) are the new lowenergy (large spatio-temporal) parameters that need tobe specified to define the correct finite and renormalizedtheory which includes the effects of the fluctuations.We find that one parameter is the critical dimensionfor the behavior of the reaction rate parameter λ , whichhappens to be d c = 1. This reaction rate gets renor-malized through a sequence of elastic collisions (two-loopgraphs) that occur between the chemically relevant in-elastic collision. As a result, it acquires a momentum de-pendence, which in the critical dimension we can specify by determining the reaction rate at late times. Equiva-lently, this is also determined by representing the elasticcollisions as a composite three body state σ = φ u φ v , andthen determining its rate of decay. We also find that in d c = 1 it is possible to distinguish between the situa-tions where σ is a bonafide composite state and the casewhere instead, it is an elementary chemical with its ownkinetic energy and decay terms. This is done by investi-gating the large momentum (short distance) behavior ofthe momentum dependent reaction rate. The decay ratein the version with an elementary σ goes to zero in thehigh momentum UV limit faster than in the case where σ is composite. The point where the wave function renor-malization of the elementary σ goes to zero, is where thetwo versions of the model are identical.Starting in two dimensions, two new parameters areneeded to describe the system. These can be thoughtof as being the decay rates ( masses ) for the compositesystems φ v and φ u φ v . To obtain these parameters onewould need to experimentally measure the reaction ratesfor two inelastic reactions such as 2 U +2 V → U +3 V and3 U + 2 V → U + 3 V . Once that is done, all the induced n → n reaction rates can be calculated. The renormal-ized equation for the σ field in two dimensions includesan induced free field kinetic term. In 3 dimensions thefluctuations have a further effect of changing the renor-malized equation for the σ field to one having higherspatio-temporal derivatives. Thus in terms of the run-ning coupling constant g R (˜ s ) as well as two measurableinduced coupling constants, we have determined the ef-fective field theory which results from the intrinsic noiseinherent in this chemical reaction diffusion model.We did not discuss the infrared properties of the theorywith broken symmetry, which is the sector that relates di-rectly to the chemistry. To do so, one would follow a twostep approach. First, one needs to determine the classicaldensities as a function of time and the classical responsefunction which depends on both momentum and time.Then we would make use of the running coupling con-stant found in this paper to determine the asymptoticbehavior of the momentum dependent densities includ-ing fluctuations by solving the Callan-Symanzik equa-tion [21, 22]. This approach is worked out in detail forthe example of the process 3 A → A in [11] and extendsthe arguments standard in relativistic quantum field the-ory to this class of non-equilibrium models. Appendix A: Chemical master equation and manybody formalism
In order to develop the master equation formalismfor our system of chemical reactions (1), we first di-vide the space in which the reactions take place into a d − dimensional hyper-cubic lattice of cells and assumethat we can treat each cell as a coherent entity. We as-sume the interactions occur locally at a single cell site andthat there is also diffusion modeled as hopping between1nearest neighbors. Assuming that the underlying pro-cesses are Markovian, they can be described by a prob-ability distribution function P ( n v , n u , t ) which gives theprobability to find the particle configuration ( n v , n u ) attime t . Here n i ( t ) = ( { n i ( t ) } ) is a vector composition variable where n i represents the number of molecules ofa species at site i . Following standard methods, one ob-tains for the master equation for the chemical reactionsin (1) including diffusion ddt P ( n v , n u , t ) = D v l (cid:88) (cid:104) i,j (cid:105) [( n v,j + 1) P ( . . . , n v,i − , n v,j + 1 , . . . , t ) − n v,i P ]+ D u l (cid:88) (cid:104) i,j (cid:105) [( n u,j + 1) P ( . . . , n u,i − , n u,j + 1 , . . . , t ) − n u,i P ]+ λ (cid:88) i [( n v,i − n v,i − n u,i + 1) P ( . . . , n v,i − , . . . , n u,i + 1 , . . . , t ) − n v,i ( n v,i − P ]+ µ (cid:88) i [( n v,i + 1) P ( . . . , n v,i + 1 , . . . , t ) − n v,i P ] + ν (cid:88) i [( n u,i + 1) P ( . . . , n u,i + 1 , . . . , t ) − n u,i P ]+ f (cid:88) i [ P ( . . . , n v,i + 1 , . . . , t ) − P ] , (A1)where l is the characteristic length of the cell and (cid:104) . . . (cid:105) denotes the sum over nearest neighbors.The master equation (A1) along with the sextic inter-action shown in Fig. 1 lends itself to a many body de-scription [13], accomplished by the introduction of an oc-cupation number algebra with annihilation/creation op-erators ˆ a i , ˆ a † i for v and ˆ b i , ˆ b † i for u at each site i . Theseoperators obey the Bosonic commutation relations (cid:104) ˆ a i , ˆ a † j (cid:105) = δ ij , (cid:104) ˆ b i , ˆ b † j (cid:105) = δ ij , [ˆ a i , ˆ a j ] = 0 , (cid:104) ˆ a † i , ˆ a † j (cid:105) = 0 , (A2)and define the occupation number operators ˆ n i,v = ˆ a † i ˆ a i and ˆ n i,u = ˆ b † i b i satisfying the following eigenvalue equa-tions:ˆ n i,v | n i,v (cid:105) = n i,v | n i,v (cid:105) , ˆ n i,u | n i,u (cid:105) = n i,u | n i,u (cid:105) . (A3)We next construct the state vector | Ψ( t ) (cid:105) = (cid:88) n v , n u P ( n v , n u , t ) × (cid:89) i (ˆ a † i ) n iv (ˆ b † i ) n iu | (cid:105) , (A4)which upon differentiating with respect to time t , can bewritten in the suggestive form − ∂ | Ψ( t ) (cid:105) ∂t = H [ ˆa † , ˆa , ˆb † , b ] | ψ ( t ) (cid:105) , (A5)resembling the Schr¨odinger equation. Finally, taking thetime derivative of Eq. (A4) and comparing terms with the Hamiltonian in (A5) we make the identification H = D v l (cid:88) (cid:104) i,j (cid:105) (ˆ a † i − ˆ a † j )(ˆ a i − ˆ a j ) + µ (cid:88) i (ˆ a † i − a i + D u l (cid:88) (cid:104) i,j (cid:105) (ˆ b † i − ˆ b † i )(ˆ b i − ˆ b j ) + ν (cid:88) i (ˆ b † i − b i − λ (cid:88) i (cid:2) ˆ a † i − ˆ a † i ˆ b † i (cid:3) ˆ a i ˆ b i − f (cid:88) i (ˆ b † i − . (A6)Having defined the space, the appropriate wave functionand the Hamiltonian, we next seek to evaluate the oper-ator e − ˜ Ht using the path integral formulation. Followingthe standard procedure for obtaining the coherent statepath integral [12, 15, 16] to the GS system, letting thecoherent state φ v (related to the operator a ) represent v and φ u (related to the operator b ) represent u we obtain e − ˜ Ht = (cid:90) D φ v D φ (cid:63)v D φ u D φ (cid:63)u e − S [ φ v ,φ (cid:63)v ,φ u ,φ (cid:63)u ] , (A7)where the action S is given by S = (cid:90) dx (cid:90) τ dt (cid:2) φ (cid:63)v ∂ t φ v + D v ∇ φ (cid:63)v ∇ φ v + φ (cid:63)u ∂ t φ u + D u ∇ φ (cid:63)u ∇ φ u + µ ( φ (cid:63)v − φ v + ν ( φ (cid:63)u − φ u − f ( φ (cid:63)u − − λ φ (cid:63)v − φ (cid:63)u ) φ (cid:63) v φ v φ u (cid:3) . (A8)This is Eq. (2) in the body of the text.2 [1] M. Cross and H. Greenside, Pattern Formation and Dy-namics in Nonequilibrium Systems (Cambridge Univ.Press, Cambridge, 2009).[2] D. Walgraef,
Spatio-Temporal Pattern Formation (Springer, New York, 1997).[3] A. S. Mikhailov,
Foundations of Synergetics I , 2 ed.(Springer, Berlin, 1994).[4] B. A. Grzybowski,
Chemistry in Motion: Reaction-Diffusion Systems for Micro- and Nanotechnology (Wi-ley, Chichester, 2009).[5] J. E. Pearson, Science , 189 (1993).[6] K.-J. Lee, W. D. McCormick, J. E. Pearson, and H. L.Swinney, Nature , 215 (1994).[7] F. Lesmes, D. Hochberg, F. Mor´an, and J. P´erez-Mercader, Phys. Rev. Lett. , 238301 (2003).[8] D. Hochberg, F. Lesmes, F. Moran, and J. P´erez-Mercader, Phys. Rev. E. , 066114 (2003).[9] F. Cooper, G. Ghoshal, and J. P´erez-Mercader, Phys.Rev. E , 042926 (2013).[10] F. Cooper, G. Ghoshal, A. Pawling, and J. P´erez-Mercader, Phys. Rev. Lett. , 044101 (2013). [11] B. P. Vollmayr-Lee, J. Phys. A: Math. Gen. , 2633(1994).[12] U. C. T¨auber, M. Howard, and B. P. Vollmayr-Lee, J.Phys. A: Math. Gen. , R79 (2005).[13] M. Doi, J. Phys. A: Math. Gen. , 1465 (1976).[14] P. Grassberger and M. Scheunert, Fortschritte der Physik , 547 (1980).[15] J. W. Negele and H. Orland, Quantum Many-ParticleSystems (Perseus Publishing, Cambridge, 1998).[16] U. C. T¨auber, in
Aging and the Glass Transition , Vol. 716of
Springer Lecture Notes in Physics , edited by M.Henkel, M. Pleimling, and R. Sanctuary (Springer-Verlag, Berlin, 2007), pp. 295–348.[17] M.-P. Zorzano, D. Hochberg, and F. Moran, Phys. Rev.E , 057102 (2006).[18] R. Stratonovich, Doklady , 416 (1958).[19] J. Hubbard, Phys. Rev. Lett. , 77 (1959).[20] M. Gell-Mann and F. Low, Phys. Rev. , 350 (1951).[21] C. G. Callan, Phys. Rev. D , 1541 (1970).[22] K. Symanzik, Commun. Math. Phys.18