{'fraction_non_alphanumeric': 0.05048648163924136, 'fraction_numerical': 0.01624220956821952, 'mean_word_length': 4.151163462093654, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Field Theory of Branching and Annihilating Random Walks

We develop a systematic analytic approach to the problem of branching and annihilating random walks, equivalent to the diffusion-limited reaction processes 2A->0 and A->(m+1)A, where m>=1. Starting from the master equation, a field-theoretic representation of the problem is derived, and fluctuation effects are taken into account via diagrammatic and renormalization group methods. For d>2, the mean-field rate equation, which predicts an active phase as soon as the branching process is switched on, applies qualitatively for both even and odd m, but the behavior in lower dimensions is shown to be quite different for these two cases. For even m, and d~2, the active phase still appears immediately, but with non-trivial crossover exponents which we compute in an expansion in eps=2-d, and with logarithmic corrections in d=2. However, there exists a second critical dimension d_c'~4/3 below which a non-trivial inactive phase emerges, with asymptotic behavior characteristic of the pure annihilation process. This is confirmed by an exact calculation in d=1. The subsequent transition to the active phase, which represents a new non-trivial dynamic universality class, is then investigated within a truncated loop expansion. For odd m, we show that the fluctuations of the annihilation process are strong enough to create a non-trivial inactive phase for all d<=2. In this case, the transition to the active phase is in the directed percolation universality class.


I. INTRODUCTION
Reaction-diffusion processes, in addition to being of direct physical interest, present a relatively simple class of nonequilibrium stochastic processes in which the role of fluctuations, so important to the study of equilibrium statistical physics, may be investigated [1]. Although such processes have been studied for many years through various kinds of mean-field or self-consistent approximations, more recently there has been some progress in taking into account the effects of fluctuations through the use either of exact methods or of systematic renormalization group theories [2], backed up by extensive numerical simulations.
It has become clear that, once the constraints of detailed balance are abandoned, these systems may exhibit a greater richness of phenomena than their counterparts close to equilibrium, leading, for example, to new types of critical behavior characterized by new universality classes. It might appear that the freedom allowed by giving up detailed balance would lead to such an unacceptably large number of possibilities as to make the concept of universality of limited usefulness, and indeed, from the point of view of the proliferation of microscopic cellular automaton models this would seem to be true. However, as with equilibrium phenomena, the introduction of a coarse-grained description through a field-theoretic formalism makes it possible, at least in principle, to make a universal classification feasible.
All such coarse-grained descriptions of stochastic processes, whether or not they satisfy detailed balance, divide the dynamics into two parts: a deterministic evolution corresponding to the mean-field dynamics, and the noisy effect of fluctuations ignored in such a description. When detailed balance is present, the form of the noise is prescribed, so that the field-theoretic version of the coarse-grained description is narrowly determined. In the absence of detailed balance this is not so, yet it may be shown that the particular form of the noise chosen is crucial to determining the universal behavior near a non-equilibrium critical point. One way to avoid this ambiguity, which works very well for reaction-diffusion processes, is to develop a field-theoretic description of the problem which is exact, being directly derived from the underlying master equation, and to perform coarse-graining, in terms of a derivative expansion, only at a later stage. This method, developed some time ago by Doi and by Peliti [3,4] (see also Ref. [5]), relies on the use of a Fock space formalism very similar to that employed for quantum many-body systems. More recently, field-theoretic renormalization group methods have been systematically applied to the dynamic field theories describing some simple processes to obtain what is, in principle, a complete characterization of their universal properties [6][7][8].
It is the purpose of this paper to describe such a field-theoretic investigation, largely but not entirely based on renormalization group methods, into a class of reaction-diffusion processes which are characterized by their description in terms of branching and annihilating random walks. (A very brief account of this work has been presented in Ref. [9].) These processes were originally studied in the 1980s by Grassberger et al. [10] as exceptions to the apparent rule that systems exhibiting non-equilibrium transitions from a trivial absorbing state into a non-trivial noisy 'active' state should be in the same universality class as directed percolation (DP). While this terminology originally refers to a realization of this universality class describing percolation processes which are asymmetric with respect to a privileged space dimension [12], more commonly it now applies to a dynamic universality class in which this dimension is interpreted as time. In terms of reaction-diffusion processes, an example of such a system would be one with a single species of particle A, undergoing diffusive behavior, single-particle annihilation A → ∅, and branching A → 2A. There is always a trivial absorbing state, with no particles and consequently no fluctuations. For sufficiently low branching rate, this is the only stationary state, but for larger values of this rate, another, non-trivial 'active' stationary state appears. The DP transition occurs at the point where this state first appears. The robustness of the DP universality class has been tested by measuring its critical exponents in numerical simulations on many different types of system, and may be understood in terms of the simplicity and irreducibility of the dynamic field theory which describes it [13]. Indeed, this 'Reggeon' field theory was first studied in a completely different context of particle physics. Among the non-equilibrium phase transitions described by this universality class are those occurring in the contact process [14], the dimer poisoning problem in the Ziff-Gulari-Barshad model [15], and in certain autocatalytic reaction models [16].
The model studied by Grassberger et al. [10] was a rather complicated probabilistic cellular automaton in one space and one time dimension, with antiferromagnetic interactions in which the 'particles' actually appear as kinks or domain walls. It was observed that, while this model appears to exhibit a transition from an inactive, absorbing state to a non-trivial active state, the various critical exponents are far from those of directed percolation. The authors of Ref. [10] already realized that this may be a result of the fact that these kinks are conserved modulo 2 in the basic dynamics, but a subsequent field-theoretic analysis of this problem appeared to shed no light on this discrepancy: rather it suggested that the system should, once again, be in the DP universality class, indicating that there was apparently no way of taking into account the constraint of local 'parity' conservation in a field-theoretic treatment [11]. One of the main points of our analysis is to show where this argument went wrong, and why such systems, at least in low dimensions, should be in a new universality class.
Subsequent work showed that similar effects occur in other systems with such a conservation law, and suggested that, at least in one dimension, they all fall into this new universality class, characterized by branching and annihilating walks with an even number of offspring [17][18][19][20]. These are reaction-diffusion systems with the underlying reaction processes 2A → ∅ and A → (m + 1)A, with m even. In one dimension, when m = 2, another realization of this universality class occurs in certain classes of dynamic Ising models which violate detailed balance [10,11,21], where once again the 'particles' correspond to domain walls. In this case, the transition from the inactive phase (no domain walls in the asymptotic stationary state) to the active phase corresponds to the disordering transition of the Ising spins [22,23]. In addition, other realizations of this universality class have been found and also been investigated numerically [24,25].
On the other hand, when m is odd and there is no conservation law, one expects to recover the DP exponents, and this appears to be the case, at least in one dimension [10,17,26,27]. In higher dimensions, the numerical situation is less clear. Takayasu and Tretyakov, in rather short simulations, found no evidence for a non-trivial transition for d ≥ 2, for either m even or odd [17]. This is certainly consistent with the picture obtained from mean-field theory, as we shall see. But the situation in two dimensions is, we believe, more subtle. It is just this dimension at which the fluctuation effects become important in the pure annihilation process 2A → ∅, giving logarithmic corrections to scaling, and they cannot therefore be ignored there. In fact we find in our analysis that these fluctuation effects are sufficient to drive the existence of an inactive phase for small but finite values of the branching rate, at least for odd m, and hence the existence of a subsequent transition in the active phase, which is in the two-dimensional DP universality class.
For m even, on the other hand, the fluctuations are not effective until a second, lower critical dimension which we estimate to be ≈ 4 3 , and which is certainly above d = 1. The existence of two critical dimensionalities make a systematic ǫ-expansion type of approach to this problem infeasible. However, we have performed a truncated loop expansion which has the correct fixed point structure to describe the new d = 1 universality class. One striking feature of this is the existence of dangerous irrelevant variables which make the identification of the various critical exponents in terms of RG eigenvalues slightly delicate. In fact, we may only make this identification strictly in the active phase, and our analysis leaves open the possibility of violations of the scaling relations between these exponents and those at the critical point and in the inactive phase, which have been written down on the basis of simple scaling arguments.
Our aim, therefore, is to investigate the following competing reaction processes of identical particles A, which perform a random walk with diffusion constant D: with k, m integers, k ≥ 2, and m ≥ 0; '∅' denotes an inert state. Occasionally, we shall in addition allow for a spontaneous creation of particles in pairs, The corresponding mean-field rate equation for the average density n(t), neglecting any local fluctuation effects, reads with an initial density n(0) = n 0 . If we set the pair production rate, which acts as an 'external field', to zero, Eq. (1.4) has the stationary solutions n = 0 ('inactive state') and ('active state'). Provided the 'mass' mσ m is positive (and n 0 > 0), the solution of Eq. (1.4) approaches the asymptotic density n s exponentially for t → ∞, On the other hand, if σ m = 0, this exponential long-time behavior is replaced by a power-law decay, In this sense, σ m = σ c = 0 may be viewed as a critical point, described by the following set of critical exponents characterizing the long-time decay at the critical point, and the growth of the asymptotic density in the active phase, respectively; according to Eqs. (1.5) and (1.7) their mean-field values are α 0 = β 0 = 1/(k − 1). Taking into account the diffusive propagation of local density fluctuations, two more independent exponents may be introduced which describe the divergence of the characteristic length and time scales as the critical point is approached, with the mean-field values ν 0 = 1/2, ν τ 0 = 1/(d + 2), and z 0 = 2. These mean-field exponents of course originate in the scaling dimensions of the parameters D, λ k , σ m , and τ of the model. Introducing a momentum scale κ, and measuring times in terms of lengths squared, (1.15) Therefore, both particle production processes appear to be relevant perturbations to the pure annihilation reaction, and for σ m = τ = 0 the annihilation rate becomes dimensionless at the (upper) critical dimension d c = 2/(k − 1). Thus, one expects the mean-field rate equation (1.4) to provide an at least qualitatively correct description for d > 2 in the case k = 2, and for d > 1 in the case k = 3. Because of the relevant perturbation σ m , the situation at the critical dimension already requires a careful analysis, and certainly for d < d c fluctuation effects need to be taken into account seriously. E.g., for the pure annihilation process with k = 2 it is known exactly that Eq. (1.7) is replaced by the slower asymptotic decay 16) due to strong particle anticorrelations which cannot be smoothened out by the diffusion in low dimensions [3,4,6]. At d = 2 for k = 2 and d = 1 for k = 3, one finds logarithmic corrections to Eq. (1.7). However, for k ≥ 4 the above rate equation should be essentially correct in all physical dimensions. In the bulk of this paper, we shall therefore be concerned with the case k = 2, and only briefly discuss the three-particle annihilation reaction (k = 3) combined with branching in Sec. VIII. The outline of the rest of this paper is as follows. After introducing the field-theoretic description of these processes, with an emphasis on the subtleties which occur when correctly taking into account processes in which a conservation law is satisfied, in Sec. III we look close to the critical dimension d = 2 of the pure annihilation reactions to discover whether the branching processes are relevant (thus leading immediately to the active phase) or otherwise. For m even we find that they are in fact irrelevant close to two dimensions, leading to a transition into the active phase at zero branching rate, with, however, non-mean-field exponents which may be evaluated within an ǫ expansion. This calculation suggests, however, that this situation might be reversed in sufficiently low dimensions, and in Sec. IV we confirm this by an exact calculation in one dimension, using both field-theoretic and exact lattice methods. The subsequent non-trivial transition into the active phase is described within the truncated loop expansion in Sec. V, and the various subtleties connected with the unusual fixed point which emerges are investigated. Comments are given on the feasibility of extending this calculation to higher orders. In Sec. VI a generalization to N species is introduced and solved, in the hope of elucidating the case N = 1. However, it turns out that the behavior as soon as N ≥ 2 is quite different, leading to yet another universality class. In the next Section we turn to the case of m odd, and use diagram resummation and RG methods to argue that for d ≤ 2 the fluctuation effects drive the existence of a non-trivial inactive phase, counter to the predictions of mean-field theory. After a brief discussion in Sec. VIII of the corresponding set of problems when branching is added to the three-particle annihilation process 3A → ∅, we conclude with a summary of our results and possible directions for future theoretical and numerical work in this area.

II. MASTER EQUATION AND FIELD THEORY
In order to include fluctuation effects systematically in a mathematical description of the diffusion-limited reactions (1.1)-(1.3), we first write down the corresponding master equation and then derive a field-theoretic representation for it, following standard procedures [3][4][5]10,6,2]. The latter may then be treated with standard methods like perturbation theory, and the long-time and long-distance scaling behavior may be inferred via the application of the renormalization group.
Master equations for the temporal evolution of the probability distribution P (α; t) for a configuration α have the general form where R α→β denotes the transition rate for the process α → β. For random walkers on a lattice (with lattice constant b 0 ) subject to the reactions (1.1)-(1.3), any configuration is characterized by the integer site occupation numbers n i , i = 1, . . . , N , with i n i = N A (t), and the master equation may be decomposed into the balance equations for the diffusion, annihilation, branching, and pair production processes at each site i, (n e + 1) P (. . . , n i − 1, n e + 1, . . . ; t) − n i P (. . . , n i , n e , . . . ; t) , where {e} denotes the set of nearest-neighbor sites adjacent to i, and finally The precise form of the initial state is not important as we shall largely be considering stationary state properties. For convenience it is simplest to choose an uncorrelated Poisson distribution As the configurations are given entirely in terms of the occupation numbers n i , this calls for a representation in terms of second-quantized bosonic operators whose effect on the state is to raise or lower the site occupation number n i by one, respectively: Then, upon defining the time-dependent state vector the master equation (2.2)-(2.6) may be recast into an 'imaginary-time' Schrödinger equation with the normal-ordered 'hamiltonian' (2.14) Using the initial state (2.7), the formal solution to Eq. (2. Our aim is of course to compute time-dependent expectation values of observables A, which may be defined in terms of the configuration probability according to (2. 16) for this, we need the projection state For then, as a consequence of 0| e ai being a left eigenstate of a † i with eigenvalue 1, one has P | Q({a † i }, {a i }) = P | Q({1}, {a i }) for any normal-ordered polynomial Q of the ladder operators. Therefore, the expectation value (2.16) may be written as Notice also that probability conservation requires and therefore the hamiltonian generally has to vanish when the creation operators a † i are formally set to 1. We remark that in order to calculate 'inclusive' probabilities, such as the expectation value of the local density n i (t) = a † i (t)a i (t) , it is convenient to commute the factor e i ai(t) through the hamiltonian in Eq. (2.18), which is equivalent to a shift of all a † i → 1 + a † i , because then one just has to compute a vacuum expectation value of a normal-ordered product, for which Wick's theorem applies. Yet, this is not necessary for the evaluation of 'exclusive' quantities, as, e.g., the probability that site i is occupied by one particle, while all the other sites are empty, δ ni,1 j =i δ nj ,0 . Certainly this shift is not required for the renormalization of the model, and may in fact be even dangerous (as we shall see below).
In order to arrive at a field-theoretic representation, we now view the computation of (2.18) as a bosonic quantum many-particle problem with the hamiltonian (2.14), and again follow standard procedures [28] to write the expectation value as a coherent-state path integral , (2.20) with the effective action Finally, we perform the formal continuum limit, using the (assumed) inversion symmetry of the lattice, which results in the action According to Eq. (2.20), S has to be dimensionless, and using Eq. (1.13) we find the scaling dimensions of the fields to be and thus recover Eqs. (1.14) and (1.15). As mentioned in Sec. I, fluctuation effects are therefore expected to be important in physical dimensions in the cases k = 2 (d ≤ 2) and k = 3 (d ≤ 1). It is important to realize that in the above derivation of Eq. (2.23), the form of the nonlinear terms followed directly from the master equation, and no further assumptions had to be made about these fluctuation contributions. On the other hand, the mean-field rate equation may be recovered from this action by considering the 'classical' equations of motion for the fieldsψ and ψ, as given by the stationarity conditions Obviously, Eq. (2.25) has the stationary homogeneous solutionψ = 1, and inserting this into (2.26) yields which, upon identifying ψ(x, t) with a coarse-grained local density (which is only possible on this mean-field level), is the natural generalization of Eq. (1.4).
Notice that for even k and even m the hamiltonian part (in brackets) of the action (2.23), which does not explicitly depend on the final or initial state, is invariant under simultaneous 'parity' transformations of the fields, (2.28) Physically, this corresponds to the fact that for even k and m the particle number is locally conserved modulo 2 by both the annihilation and branching processes. However, if one expands the action (2.23) about the stationary solution according toψ (x, t) = 1 +ψ(x, t) , (2.29) corresponding to commuting the factor e i ai through the hamiltonian in Eq. (2.18), one arrives at the new action in which the symmetry (2.28) for even k and m is obscured. Even worse, if in accordance with usual naive power counting arguments only the three-point vertices are kept and all the higher nonlinearities are disregarded, this hidden symmetry is lost completely. (Reasoning along these lines led a previous investigation to the erroneous conclusion that a field-theoretic treatment would predict that the dynamic phase transition for BARW with k = 2 should be in the directed percolation universality class irrespective of the parity of m [11].) For a consistent renormalization group analysis, it is of course imperative to preserve all the symmetries of the problem. In the present case, this may be done by observing that the RG equations themselves (as opposed to the calculations of observables such as the density) should be independent of which basis is used, and it is therefore possible, and, indeed, necessary, to perform the computations in the representation of the model in which the symmetry for even k and m is manifest, namely the unshifted action (2.23).

III. BARW WITH TWO-PARTICLE ANNIHILATION NEAR TWO DIMENSIONS
A. RG eigenvalues at the annihilation fixed point and generation of new processes We now begin the investigation of fluctuation effects for BARW with two-particle annihilation (k = 2, λ 2 ≡ λ), based on the action (2.23). For 0 < t < t 0 , this field theory is characterized by a propagator Θ(t) exp −(Dq 2 + σ)t (in the time domain) or (−iω + Dq 2 + σ m ) −1 (in the frequency domain), and is graphically represented by a directed line, see Fig. 1(a), and by the vertices corresponding to the annihilation [ Fig. 1(b,c)], branching [ Fig. 1(d)], and pair production [ Fig. 1(e)] reactions, respectively. With these elements, a systematic perturbation expansion may be constructed applying the usual techniques of quantum field theory [29]. Taking into account fluctuation effects, we then define the renormalization constants Z λ , Z σm , and Z τ according to The dimensionless renormalized couplings ℓ and s m introduced here are to be determined from the vertex functions Γψψ ψψ (−q/2, −ω/2; −q/2, −ω/2; q/2, ω/2; q/2, ω/2) and Γψψψ ψ (−q/2, −ω/2; −q/2, −ω/2; −q/2, −ω/2; 3q/2, 3ω/2), respectively, at some specified symmetric normalization point setting a momentum scale κ, e.g., q 2 /4 = κ 2 , ω = 0, or q = 0, iω = 2Dκ 2 . In Eq. (3.1), C d ≡ Γ(2 − d/2)/2 d−1 π d/2 denotes a d-dependent geometric factor (C 1 = 1/2, C 2 = 1/2π). The renormalized pair production rate T may be obtained from the vertex function Γψψ in the time domain, by introducing a lower cutoff t i and taking the Laplace transform with respect to t i at iω = 2Dκ 2 . Similarly, the renormalized annihilation and branching rates ℓ and s m may also be computed from the vertex functions Γ ψψ and Γψ ψ , respectively, in the time domain with an upper cutoff t f . Near the upper critical dimension d c = 2, we infer the scaling behavior by studying first the ultraviolet (UV) divergences of the perturbation theory, to be absorbed in the renormalization constants Z λ , Z σm , and Z τ defined in Eqs. (3.1)- (3.3). From the solution of the RG (Callan-Symanzik) equations for the correlation functions via the method of characteristics, we derive running couplings given by the differential RG flow equations A scale-invariant asymptotic regime is then described by a zero ℓ * of the beta function, and the corresponding power laws are given by the anomalous dimensions (RG eigenvalues) at such a fixed point of the RG transformation, which is infrared-stable if β ′ ℓ (ℓ * ) > 0. The obvious first question to be addressed is whether the particle creation reactions with rates σ m and τ remain relevant perturbations to the pure annihilation process, i.e., whether the critical point remains at σ c = τ c = 0, once fluctuation effects are taken into account for d ≤ 2 dimensions; and if so, what the corresponding critical exponents are.
To one-loop order, the only UV-divergent Feynman diagrams at d c = 2 which contribute to the renormalization of the annihilation rate λ, the branching rate σ m , and the pair production rate τ are depicted in Fig. 2(a-c). These diagrams are identical with the lowest-order contributions at the annihilation fixed point, where σ m = τ = 0, provided the mass term ∝ σ m in the propagator is set to zero. Upon using the dimensional regularization and minimal subtraction scheme, one then finds the following renormalization constants in d = 2 − ǫ dimensions, The factor m(m + 1)/2 for Z σm simply originates in the number of possibilities to select the external legs to the right of the bubble in Fig. 2 , and (3.13), we infer the beta function In the pure annihilation model (σ m = τ = 0), this is actually an exact result, which may be obtained either by summing a geometric series of the bubble diagrams in Fig. 2(a), or by writing down the equivalent Bethe-Salpeter equation for the annihilation vertex [6], leading to and hence (3.16). For d > 2 the stable fixed point is just the 'Gaussian' fixed point ℓ * = 0, while for d ≤ 2 the non-trivial annihilation fixed point becomes stable and governs the asymptotic behavior. E.g., as there is no renormalization of the propagator, and hence neither field nor diffusion constant renormalization in the pure annihilation model, the solution of the RG equation for the average density may be written as Near the fixed point (3.18), ℓ(l) → ℓ * , we then match e 2l = t κ 2 , and arrive at the already cited exact result (1. 16).
Notice that in terms of the fixed point (3.18) the renormalization constant can be written as and hence ℓ = ℓ * corresponds to an infinite bare annihilation rate λ, c.f. Eq. (3.1) [6]. Moreover, Eqs. (3.8), (3.9), (3.14), and (3.15) yield the zeta functions and according to Eqs. (3.11) and (3.12) the RG eigenvalues for the branching and pair production processes at the pure annihilation fixed point (3.18) become for d ≤ 2 Thus, both these couplings remain indeed relevant near d c = 2. However, it is important to note that in addition to the original processes all the lower branching reactions with m − 2, m − 4, . . . offspring particles become generated via fluctuations involving combinations of branching and annihilation processes, see Fig. 2(d). In a consistent RG treatment all these additional reactions must be incorporated in the effective action from the start. In the case of odd m, the most important of these new processes is the single-particle annihilation A → ∅. Upon assigning to that the decay rate µ, decay : A → ∅ , rate µ , (3.25) which has the scaling dimension [compare Eq. (1.15)], the complete effective action reads The emergence of the spontaneous decay (3.25) alters the behavior of the system drastically, and in fact leads to a shift of the critical point to σ c > 0 for d ≤ 2, with the ensuing dynamic phase transition being in the directed-percolation universality class, and n(t) approaching its asymptotic value exponentially in both the active and inactive (absorbing) phases. We shall defer the discussion of the action (3.27) and the derivation of these results to Sec. VII. For even m, on the other hand, just the lower even-offspring branching processes are generated, and the effective action becomes However, Eq. (3.23) implies that the most relevant of these branching reactions is actually the one with smallest m, i.e., m = 2, with rate σ 2 ≡ σ and RG eigenvalue The branching process with m = 2 will therefore describe the entire universality class of BARW with even offspring number m, clearly distinct from the universality class of BARW with odd m as a consequence of the underlying 'parity' symmetry (2.28), i.e., the conservation of particle number modulo 2. In the following chapter, we shall derive the critical exponents for the BARW with m = 2 near two dimensions, and furthermore discuss the logarithmic corrections induced by the marginality of the annihilation vertex in d c = 2.
B. The case m = 2: Critical exponents near d = 2 and logarithmic corrections We have just argued that the generic universality class of the dynamic phase transition for BARW with even offspring will be characterized by the critical exponents for m = 2. To one-loop order [see Figs. 2(a-c)], i.e., to first order in ǫ = 2 − d, we found the RG eigenvalues (3.29) and (3.24) at the annihilation fixed point (3.18). Hence both σ and τ are relevant near d c = 2, and the critical point remains at its mean-field location σ c = τ c = 0. Upon utilizing the mean-field result (1.5), we write the general solution of the RG equation for the density as n(κ, ℓ, s, T, t) = κ d e −d l s(l) ℓ(l)n t κ 2 e −2l ℓ(l) , s(l)/ℓ(l) , T (l)/ℓ(l) ; (3.30) the pure annihilation result (3.19) is then recovered as s → 0 and T → 0, if we demand thatn(x, y, 0) → y −1 as y → 0 andñ(x) = lim y→0 [yn(x, y, 0)]. Therefore at the critical point the density decays as n(t) ∝ t −d/2 , i.e., using Eq. (1.8), For nonzero branching and pair production rate, but in the vicinity of the annihilation fixed point ℓ(l) → ℓ * , Eq. (3.30) reduces to n(κ, ℓ, s, T, t) ∝ s e (yσ−d) ln t κ 2 e −2l ℓ * , s e yσl , T e yτ l . (3.32) As lengths scale as e l , the correlation length diverges according to ξ ∝ s −1/yσ as s → 0 and ξ ∝ s −1/yτ as T → 0; therefore we identify [see Eqs. (1.10), (1.12)] Finally, the matching condition e 2l = tκ 2 leads to which, upon comparison with (3.33) and (1.11) immediately implies Furthermore, for T = 0 and t → ∞ the t dependence has to cancel, i.e.,n(ℓ * , x, 0) ∝ x −1+d/yσ as x → ∞. Hence with the definition (1.9) we conclude that Notice that the relations (3.31), (3.33), (3.35) and (3.36) are exact provided that σ c = τ c = 0. Only the RG eigenvalues y σ and y τ need to be computed perturbatively, e.g., in an expansion near d c = 2. Yet, assuming that the one-loop result (3.29) may be extrapolated with reasonable accuracy to some finite value of ǫ = 2 − d, one finds that the branching reaction actually appears to become irrelevant for ǫ > 2/3, i.e., below a supposedly new critical dimension Notice also that the exponents ν and β diverge as d ′ c is approached from above. This would suggest the possibility of an entire inactive phase to emerge for d < d ′ c and σ < σ c , i.e., the dynamic phase transition would be shifted to some positive critical value of the branching rate. For σ < σ c the effective branching rate would then scale to zero under renormalization, and the power laws of the pure annihilation theory would apply. In such a situation the above scaling relations and values of the critical exponents would have to be modified considerably.
In order to pursue this issue further, we now determine the RG eigenvalue y σ and y τ at the annihilation fixed point to second order in ǫ. For that purpose, we need the UV-divergent two-loop diagrams for the annihilation and branching vertices, as depicted in Fig. 3(a,b). Taking the external momenta or frequencies at a symmetry point and using either of them to fix the normalization scale, one may readily compute the Z factors for λ and σ, as defined in Eqs. (3.1) and (3.2); using dimensional regularization and minimal subtraction, one eventually finds Furthermore, by introducing a lower cutoff t i for the diagrams contributing to the renormalization of the pair production rate, see Fig. 3(c), one finds As to be expected, Eqs. (3.7) and (3.38) result in the same beta function (3.16) as to one-loop order, for the diagrams in Fig. 3(a) are precisely those of the pure annihilation theory, for which (3.16) holds to all orders in ℓ, and therefore the fixed point is again ℓ * = ǫ. With Eqs. (3.39) and (3.40), the zeta functions for the branching vertex (3.8) and for the pair production rate (3.9) become and thus We now employ Eq. (3.30) to infer the scaling behavior in the asymptotic limit l → ∞, where ℓ(l) ∝ l −1 , s(l) ∝ e 2l l −3 , and T (l) ∝ e 4l l −1 , and thus Consequently, upon identifying ξ with e l again, and by matching the second and third argument ofn to one, respectively, we find in addition, from the prefactor ofn we infer, using 2l ≈ ln(1/s),  Matching the argument ofñ to a constant then leads to which is consistent with the solution of the mean-field rate equationṅ = −λ(t) n 2 with a time-dependent annihilation rate λ(t) ∝ 1/ ln t, see Eq. (3.46) [6]. For large initial densities n 0 ≫ n s , this logarithmic correction to the meanfield power-law relaxation may be observable even for nonzero branching rate, in an intermediate time window Similarly, the non-trivial logarithmic corrections corresponding to Eq. (3.47) should appear in the opposite limit n 0 ≪ n s for times (2 σ) −1 ≪ t ≪ (2 λ n 0 ) −1 , inducing a subexponential increase (3.55) At very long times, however, any positive branching rate will eventually lead to a purely exponential approach to n s . This crossover to an asymptotically Gaussian theory may actually be described by simply keeping the mass term σ in the propagator for the one-loop diagrams depicted in Fig. 2(a-c), which then results in the following beta and zeta functions, (for a related discussion of the crossover from a critical point to a Gaussian theory, see Ref. [30]). The effective coupling emerging here is and its flowing counterpart g(l) will vanish as l → ∞ because of the divergence of s(l), as long as s remains a relevant coupling. In this situation, i.e., for d > d ′ c , the Gaussian fixed point g * = 0 describing the active phase is approached, and the above flow functions (3.56)-(3.58) assume their mean-field values, see Sec. V.
In summary, we have found that near two dimensions the BARW with even offspring, described by the universality class of m = 2, are characterized by a critical point at vanishing branching and pair production rate, with critical exponents α = d/2, z = 2, and only two non-trivial exponents ν and ν τ , which we have computed to second and first order in ǫ = 2 − d, respectively. We have also calculated the leading terms for the logarithmic corrections appearing directly in two dimensions, and identified the active phase with a Gaussian fixed point in terms of the effective coupling (3.59). However, there remains the possibility that in one dimension the branching process may actually be irrelevant near the annihilation fixed point, opening up an entire inactive phase for 0 ≤ σ < σ c .

IV. BARW WITH TWO-PARTICLE ANNIHILATION AND EVEN NUMBER OF OFFSPRING IN ONE DIMENSION
A. Computation of yσ m directly in d = 1 In Sec. III it was shown that the RG eigenvalue y σm which characterizes the relevance of the branching rate at σ = 0 may be computed within the ǫ = 2 − d expansion, and, although it is relevant in d = 2, y σm starts to decrease as ǫ increases, indicating that it may eventually change sign. In this section we compute its value exactly in d = 1, by various methods, and show that this is indeed the case.
In order to do this, we need to compute some dimensionless physical quantity, at some particular length or time scale characterized by the normalization scale κ, as a function of the bare parameters σ m and λ, and understand how these latter should change as κ is varied, keeping the physical quantities fixed. Since we are only interested in the O(s m ) terms in the RG equations (where s m denotes the dimensionless counterpart of σ m , see Sec. III), we need to compute these physical quantities only to first order in σ m , although the computation should be carried out, in principle, to all orders in the annihilation rate λ, since for d = 1 λ is not small at the annihilation fixed point. Rather than the renormalized vertex functions considered in Sec. III, we shall consider a physically more accessible quantity, namely the probability P m+1 (t) that, given there is just one particle at t = 0, there are exactly m + 1 particles found at time t. For sufficiently small t, we have P m+1 (t) ∼ σ m t. Thus, we may define a renormalized version of σ m in terms of the derivative dP m+1 (t)/dt, evaluated at some later time t = 1/Dκ 2 . This is clearly quite different from the definition used in Sec. III; however, for universal quantities such as y σm , it should yield identical results. In order to make this renormalized branching rate dimensionless, we then definẽ where the tilde emphasizes the difference from the definition used in Sec. III. P m+1 (t) is given by the sum of all diagrams in the unshifted theory given by Eq. (2.23), which, in time-ordered perturbation theory, have exactly one particle at t = 0 and m + 1 particles when cut at time t. To first order in σ m , there is only one branching vertex. Even so, the set of diagrams which may contribute to this may be very complicated. An important set is illustrated in Fig. 4. It is clear that, once the initial branching has taken place, the remaining portion of any such diagram corresponds to a solution of the quantum (m + 1)-body problem for non-relativistic bosons interacting with a short-range repulsive potential of strength λ. For general values of d and λ this is intractable. However, in one dimension, and in the limit λ → ∞, these bosons should behave like free fermions and therefore the problem be solvable.
Fortunately, it is just this limit which is appropriate to discuss the behavior at the annihilation fixed point ℓ = ℓ * , for this corresponds to a divergent bare coupling λ. For d < 2 there are no ultraviolet divergences in the continuum theory, so that, on dimensional grounds, in one dimension P m+1 (t) has the form (σ In this expression λ is supposed to be expressed in terms of the renormalized coupling by c.f. Eqs. (3.1) and (3.20). From these equations we see that, as ℓ → ℓ * , λ → ∞, and that this is the same as the limit t → ∞. This limit is more easily understood in the model defined on a lattice. For then two particles annihilate with probability one when they land on the same site. In this case, the branching process does not make sense unless the offspring are placed on different sites. Thus, we consider a slightly different version of the model in which the dimensionful parameter λ is replaced by the lattice spacing b 0 . Once the m particles are placed on neighboring sites, they execute independent random walks until two or more of them fall on the same site and they annihilate. Let S m+1 (t) be the probability that all m + 1 particles, initially placed on neighboring sites, survive until time t. Then we have the simple relation that so thats m = σ m S m+1 (t = 1/Dκ 2 ) + O(σ 2 m ). The asymptotic behavior of S m+1 is simple to compute. This problem has been studied in d = 1 in several contexts before [31]. Denoting the coordinates of the particles by (x 1 , x 2 , . . . , x m+1 ), we observe that this vector undergoes an isotropic random walk in the region x j+1 − x j > 0 of an (m + 1)-dimensional space, with absorbing boundaries along the hyperplanes x j+1 = x j , beginning at the point x 0 j = jb 0 . S m+1 (t) is just the survival probability of this walk. The Green's function for this problem is where A m+1 is the completely antisymmetrizing operator on the m + 1 coordinates x j , which ensures that the boundary condition is satisfied, and G (m+1) is the usual (m + 1)-dimensional lattice Green's function on Z m+1 . In the long-time limit we may use a continuum approximation to this, with care: The survival probability is then given by integrating the full Green's function over the region x j+1 > x j . The result must be a function of only the scaled variables (x 0 i − x 0 j )/(Dt) 1/2 , each of which is O(b 0 /(Dt) 1/2 ) ≪ 1. The antisymmetry requires that the lowest term in an expansion in these variables must have the form This gives the result for finite lattice spacing, in the limit λ → ∞. The corresponding result in the continuum limit for finite λ is obtained by the simple replacement b 0 → (2/λ). The precise coefficient (which is not, in any case, relevant to the present calculation) may be obtained by studying the case m = 2, which is solvable in both limits.
Putting these results together, we then find that, in d = 1, exactly.
The form of the result in Eq. (4.5) reflects the fermionic nature of the particles for infinite λ. In fact, Eq. (4.6) may be derived more simply, if more formally, using this. Once the branching event has occurred, the particles propagate as an (m + 1)-body fermionic state. If we express this in terms of anticommuting annihilation and creation operators c i and c † i on each lattice site, the branching term in the hamiltonian has the form [for m even; in the case of m odd, the product runs from −(m + 1)/2 to (m + 1)/2, omitting j = 0.] The continuum limit of this expression is found by expanding in powers of jb 0 . This will be different from the bosonic case because the anticommuting nature of the c † allows each derivative to appear only once. Thus the lowest order, most relevant, term has the form so that this term is multiplied by an effective expansion parameter b m(m+1)/2 0 σ m , and its dimension is modified in a way corresponding precisely to the second term in Eq. (4.6).

B. Results from an exactly solvable case
Finally, we show how this result may, in the case m = 2, be derived from an exactly solvable lattice version of the problem. This is the model discussed by Takayasu and Tretyakov [17], in which a site is chosen at random, and, if there is a particle there, it is either moved to the left or right (with probability p), or two offspring are placed on the neighboring sites (with probability 1 − p). If the particles are interpreted as Ising domain walls, this model is a discrete time version of a kinetic Ising model solved by Droz, Rácz, and Schmidt [21], with Glauber dynamics at effectively zero temperature, and Kawasaki dynamics at effectively infinite temperature. In fact this model is always in an inactive phase (except for p = 0), and does not exhibit the non-trivial transition to an active phase which is one of our main concerns in this paper. As shown in Ref. [22], it is necessary to include another parameter in the model, which enhances the branching rate relative to the diffusion process, in order to find such a transition. Nevertheless, for the purposes of extracting the universal eigenvalue y σ , at the pure annihilation fixed point, this deficiency is unimportant.
Takayasu and Tretyakov [17] consider the probability Q r (t) that a given interval of length r contains an even number of particles. Since they assume translationally invariant initial conditions, it does not matter which interval. In the Ising spin language, this is simply related to the correlation function s i s i+r . They show that this satisfies the linear system of equations 10) where N is the number of sites. Although Takayasu and Tretyakov did not analyze this equation in detail, it is not difficult to take a continuum limit and extract the scaling behavior of the parameter 1 − p, which is the analog of the branching rate σ ≡ σ 2 . Introduce Q 0 so that the second equation above has the same form as the first, extended down to r = 1. This means taking On rescaling t →t = t/N and r →r = b 0 r, the continuum limit 1/N → 0, b 0 → 0 now gives a simple diffusion equation for all real r > 0, together with boundary condition at r = 0. This corresponds then to a heat diffusion problem on the half-line r > 0, with a radiating boundary at r = 0. There is a length scale L b = (1 − p)b 0 implicit in the boundary condition. For r ≫ L b , the solution at late times will approach that with the Dirichlet condition Q(r = 0) = 0. However this will be modified for r ≪ L b . The fact that 1 − p ∼ σ scales with length in this way immediately implies that y σ = −1, consistent with Eq. (4.6). In fact, for this continuum model, this is true to all orders in σ. This is because the non-trivial fixed point, determining the transition at p = 0, is infinitely far away in this scheme. Within this model, it is possible to determine P 3 (t) and hence check the calculation given earlier in this section. If there is initially only one particle in the system, then Q r (t = 0) = r/N . If p = 1 (no branching), this is an exact steady state solution, corresponding to a uniform heat current. Let us find the solution in perturbation theory in (4.14) To first order in σ, there are only one or three particles in the system, so that the density is , and the O(1 − p) contribution to P ′ 3 (t) therefore scales as t −3/2 . Following through the RG procedure described earlier then leads to y σ2 = −1 as before.
Within this model, it is possible to investigate the nature of the long-time behavior to all orders in σ 2 ∼ 1 − p. If p < 1, the boundary temperature Q 0 (t) will increase to a point where the heat current crossing the boundary equals that coming from large r. Thus Q 0 (t) ∼ 2(1 − p)/N , so that the density of particles, Q 1 (t) ∼ [1 + 2(1 − p)/p]/N . There is therefore a finite average number of particles in the whole system, which diverges as p → 0. It is easy to check that this is the asymptotic solution for any finite odd initial number of particles. On the other hand, if this number is even, then initially Q r (0) = r(N − r)/N . In this case, the only steady state solution is Q r = 0, and it is straightforward to check that the particle density Q 1 (t) ∼ t −1/2 , as expected in the inactive phase.

A. RG flows to one-loop order
In Sec. III B we saw that near d c = 2 dimensions, the branching rate σ for BARW with two offspring particles remains a relevant perturbation. Thus the critical point is at σ c = 0, as in mean-field theory, and for any positive branching rate there exists only the active phase, which is governed by a Gaussian fixed point. On the other hand, the direct analysis of the one-dimensional model in Sec. IV unambiguously established the existence of an inactive phase described by the power laws of the pure annihilation model. Hence, in the inactive phase, σ must constitute an irrelevant operator in the RG sense. Clearly, this behavior cannot be constructed within an expansion about the upper critical dimension d c = 2. Yet, by employing a truncated loop expansion at fixed dimension we can establish a unifying RG framework for the above results. To one-loop order, we shall indeed find a new critical dimension d ′ c = 4/3, below which the inactive phase emerges. The complete set of tree and one-loop diagrams for the two-and four-point vertex functions is displayed in Fig. 5. As opposed to the analysis in Sec. III, we now retain the full dependence of both the vertices and the propagators on the branching rate σ. Upon applying an upper cutoff t f in time for the processes depicted in Figs. 5(a,b), taking the Laplace transform with respect to t f at the normalization point iω = 2Dκ 2 , and employing our previous definitions (3.1) and (3.2), one readily arrives at the renormalization constants Consequently, with Eqs. (3.7), (3.8) and ζ ℓ = β ℓ /ℓ, we find the zeta functions For s = 0, of course, these reduce to Eqs. (3.16) and (3.21) (for m = 2) as computed at the annihilation fixed point (3.18). Within this truncated one-loop expansion, we find for the RG eigenvalue of σ at this fixed point which indicates that indeed for d < d ′ c = 4/3 the branching processes become irrelevant. Note that in one dimension, this happens to coincide with the exact result (4.6) for m = 2. For d > d ′ c , however, the annihilation fixed point is unstable with respect to branching, and s(l) will diverge as l → ∞. In this case we expect the flow of ℓ and s to be described by their naive scaling dimensions. Indeed, with ℓ(l) ∼ e −(d−2)l and s(l) ∼ e 2l , we have ℓ(l)/s(l) 2−d/2 ∼ e −2l → 0, and the fluctuation contributions to Eqs. (5.3), (5.4) vanish. In Sec. III B, we had already anticipated these flow functions describing the crossover to the asymptotically Gaussian model that characterizes the active phase.
Before discussing the ensuing flow equations in more detail, we compute the renormalization of the pair production rate as well. Notice that in this calculation at fixed dimension, we cannot just restrict ourselves to those diagrams that become ultraviolet divergent at d c = 2. The three graphs in Fig. 5(c), computed with a lower cutoff t i , lead to the Z factor For σ = 0, this obviously reduces to Eq. (3.15); on the other hand, as σ → ∞, we may expand the term in square brackets with the result With the definition (3.9) this leads to which is identical to Eq. (3.58). Note that Eq. (5.8), which has been obtained to leading order in an expansion for large σ/Dκ 2 , correctly incorporates both the inactive (s = 0) and active (s → ∞) phases, as do Eqs. (5.3) and (5.4). It is instructive to realize that one arrives at the same results for ζ ℓ and ζ s by considering the four-point vertex functions Γψψ ψψ (−q/2, −ω/2; −q/2, −ω/2; q/2, ω/2; q/2, ω/2) and Γψψψ ψ (−q/2, −ω/2; −q/2, −ω/2; −q/2, −ω/2; 3q/2, 3ω/2), respectively, at either the normalization point q 2 /4 = κ 2 , ω = 0, or at q = 0, iω = 2Dκ 2 . This comes about because the expressions from the triangular loops in Figs. 5(d,e) may be written in the form of the last term of Eq. (5.6), and hence do neither contribute at the annihilation nor at the Gaussian active fixed point. Finally, the loop in the fourth diagram in Fig. 5(d) carries no external momentum or frequency, and thus has no influence on the zeta function. As in the one-loop ǫ expansion near d c = 2, there is neither field renormalization nor renormalization of the diffusion constant at this level, because the loop diagram for Γψ ψ (q, ω) (with no upper cutoff in the time domain), see Fig. 5(b), carries no wave vector or frequency dependence.
Returning to Eqs. (5.3), (5.4), and (5.8), it is clear that the relevant effective coupling is indeed g = ℓ/(1 + s) 2−d/2 , as defined in Eq. (3.59). In the prospective inactive phase, with s → 0, we have g → ℓ, the coupling of the pure annihilation model. On the other hand, if s diverges, we may construct the beta function for g as follows: Thus, as explained above, despite both ℓ(l) and s(l) diverging according to their naive scaling dimensions, there is a stable Gaussian fixed point g * = 0 describing the active phase. Yet, Eq. (5.9) yields the additional non-trivial fixed point assumes the value −2 at g = g * , this new fixed point turns out to be infrared-unstable. As the fixed-point values for the new effective coupling are limited by those of the annihilation rate, which in turn has as upper bound ℓ * = 2 − d (for according to Eq. The ensuing phase diagram as function of space dimension is displayed in Fig. 7. For d > d ′ c and any nonzero branching rate there exists only the active phase. For d > d c = 2, the transition at σ c = 0 is described by the meanfield exponents, see Sec. I, while for d ′ c ≤ d ≤ 2 one finds the non-trivial critical exponents (logarithmic corrections at d c = 2) evaluated in Sec. III B. Only for d < d ′ c , when g * < ℓ * = 2 − d and the annihilation fixed point becomes stable for sufficiently small values of s, the inactive phase which is governed by the power laws of the pure annihilation model emerges. The critical behavior at the dynamic phase transition with σ c > 0 is characterized by exponents belonging to a novel universality class. In the following subsection, we shall analyze the scaling behavior at this transition, and compute the critical exponents to one-loop order. Although the actual values we find for these exponents turn out to be rather poor estimates, as a consequence of the absence of any small expansion parameter, we believe that the qualitative features of the phase diagram, and specifically the mechanism for how the inactive phase becomes possible, are correctly encoded in the truncated loop expansion.

B. Scaling analysis and critical exponents
We now want to explore the critical behavior in the vicinity of the transition described by the unstable fixed point g = g * , Eq. (5.10). We thus introduce 12) and note that because of g * dε(l)/dl = β g (g(l)) = −g * ε(l) dβ g (g)/dg g=g * , we identify the corresponding zeta function as Collecting the results (5.3), (5.4), (5.8), and (5.11), we then find the following RG eigenvalues at the unstable fixed point g * , y ε = −ζ * ε = 2 , (5.14) We remark that y λ and y σ are not independent, but fixed by the condition that β g (g * > 0) = 0, and hence y λ = (2 − d/2)y σ . Notice also that to one-loop order there are no anomalous dimensions for the diffusion constant or the fields themselves.
In identifying these RG eigenvalues with the critical exponents, however, we have to be careful as now the critical point is shifted away from σ c = 0 to g c = g * , and we therefore cannot apply Eqs. (3.33) and (3.36). In general, the critical exponents will rather depend on y ε , as well as on y λ and y σ . Only the exponent ν τ is not affected by these modifications, and therefore By itself, this is not sufficient to determine the exponents, since s(l) and ℓ(l) do not themselves flow to fixed points and the dependence of n on these, even at early times, is not known a priori. In the active phase, however, we know that they do flow to a region in which mean-field theory is valid. This predicts that the right hand side of (5.19) has the form s(l)/ℓ(l) times a function of the combination t(κ 2 /D)e −2l s(l). In the vicinity of g * , then, (5. Notice that lengths now scale as e (1−yσ/2)l . The mean-field value for the dynamic exponent z naturally follows from the absence of diffusion constant renormalization. We emphasize that as a consequence of the appearance of the dangerous irrelevant variable 1/s, we cannot unambiguously extract the exponents defining the power laws at the critical point itself, such as α. In addition, this also precludes us to provide a sound foundation based on the renormalization group for scaling relations like β = zνα [11,19,22] for the dynamic phase transition for d < d ′ c . Such scaling relations would follow only if the mean-field dependence s(l)/ℓ(l) were also to be valid as both s(l) and ℓ(l) grow without bound along the critical line g = g * . If the fixed point occurred at finite values of these quantities, it would be reasonable to assume that the right hand side of Eq. (5.19) does have this mean-field form, as it is to be evaluated at some finite rescaled time. It would be interesting to look carefully for possible violations of these scaling laws within simulations.
Upon inserting Eqs. (5.14)-(5.17), we find the following results as function of dimension, These numerical values, compared to the actual simulation results [11,17,22,26] are generally rather poor, with the remarkable exception of ν τ . Keeping in mind that the truncated loop expansion, with no small parameter at hand, constitutes an uncontrolled approximation scheme, the unsatisfying accuracy may not be too surprising.

C. Some remarks regarding higher orders in the truncated loop expansion
It is by no means clear that including higher loop orders in the above analysis would yield considerably better results for the critical exponents, given that there is no a priori small expansion parameter present at all. Yet some general remarks about a possible extension, say, to two-loop order, are in place here.
To two-loop order, one would obviously expect both field renormalization and diffusion constant renormalization to appear, e.g., from the frequency and wave vector dependence of the loop in the fourth diagram of Fig. 3(b), or more generally, from the set of graphs depicted in Fig. 4. Defining this will imply two additional RG eigenvalues y ψ = −ζ * ψ and y D = −ζ * D . As may be inferred from Eqs. (5.19) and (5.20), their appearance modifies the previous scaling relations (5.21)-(5.23) to and would thus allows for a non-trivial value of the dynamic exponent as well. It is, however, far from obvious that an extension of the above one-loop analysis to higher orders in the perturbation expansion is feasible. Certainly, at the annihilation fixed point describing the inactive phase, a loop expansion is welldefined to any order. Also, the Gaussian fixed point is trivially described by the naive scaling dimensions of the model parameters. Yet the unstable fixed point governing the dynamic phase transition requires a subtle balancing of the divergences of the couplings s(l) and ℓ(l), as encoded in the appearance of the new effective coupling g(l). To one-loop order, the expansion parameter in the bare theory is (λ/D)/(κ 2 + σ/D) 1−d/2 , and hence the effective coupling in the flow function becomes g = ℓ/(1 + s) 2−d/2 . The two-loop diagrams, then, are proportional to (λ/D) 2 /(κ 2 + σ/D) 2−d , and thus their contributions to the zeta functions proportional to ℓ 2 /(1 + s) 3−d = g 2 (1 + s) -which diverges at the transition. It appears, therefore, that in an expansion in 1/s, the leading terms would have to cancel in order to render the loop expansion well-defined. However, although it may be shown for ζ ℓ that this is fact true for those two-loop diagrams that in the frequency domain are product of the one-loop terms, this appears not to be the case for the nested loop diagrams such as the fourth graph in Fig. 3(b), and perhaps the summation of an entire series of such diagrams, see Fig. 4, might be indispensable. Despite some effort, we therefore did not succeed in extending our one-loop analysis further, or even to demonstrate that there exists a meaningful truncated two-loop theory at all.

VI. GENERALIZATION TO N SPECIES OF PARTICLES
We now generalize the BARW with k = m = 2 to N species of particles α = 1, . . . , N , with only identical particles annihilating each other, and two different types of branching processes,

2)
3) pair production : The motivation behind such a generalization is that we found in Sec. V that in the interesting physical case of N = 1 in d = 1 the loop expansion was uncontrolled and does not lead to reliable values for the exponents, even though it correctly predicts the existence of a transition. As in the case of N -component magnets, then, one might hope that there exists an N -species generalization of the model with m = 2 introduced in Sec. II, which becomes exactly solvable in the limit N → ∞, yet which still retains the essential physics. Unfortunately this does not seem to be so. There appears to be no simple way of generalizing the action (3.28) (with m = 2) so that it possesses an O(N ) symmetry as in the magnetic case. The highest symmetry that can be embedded while at the same time retaining the conservation modulo 2 of each species is S N , the permutation symmetry of the N species. This is clearly exhibited by the allowed reactions above. Note, however, that the permutation symmetry allows different possible rates for the reactions in (6.2) and (6.3). The symmetries of this model are most clearly seen in the form of hamiltonian density in the second-quantized formalism: When σ ′ = (N − 1)σ, the branching terms do in fact exhibit an O(N ) symmetry which generalizes the Z 2 symmetry of the original model. However, the annihilation terms inevitably violate it, and any attempt to rectify this falls foul of the requirement of probability conservation. This has the consequence that the symmetry of the branching terms is not preserved under renormalization. Physically, this is because after the first process (6.2), any pair of the three outgoing particles may annihilate, so the chances of observing three particles after some amount of time are reduced as compared with process (6.3), where only one pair may annihilate. The two-loop calculation of Sec. III B is readily generalized to N particle 'flavors' by modifying the appropriate combinatorial factors in the evaluation of the Feynman diagrams depicted in Fig. 3. Near two dimensions, this yields the Z factors (3.38), (3.39) as before, while the graphs in Fig. 3(b) for the new branching process (6.3) give Consequently, with the definitions used in Sec. III, the beta function (3.16) and fixed point (3.18), as well as Eq. (3.41) remain unaltered, while the additional zeta function reads and therefore the RG eigenvalue of the branching rate σ ′ at the annihilation fixed point is Upon comparing this result with (3.43), we see that the branching with identical offspring (6.2) is irrelevant as compared to the production of different particle 'flavors' (6.3): for l → ∞, one has s(l)/s ′ (l) ∼ exp{(y σ − y σ ′ ) l} = exp{[−2ǫ + 3 ln(4/3)ǫ 2 ] l} → 0 . (6.9) Therefore the process (6.3) is going to dominate the long-time behavior, and we may effectively set σ = 0, but retain σ ′ > 0, for any N > 1. The remaining Feynman diagrams are then identical to those that would be obtained in the limit N → ∞, and one therefore expects that actually all reactions with N > 1 are asymptotically described by the model with infinitely many particle species and sole branching process (6.3). The above model in the limit σ = 0, or equivalently, for N → ∞, has a considerably simpler structure as compared to the N = 1 case, because 'nested' diagrams as, e.g., the last graph in Fig. 3(b) do not appear any more. Consequently, all fluctuation contributions to the propagator may be absorbed into a renormalization of the 'mass' σ. Therefore one is left with the task of summing an infinite series of branching 'bubble' diagrams ('cactus' graphs). The most convenient way of achieving this is to write down self-consistent Bethe-Salpeter type equations for the annihilation, branching, and pair production rates. Graphically, these are represented in Fig. 8(a-c), respectively. The corresponding analytic expressions read 11) which are readily solved by If we now evaluate the integral at the normalization point σ ′ R /D = κ 2 , we arrive at the exact result which immediately implies that Eqs. (3.16), (3.22), and (6.7) hold to all orders in ℓ, Thus in this model, σ ′ c = τ c ≡ 0, and the only non-trivial RG eigenvalues are which determine the divergence of the correlation length as either σ ′ → 0 or τ → 0.
We have thereby demonstrated that the above N -species generalization of the BARW with m = 2 leads to a new universality class, which for all N > 1 is asymptotically characterized by the model with N → ∞. In this model, the critical point remains at vanishing branching and pair production rate, and for d < 2 its critical behavior is governed by the exponents ν = 1/d , ν τ = 1/2d , z = 2 , α = d/2 , β = 1 . (6.17) In d c = 2 dimensions, one finds a logarithmic correction for the density decay at the critical point, just as in Eq. (3.54).
The N → ∞ limit of the model defined by the processes (6.1)-(6.4) may in fact be solved directly, without going through the renormalization group machinery. This is based on the simple observation that in the reaction A α → A α + A β + A β , the two products A β are far more likely to annihilate against each other, since they have the same flavor index, than against the remaining A α , or indeed any other particle in the system, which is unlikely to have a matching flavor index. In the limit N → ∞, this becomes overwhelmingly the case. Now the average particle density can change by branching processes, which will always be proportional to the existing local density, and annihilation processes. These latter are either between particles which are siblings, that is they are products of the same branching event in the past, and which will not, in this limit, annihilate against any other particle; or they are between particles which are from independent branching events. In this case they are unlikely to have the same flavor index, so that this term is suppressed by a factor 1/N ; however, it is proportional to the square of the mean density, which is itself O(N ), and therefore it also gains a factor of N relative to the linear terms. If the particles come from independent branching events, we may also, to leading order in 1/N , neglect the correlation between these events.
One may therefore write down an integral equation for the time evolution of the mean density n(t): where L(t − t ′ )dt is the probability that two particles, created at time t ′ , annihilate in the interval (t, t + dt). Note that L(t) = −S ′ 2 (t), where, as in Sec. IV, S 2 (t − t ′ ) is the survival probability at time t for two isolated particles produced at time t ′ . The factor of 1 − L(t − t ′ ) in the last term ensures that the particles have not annihilated at some earlier time.
Eq. (6.18) may be solved simply by Laplace transform. Definingñ(p) = ∞ 0 n(t)e −pt dt, and similarly for L(t), one finds pñ(p) − n(0) = 1 −L(p) 2σ ′ñ (p) − λ n 2 (p) , (6.19) whereL(p) is given by a simple sum of bubble diagrams as L(p) = 2λI(p) 1 + 2λI(p) , (6.20) and I(p) = k 1/(p + 2Dk 2 ). Although (6.18) cannot be solved in closed form, the various critical exponents may easily be deduced. If we begin from an initial state when n(0) is small, we may initially neglect the nonlinear term to find that, whenever σ ′ > 0,ñ(p) has a pole for positive real p at the solution p = p 0 (σ ′ ) of p [1 + 2λI(p)] = 2σ ′ . (6.21) The density therefore increases exponentially at a rate e p0t . When d > 2, I(p) is analytic at p = 0 and we see that p 0 ∝ σ ′ . When d < 2, p 0 ∝ σ ′ zν , where zν = 2/d, and for d = 2, p 0 ∼ σ ′ / ln(1/σ ′ ), all in accord with the renormalization group analysis presented above. For later times, this exponential growth is capped by the nonlinear term and the system reaches a steady state. Since the factor 1 −L(p) is common to both terms on the right hand side of (6.19), the steady state density behaves simply as σ β , where β = 1 irrespective of the value of d.

VII. BARW WITH ODD NUMBER OF OFFSPRING
A. Mapping to Reggeon field theory for directed percolation In this section, we consider the case of m odd, which turns out to be quite different from the case of even m. For odd m, there is no 'parity' symmetry and one might expect the transition from an inactive state, provided it occurs at some non-trivial value of the branching rate σ m , to be in the directed percolation (DP) universality class, as a consequence of the generated particle decay processes A → ∅. This in fact is what we find in low dimensions, but there are subtle issues connected with the fact that these decay processes are themselves fluctuation-induced, and thus proportional to the branching rate σ m . Therefore the very existence of a non-trivial transition must depend on fluctuation effects, which are expected to be important only when d ≤ 2. As pointed out in Sec. I, the mean-field equation always predicts an active state for all σ m > 0, and is expected to be valid, at least in the pure annihilation problem, for d > 2. Indeed we find that only for dimensions d ≤ 2, fluctuations are strong enough to produce a nonzero critical branching rate σ c > 0, and consequently a DP phase transition separating the active from the inactive phase. This is despite the fact that the upper critical dimension for the DP transition itself is d c = 4.
For these reasons we investigate this problem in the vicinity of two dimensions. As was shown in Sec. III, from an initial model with a given value of m > 1, under renormalization the other processes with numbers of offspring m − 2, m − 4, . . . , 1, −1 are generated, compare Fig. 2(d). Therefore we shall begin by considering the simplest case m = 1. We start by writing the interaction terms appearing under the integral in the effective action Eq. (3.27): The last term corresponds to the process A → ∅, which, although not present in the original model, is generated by the combined processes A → 2A, 2A → ∅. It is also important to record the form of the interaction in terms of the shifted fieldψ =ψ − 1. This is which, apart from the last term, has the form of the effective interaction for the (Reggeon) field theory of directed percolation [13]. With the simple rescaling φ = 2λ/σ 1 ,φ = σ 1 /2λ (7.3) one arrives at with the new effective coupling for the three-point vertices Hence u becomes marginal in d = 4, the upper critical dimension for DP, and indeed one would argue that the last term in Eq. (7.4) should then be irrelevant close to the DP transition by power counting. Thus, provided that σ 1c > 0, we have shown that the ensuing dynamic phase transition from the inactive to the active state is characterized by the critical exponents of the DP universality class. Notice that both the actice and the inactive phase are characterized by exponential long-time behavior, as opposed to the power-law inactive phase for BARW with even m. At first sight counterintuitively, the branching process accelerates the particle density decay in the inactive phase, due to the fluctuation-induced generation of spontaneous particle decay processes. Note, however, that in the bare theory the 'mass' term (the coefficient ofψψ) is negative, since µ = 0 in that case. This would seem to imply that the model is immediately in its active phase as soon as σ 1 > 0, and hence the critical point at σ 1c = 0, which would render the rescaling (7.3) obsolete. However, this argument neglects renormalization effects. It may well be that the renormalized quantity µ R − σ 1R is positive. As we shall argue, this is in fact the case for d ≤ 2 and sufficiently small σ 1 .

B. Summation of the leading singularities near two dimensions
We are free to study the renormalization effects in either version of the theory Eq. (7.1) or Eq. (7.2). In fact, this is much simpler in the former, because it will turn out that the effective expansion parameter is λ, and also because the mass term in that case is µ + σ 1 , which is always positive, allowing for a perturbative expansion with a meaningful propagator. Even so, it is impossible to compute the renormalization to all orders. Instead, close to two dimensions we adopt the strategy of retaining only the most singular diagrams as ǫ = 2 − d → 0, at each order in λ. It is not difficult to see that these are given by the iterated bubble diagrams shown in Fig. 9(a,b). Notice that all the propagators are renormalized by a series of processes as depicted in Fig. 9(c), but all these side branches may be resummed with the effect of replacing the mass in the internal propagator by its renormalized value. In the rest of this section, we set the diffusion constant D = 1, since it is not renormalized at the order we consider. The ensuing geometric sums then give the following simple self-consistent equations, to leading order in ǫ, Note that µ R + σ 1R = µ + σ 1 in this approximation. With µ = 0 originally, this amounts precisely to the summation of branched diagrams [see Fig. 9(c)] mentioned above. However, we are interested in the renormalized value of the DP mass (c.f. (7.4) where we have used the above result, with µ = 0. From Eq. (7.8) we see that, for sufficiently small σ 1 , the DP mass is in fact positive, indicating that the system is in fact in the inactive phase. The transition to the active phase does not happen until σ 1 = σ 1c , where, in this approximation, Although we would not expect this result to be quantitatively correct for ǫ = 1, the power of λ is exact and is dictated by dimensional analysis. Thus we expect a DP transition at a non-trivial value of σ 1 for all d < 2. This transition is driven by the fluctuations in the annihilation process.
Since the annihilation rate is (marginally) irrelevant exactly in two dimensions, it is an interesting question whether it is able to drive a DP transition in this case. The single bubbles are then logarithmically divergent, and Eq. (7.8) is replaced by where Λ is a wave number cutoff (of the order of the inverse lattice spacing b −1 0 ). We see that indeed there is a non-trivial transition, at σ 1c ∼ Λ 2 e −4π/λ . (7.11) This is one of the important results of this paper. It shows how essential it is to take into account fluctuation effects in studying this and similar processes. It appears to be in accord with recent detailed simulations [27], even though initial work [17] seemed to indicate that σ c might be zero. Finally, for d > 2, the same approximation (which is, however, less justified), gives where In this case, as long as λΛ d−2 ≪ 1, the numerator will be negative even at σ 1 = 0, and the system will be in the active phase for all σ 1 > 0. There does appear to be the possibility of a non-trivial transition at larger values of the annihilation rate, when λΛ d−2 ≫ 1, but it should be noted that our approximation of summing the leading bubble diagrams breaks down long before then, and that it may be argued that infinite annihilation rate on the lattice already corresponds to λΛ d−2 ∼ 1 in the continuum theory. There is no evidence for a non-trivial transition in d = 3 from the simulations [17].

C. RG approach
The same problem, with m = 1, may be attacked using the RG loop expansion used in Sec. III. Once again, we work in the unshifted theory given by Eq. (7.1). The one-loop diagrams contributing to the renormalization of λ, σ 1 , and µ are shown in Fig. 2(a,b,d), and in fact only the combinatorial factors differ from the case m = 2 studied above. Defining dimensionless renormalized couplings as in Eqs. (3.1), (3.2), and in addition M = Z µ µκ −2 /D , (7.14) we find the flow equations Note that the effective dimensionless mass appearing in the denominators is now M + s 1 . Adding the last two equations, we see that M (l) + s 1 (l) = σ 1 (0) e 2l (equivalent to our earlier result that the dimensionful mass does not renormalize). Although the remaining equations cannot be solved by quadrature, we may approximate their solution by setting 1 + M + s 1 ≈ 1 for l < l 0 ∼ 1 2 ln(1/σ 1 ). Beyond that point, the couplings simply evolve according to their simple scaling behavior at the Gaussian fixed point. The equations may now be integrated, with the result for the dimensionless renormalized DP mass We may now distinguish various cases. For d < 2 and sufficiently small σ 1 (large l 0 ), the exponent in Eq. (7.18) behaves like ℓ * l 0 ≪ 1, so that the DP mass is positive, consistent with the system being in the inactive phase. For large σ 1 , however, the exponent is small and the DP mass is negative. The critical point, in this approximation, occurs when For small initial annihilation rate λ we have ℓ(l ′ ) ∼ C d λe ǫl ′ . This gives σ 1c ∼ e −2l0 ∼ (C d λ/ǫ ln 2) 2/ǫ , in agreement with Eq. (7.9) apart from the factor of ln 2, which may be ascribed to the crudeness of the truncation procedure described above. Thus, as expected, the one-loop RG equations correctly sum the most singular diagrams in 1/ǫ in the weak coupling limit. However, from Eq. (7.19) we may also extract the behavior at strong coupling, because then ℓ(l ′ ) ∼ ℓ * = ǫ. In that case we find σ 1c ∼ e −2 ln 2/ǫ , (7.20) where, once again, the precise numerical factor is not to be taken too seriously.
In exactly two dimensions, on the other hand, l0 0 ℓ(l ′ ) dl ′ ∼ ln(1 + C d λl 0 ), c.f. Eq. (3.46), so that σ 1c ∼ e −4π/λ , in precise agreement with the summed bubbles calculation, Eq. (7.11). The only corrections to this formula are seen to come from higher order loops, but, since ℓ is now marginally irrelevant, it should asymptotically become exact as λ → 0. Finally, we may use the criterion Eq. (7.19) to explore the case d > 2. In that case we may write Since now ℓ(l 0 ) → 0 as σ 1 → 0, we see that there is a possibility of finding the inactive phase for sufficiently small, but nonzero σ 1 only if ln[(|ǫ| + C d λ)/|ǫ|] > ln 2, i.e., C d λ > |ǫ|. Once again, of course, one must question the validity of the loop expansion in this region. However, the result does indicate that the inactive phase does not exist for nonzero σ 1 at weak coupling.

D. Larger odd values of m
The advantage of the RG loop expansion approach outlined above for the case m = 1 is that it may readily be extended to larger values of m, where the identification and resummation of the most singular diagrams is somewhat more difficult. The only case that we shall treat explicitly here is m = 3. In that case we have a branching rate σ 3 , corresponding to a dimensionless renormalized coupling s 3 , but, as before, couplings σ 1 and µ are generated. (Note that σ 2 is also generated, but only at higher orders, and in any case this is irrelevant near d = 2.) Once again, the one-loop diagrams are simple, see Fig. 2(a,b,d) and differ from those for m = 1 only by combinatorial factors. The RG flow equations are Notice that the mass term is now M + s 1 + s 3 . As before, we see that this combination obeys a simple equation, and that it grows exponentially like e 2l . Thus, once again, we may make the same approximation and ignore it in the denominators up to l 0 ∼ 1 2 ln(1/σ 3 ). The remaining equations may now be solved by quadrature. We give the results only for the case d = 2. If we shift the fields we now find that the DP mass is ∆ ≡ µ − σ 1 − 3σ 3 . The dimensionless renormalized version of this quantity at scale l 0 then turns out to be M (l 0 ) − s 1 (l 0 ) − 3s 3 (l 0 ) = σ 3 e 2l0 1 − 12/5 1 + λl 0 /2π − 8/5 (1 + λl 0 /2π) 6 . (7.26) Once again, this is positive for large l 0 (small σ 3 ), indicating the existence of an inactive phase, while it is clearly negative for sufficiently small large σ 3 . The zero must be found numerically, and occurs at λl 0 /2π ≈ 1.42, corresponding to σ 3c ∼ e −5.68π/λ . (7.27) This is expected to be asymptotically exact (apart from a possible prefactor) as λ → 0. Note that, compared with the case m = 1, the transition occurs at a much smaller value of the branching rate. This is because the generation of the process A → ∅ now happens only at higher order in the annihilation rate λ, which is asymptotically small. This pattern of the decrease of σ mc with increasing m continues, and is consistent with what is observed in simulations [17]. These calculations for larger m may also be extended to the cases d = 2, with similar results as for m = 1.

VIII. BARW WITH THREE-PARTICLE ANNIHILATION
We finally briefly discuss the case of BARW with k = 3, i.e., three-particle annihilation processes instead of the two-particle annihilation A + A → ∅ which has concerned us in the bulk of this paper. In physical dimensions d ≥ 1, BARW with any higher value of k will either be satisfactorily described by the mean-field rate equation (1.4), according to the fact that the annihilation vertex λ k has upper critical dimension d c = 2/(k − 1), c.f. Eq. (1.15), or simply generate the processes with lower k, see below.
We may of course largely follow the arguments given in the preceding sections. First, via a combination of branching and annihilation the branching processes with m offspring particles (rate σ m ) now generate all the lower branching processes ∝ σ m−3 , σ m−6 , etc., see Fig. 10(a). Clearly, the reactions with lowest offspring number will be the most relevant ones, and therefore we shall essentially have to discuss only the cases m = 1, m = 2, and m = 3. For the latter, and equally for any m = 3l, l = 1, 2, . . ., the branching and annihilation reactions conserve the particle number locally modulo 3, and this symmetry precludes a fluctuation-induced generation of spontaneous particle decay processes, which might lead to a phase transition of the DP universality class. In analogy with the situation for even m in the case k = 2, there might potentially arise a new non-trivial universality class below some novel critical dimension d ′ c < 1; yet in one dimension one merely expects logarithmic corrections to the mean-field behavior. This leaves us with the cases m = 2 and m = 1 (or, more generally, m = 2 mod 3 and m = 1 mod 3). Here, the above special symmetry is not present, and via the processes depicted in Fig. 10(b) and (c), respectively, the singleparticle decay A → ∅ is generated. For m = 2, this already happens to lowest order in the branching and annihilation vertices, O(σ 2 λ 3 ), see Fig. 10(b). In addition, the processes A → A + A and A + A → ∅ also become generated, but require diagrams of O(σ 2 2 λ 3 ) and O(σ 2 2 λ 2 3 ), respectively. This leads to all the terms required for Reggeon field theory to represent the effective action for the ensuing dynamic phase transition, provided σ 2c > 0 (compare Sec. VII). Yet as already the lowest-order diagram in Fig. 10(b) is (logarithmically) divergent in one dimension, following the arguments given for the odd-offspring BARW with k = 2 we expect such a non-trivial transition to occur in d ≤ 1, with the exponents of the DP universality class. In fact, as the BARW with k = 2 and m = 1 becomes generated, albeit to some high order in the original couplings, the DP transition may actually persist for any d ≤ 2. For the case m = 1, the lowest-order diagram, O(σ 2 1 λ 3 ), diverges already for d ≤ 2 dimensions, and as again the process A + A → ∅ is generated via diagrams of order O(σ 4 1 λ 2 3 ), this maps onto BARW with k = 2 and odd m and hence one has to expect a non-trivial phase transition of the DP universality class for d ≤ 2 dimensions.

IX. SUMMARY AND CONCLUSIONS
We have studied branching and annihilating random walks defined by the combined processes 2A → ∅ (k = 2), A → (m + 1)A with both even and odd offspring number m in low dimensions, where fluctuation effects and particle anticorrelations are crucial. We have employed a field-theoretic representation of the corresponding diffusion-limited reaction problem, derived from the classical master equation. Apart from the continuum limit, the derivation of the action requires no further approximations (Sec. II).
While the mean-field rate equation predicts the existence of an active phase only for any nonzero value σ m of the branching rate (Sec. I), in low dimensions fluctuations lead to the emergence of a non-trivial inactive phase as well, and a dynamic phase transition at a critical value σ c (for fixed annihilation rate λ and diffusion constant D). As a consequence of the local parity conservation in the case of even m, which is reflected in a discrete symmetry in the action, both the inactive phase and the critical properties at the nonequilibrium phase transition differ drastically for even and odd m, respectively, as was already known from numerical simulations [11,17,19,22,26,27]. Using diagrammatic summations and renormalization group methods, we believe we have provided at least a qualitatively satisfying understanding of these different universality classes [9].
For odd m and d ≤ 2, including the borderline dimension d c = 2, fluctuations generate the spontaneous decay A → ∅, as well as branching processes with all odd offspring numbers smaller than m. This single-particle decay provides a very effective mechanism counteracting the particle production through branching, and induces a nontrivial dynamic phase transition characterized by the critical exponents of the directed percolation universality class, as has been established in simulations [10,17,26,27]. In the emerging inactive phase, the spontaneous decay processes dominate, and thus the long-time decay is exponential, i.e. faster than in the pure annihilation model without the branching processes. For dimensions d > 2, we find σ c = 0, and the system is essentially described by mean-field theory. Thus the values of the critical exponents display a discontinuous jump at d c = 2 (for the upper critical dimension of directed percolation is d c = 4). Precisely in two dimensions, we still predict a non-trivial transition with DP exponents, and first indications for this have recently been found [27]. Indeed, even the tendency that σ c should decrease with increasing m is at least qualitatively explained by our results (Sec. VII).
While the case of odd m is thus well understood now, there remain a number of open questions for even offspring numbers. Similarly to the situation for odd m, the even m universality class is determined by the lowest possible value m = 2. Our analysis shows that near two dimensions, the branching rate always remains a relevant quantity (i.e., an active phase governed by exponential correlations), and thus necessarily σ c = 0. In Sec. III, we have evaluated the critical exponents to second order in ǫ = 2 − d, and have discussed the logarithmic corrections arising at the critical dimension d c = 2. On the other hand, a direct and exact computation in one dimension (Sec. IV) establishes that the branching rate is actually irrelevant in the limit λ → ∞, which implies the existence of an inactive phase only in this special situation, governed by the power laws of the pure annihilation model (see also Ref. [17]). Therefore a new critical dimension 1 < d ′ c < 2 must exist, below which there appears a dynamic phase transition to a power-law inactive phase, see Fig. 7.
Within the framework of a truncated loop expansion to one-loop order (Sec. V), we may describe the above scenario as a function of dimension d at least qualitatively; in this approximation, we find d ′ c = 4/3. Although the structure of the phase diagram seems to be correct, the ensuing one-loop critical exponents in this rather uncontrolled approximation are well off the numerically found values [17,19,22,23]. Yet unfortunately, we were unable to improve on our approximations, e.g. by extending the calculations to the two-loop level. Also, the emergence of dangerous irrelevant variables precludes the firm establishment of scaling relations [22] relating exponents in the active phase, say, and those at the critical point. Although so far there has been no evidence for this, the possibility of violations of such scaling laws should be carefully and thoroughly reinvestigated in numerical simulations. It would also certainly be highly desirable to possibly design a different renormalization scheme to address the even m universality class. The desirability of going beyond the one-loop approximation is even clearer in the inactive phase, where the density is expected to decay as A/(D R t) d/2 , for d < 2, and D R is a renormalized diffusion constant. To one-loop order, no such renormalization takes place. It is therefore not possible to establish postulated scaling laws relating the manner in which D R is supposed to behave near the critical point to the other exponents defined in the active phase. The inactive phase is also interesting in that it is expected to exhibit a kind of spontaneous symmetry breaking as a result of the parity invariance of the dynamics. That is, if we start with a state with a even number of particles (equivalent to imposing periodic boundary conditions on the corresponding Ising model), the asymptotic state is the trivial absorbing state (fully ordered in the Ising language). But if we begin with an odd number of particles this cannot happen: instead we expect the system to evolve to a state with a finite average total number of particles. Once again, one may postulate scaling laws relating the power law divergence of this number as the critical point is approached to other exponents, and it would be very useful to have a more systematic derivation of these, which is so far lacking in our analysis. The study of the nature of the inactive phase seems to be particularly interesting, as numerical studies appear to indicate the existence of unexpected simplifications [32].
As discussed in Sec. VI, an obvious generalization to N > 2 particle 'flavors' leads to qualitatively new behavior, governed by the exactly computable exponents of the limit N → ∞. To our knowledge, this case of BARW with even m and N > 2 particle types has not been studied in simulations, nor have BARW with k = 3 (see Sec. VIII).
FIG. 2. One-loop diagrams for the renormalizations of (a) the annihilation vertex λ, (b) the branching vertex σm, and (c) the pair production vertex τ ; (d) the generation of the branching process ∝ σm−2 through a combination of branching ∝ σm and two-particle annihilation.   9. Most singular (bubble) diagrams for the renormalization of (a) the particle decay rate µ, (b) the branching rate σ1, and (c) the propagator (mass) µ + σ for BARW with m = 1 near dc = 2 dimensions.  .