{'fraction_non_alphanumeric': 0.057623585467540206, 'fraction_numerical': 0.02036926742108398, 'mean_word_length': 3.922456757549106, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Persistence properties of a system of coagulating and annihilating random walkers

We study a d-dimensional system of diffusing particles that on contact either annihilate with probability 1/(q-1) or coagulate with probability (q-2)/(q-1). In 1-dimension, the system models the zero temperature Glauber dynamics of domain walls in the q-state Potts model. We calculate P(m,t), the probability that a randomly chosen lattice site contains a particle whose ancestors have undergone exactly (m-1) coagulations. Using perturbative renormalization group analysis for d<2, we show that, if the number of coagulations m is much less than the typical number M(t), then P(m,t) ~ m^(z/d) t^(-theta), with theta=d Q + Q(Q-1/2) epsilon + O(epsilon^2), z=(2Q-1) epsilon + (2 Q-1) (Q-1)(1/2+A Q) epsilon^2 +O(epsilon^3), where Q=(q-1)/q, epsilon =2-d and A =-0.006. M(t) is shown to scale as t^(d/2-delta), where delta = d (1 -Q)+(Q-1)(Q-1/2) epsilon+ O(epsilon^2). In two dimensions, we show that P(m,t) ~ ln(t)^(Q(3-2Q)) ln(m)^((2Q-1)^2) t^(-2Q) for m<<t^(2 Q-1). The 1-dimensional results corresponding to epsilon=1 are compared with results from Monte Carlo simulations.


I. INTRODUCTION
Persistence is understood as a property of an evolving system to "remember" its initial configuration for anomalously long times. A particular case of persistence that has received much attention is that of site persistence (see [1] for a review). The site persistence probability is defined as the probability that the values of a dynamical variable at a given set of sites do not change up to time t. For instance, in a spin system this could be the probability that a spin at a given site does not flip up to time t, or in a reaction-diffusion system, the probability that no reaction takes place at that site up to time t. In many cases the site persistence probability decays at large times as a power-law [2].
A natural generalization of site persistence is persistence of a pattern present in the initial configuration [3]. An instance of pattern persistence would be the survival of a test particle in a random environment. Examples of the random environment include diffusing traps [4,5], reaction-diffusion systems such as A∅ ↔ ∅A [6], A+ B → ∅ [6,7] or A i + A j → A i+j with mass dependent diffusion rates [8,9,10], and predators in predator-prey models [11,12]. The problem of calculating the survival probability of the test particle is hard, mainly because in the rest frame of the test particle, the motion of the other particles is correlated.
The 1-dimensional Potts model has been a testing ground for various concepts of persistence. The site persistence problem mentioned above, has been exactly solved for the 1-dimensional q-state Potts model evolving via zero temperature Glauber dynamics [13]. Besides site persistence, several other persistence properties of the Potts model have been studied. Among these are the probability that a domain wall has never encountered another domain wall [3,14,15], and the probability that a domain present in the initial configuration survives up to time t [16]. The former problem has been studied numerically [3,15], by mean field approximations [3] and perturbatively near q = 1 [14]. However, the results obtained by these techniques do not approximate well the numerical results in the whole range of q.
In dimensions greater than one, the dynamics of domain walls in the Potts model is difficult to treat analytically. Instead, we note that in 1-dimension, the zero temperature Glauber dynamics of the q-state Potts model is equivalent to a system of diffusing particles that on contact either annihilate with probability 1/(q − 1) or coagulate with probability (q − 2)/(q − 1) [17,18,19]. We study the persistence properties of this reaction-diffusion system in an arbitrary number of dimensions using the renormalization group method and calculate the exponents as an ǫ-expansion.
The question that we ask is, given this reactiondiffusion system, what is the fraction of particles that have never encountered another particle up to time t?
More generally, what is the fraction of particles whose ancestors have undergone m coagulations up to time t? A convenient way to keep track of the history of coagulations is to assign a mass to each particle as follows. At time t = 0, let all particles be of mass one. Each time two particles coagulate, the new particle has a mass which is the sum of the masses of the two parent particles. It is clear that the particles of mass m will be those whose ancestors have undergone exactly (m − 1) coagulations.
LetP (m, t) be the probability that a randomly chosen site at time t contains a particle of mass m. LetN (t) and ρ(t) denote the average particle density and average mass density respectively. Then, the probability distribution P (m, t) is expected to have the scaling form where l m = (m/ρ(t)) 1/d is the length scale associated with mass. The large time behavior ofρ(t) andP (m, t) is characterized by two exponents δ and θ. The mass densityρ(t) ∼ t −δ . For masses much smaller than the typical mass,P (m, t) decays as a power law in time as m ζ/d t −θ . We will call θ the persistence exponent. The exponent ζ characterizes the small x behavior of the scal- , where we have used the fact thatN (t) ∼ t −d/2 in dimensions less than two [20]. The two independent exponents θ and δ are known for some limiting cases. When q = 2, the model reduces to the reaction diffusion model A + A → ∅. All particles are of mass one and it is known thatP (1, t) ∼ t −d/2 for d < 2 and P (1, t) ∼ ln(t)/t for d = 2 [20,21]. Hence, δ = d/2, θ = d/2 and ζ = 0 for q = 2. When q = ∞, the model is equivalent to the reaction-diffusion system [10,18,22,23,24,25]. Since mass is conserved, δ = 0. It has been shown that θ = d+ǫ/2+O(ǫ 2 ), and ζ = ǫ + O(ǫ 2 ) for ǫ = 2 − d > 0 [10]. In 2-dimensions P (m, t) ∼ ln(m) ln(t)/t [10]. In 1-dimension, it is known via an exact calculation thatP (m, t) ∼ mt −3/2 [22]. When q ≈ 1, θ has been calculated perturbatively to be θ = (q − 1)3 √ 3/(2π) + O((q − 1) 2 ) [14]. However for arbitrary values of q, the only known analytical result follows from a mean field approximation [3]. But the numerically obtained value for θ is very different from the mean field value. In this paper, we address this issue by using the renormalization group formalism to systematically calculateP (m, t) for arbitrary q.
We now summarize our main results and give an outline of the rest of the paper. In Sec. II, we give a precise definition of the model and derive the stochastic partial differential equations obeyed by the mass distribution. In Sec. III we express the exponent δ in terms of the exponent θ, though at a different value of q, reducing the number of unknown exponents to one. We show that Thus, In Sec. IV we use the technique developed in [10] to calculate the persistence exponent θ as an ǫ-expansion, where ǫ = 2 − d > 0. We show that given that t → ∞ and m ≪ M (t), where M (t) is mass of a typical particle at time t. The analytical results for θ and δ in 1-dimension obtained by putting ǫ = 1 are compared with the results from numerical simulations.
In Sec. V we show that the coefficient of ǫ n in Eq. (4) is a polynomial of degree 2n in the variable Q = (q − 1)/q. This observation allows us to calculate the 2-loop corrections to the exponent ζ to be where A = −0.006 . . .. The analytical results for ζ in 1dimension obtained by putting ǫ = 1 are compared with the results from numerical simulations. Finally, we conclude with a summary and discussion in Sec. VI.

II. MODEL AND FIELD THEORETIC FORMULATION
In this section, we define the model and derive the stochastic partial differential equation obeyed by the mass distribution. Consider a d-dimensional lattice whose sites may be occupied by particles that possess a positive integral mass. Multiple occupancy of a lattice site is allowed. Given a certain configuration of particles on this lattice, the system evolves in time via the following microscopic moves. (i) With rate D, each particle hops to a nearest neighbor lattice site. (ii) With rate λ c , two particles at the same site coagulates together to form a new particle whose mass is the sum of the masses of the two parent particles. (iii) With rate λ a , two particles at the same site annihilate each other. To make connection with the model discussed in the introduction, we have to choose λ c = λ(q − 2)/(q − 1) and λ a = λ/(q − 1) where λ is a reaction rate. The limit λ → ∞ corresponds to instantaneous reactions. In dimensions d ≤ 2 and in the limit of large time, the statistical properties of a finite reaction rate particle system were shown to be equivalent to those of a system with infinite reaction rates [10]. However, from the field theoretic point of view, it is more convenient to work with finite reaction rates, and hence λ will be taken to be finite in this paper.
Starting from the master equation for the time evolution of the system, we now derive the effective field theory of the model. Let where the time dependence of P has been dropped for notational simplicity and {n i,m + 1} denotes the configuration (n i,1 , n i,2 , . . . , n i,m +1, . . .) at site i. The first term in the right hand side of Eq. (7) describes the loss and gain terms arising from particles diffusing to their nearest neighbors with rate D. The second term describes the loss and gain terms due to the coagulation of a pair of particles at a site with rate λ c to form a new particle whose mass is the sum of the constituents. The third term describes the loss and gain terms due to annihilation of a pair of particles at a site with rate λ a . The field theory corresponding to the problem can be derived from the master equation using Doi's formalism [26]. In short, regarding the master equation as a Schroedinger equation in imaginary time, the functional integral representation of the corresponding non-Hermitian evolution operator is constructed. This allows one to write down a functional integral expression for any correlation function of the problem, includingP (m, t). After taking the continuum limit, one is left with the problem of solving an interacting field theory. The application of Hubbard-Stratonovich transformation to this field theory leaves one with a stochastic partial differential equation. We refer to Refs. [25,27,28] for reviews of this procedure. Following this procedure, solving the master equation Eq. (7) is equivalent to solving the following Langevin equation for a stochastic fieldP ( x, m, t): , ξ is white noise in space and time with unit standard deviation and i 2 = −1.
The stochastic fieldP ( x, m, t) is complex and is different from the local mass distribution P ( x, m, t), which denotes the number of particles of mass m in the volume d d xdm at time t. However, the moments of P are related to the moments ofP (for instance, see [28,29]). 2 , and so on, where the bar on top denotes an averaging over noise, and ∆x and ∆m are lattice cutoffs. In this paper we only study the first moment of P ( x, m, t), and hence disregard the difference between P ( x, m, t) andP ( x, m, t) in the rest of the paper.
We will be studying the behavior of the following three quantities, where N ( x, t) is the local particle density, ρ( x, t) is the local mass density and M (t) ∼ρ(t)/N (t) is the typical mass at time t. The time evolution equations obeyed by N ( x, t) and ρ( x, t) are easily obtained from Eq. (8). As for P (m, x, t), for m ≪ M (t), we neglect the convolution term in the right hand side of Eq. (8). This approximation is justified because at large times the probability of collision of two light particles is negligible compared to the probability of collision of light and heavy particles. The resulting equations are: Note that the dependence of P on mass is no longer governed by Eq. (13). Once the time dependence ofP is calculated, its mass dependence can be restored using dimensional analysis. In the rest of the paper, for the sake of notational simplicity, we omit the dependence of P on mass, unless there is a cause for confusion. Equations (12), (13) and (14) can be simplified as follows. Let for some parameter λ. Rescaling the local particle density, local mass distribution and average density according to brings Eqs. (12)-(14) into the following form: where It can be shown that Eqs. (18) and (19) describe the two Qλ → ∅, in the limit when the concentration of A-particles is much greater than concentration of B-particles.

III. SCALING ANALYSIS OF STOCHASTIC EVOLUTION EQUATIONS.
In this section, we obtain some exact results for the model. First, the scaled density of particles N ( x, t) obeys the same equation as the particle density in the A + A → A reaction. For this reaction, the density of particles decays for large times as t −d/2 in d < 2 [20,21]. Thus, where c is a constant depending on dimension only. This is a generalization of the exact 1-dimensional result [30,31].
Secondly, there is a relation between the local mass distribution of light particles P and the local mass density ρ. Under the substitution Q → (1 − Q), or equivalently q/(q − 1) → q, Eq. (19) transforms into Eq. (20). Therefore, if F P (Q, x 1 , t 1 , x 2 , t 2 , . . .) is a correlation function of P -fields, which is independent of initial conditions, then is the correlation function of the same configuration of ρ-fields. In particular, sincē We now examine Eqs. (18) and (19) for special values of the parameter Q. When Q = 0, the non-linear term in Eq. (19) vanishes and the concentration of monomers is conserved on average. Therefore, θ = 0 for Q = 0. [10]. If d = 1, then θ = 3/2, which is a consequence of an exact solution [22], rather than the ǫ-expansion cited above. Collecting these results together, we have The question remains as to whether θ for Q < 1/2 or equivalently q < 2 has any physical implication. For the site persistence problem in 1-dimension, it is known that the site persistence probability of the q-state Potts model maps to an Ising system with an initial magnetization given by 2/q −1, evolving via zero-temperature Kawasaki dynamics [18,19,32]. The latter system is defined at any value of q > 1. The correspondence between these two models also holds for the persistence probability of a single domain in the Potts model [16]. It would be interesting to understand what quantity, if any, in the Ising model corresponds to the survival probability of domain walls in the Potts model.
In this section, we calculate the large time behavior of P (t) using the formalism of perturbative renormalization group. We closely follow the solution of the A i + A j → A i+j model presented in Ref. [10] .
The solution toP (t) as a perturbative expansion in powers of λ can be constructed from Eqs. (18) and (19) using Feynman diagrams [33]. The Feynman rules for constructing terms of the expansion are summarized in Fig. 1. Diagrammatically,P andN are the sum of all Feynman diagrams with one outgoing P and N line respectively. Clearly, there are an infinite number of diagrams contributing toP andN . These diagrams can be grouped together according to the number of loops that they contain, thus giving rise to the loop expansion. Let ǫ = 2 − d. The contribution from each diagram is a function of the dimensionless terms λN 0 t and g(t) = λt ǫ/2 and a term that gives the correct physical dimension [(λt) −1 forN and (λt) −2 forP ]. A simple combinatorial argument shows that the contribution from a diagram with n loops is proportional to g(t) n [29]. When ǫ < 0, the main contribution toP andN comes from properly renormalized tree level diagrams (diagrams without loops) [34]. When ǫ > 0, the loop expansion fails since for large times g(t) is no longer a small perturbation parameter. We therefore conclude that 2 is the upper critical dimension.
For d < 2 we will use the formalism of perturbative renormalization group to convert the loop expansion into an ǫ-expansion and calculate scaling exponents as a series in ǫ.

A. Tree level diagrams
LetN mf andP mf be mean field densities given by the sum of contributions coming from tree diagrams with a single outgoing N -line and P -line respectively. We de-noteN mf andP mf by thick solid lines and thick dashed lines respectively. The integral equations satisfied byN mf andP mf are presented in diagrammatic form in Fig. 2(a) and 2(b). After differentiating with respect to time, they can be written in analytic form as in which one can easily recognize the Smoluchowski rate equations of the model, obtained from Eqs. (18) and (19) by neglecting the noise terms in the right hand side. Equations (25) and (26) are easily solved yieldinḡ From Eq. (28), we obtain where θ mf is the mean field answer for θ.
In calculating loop corrections to any given order, we are faced with the problem of summing over infinitely many diagrams containing a given number of loops. This problem can be simplified by introducing mean field propagators which are sums of all tree diagrams with one incoming line and one outgoing line. Expressed in terms of these mean field propagators, there are only finitely many diagrams with a fixed number of loops.
Let G NN mf and G PP mf be mean field propagators. The integral equations satisfied by them are presented in diagrammatic form in Figs. 2(c) and 2(d). The solutions to these equations are where 1 = ( x 1 , t 1 ), 2 = ( x 2 , t 2 ) and G 0 is the Green's function of the linear diffusion equation.

B. One loop diagrams
Using the mean field propagators and densities, it is easy to classify all one loop diagrams contributing tō P (t). These are shown in Fig. 3. The computation of the corresponding Feynman integrals is straightforward. The contributions from one-loop diagrams in the limit where (a), (b) and (c) refer to the contributions from diagrams in Figs. 3(a), 3(b) and 3(c) respectively. Adding these one-loop contributions to the mean field answer Eq. (28), we obtain in the limit N 0 → ∞, ǫ + 6 − 2Q(ǫ + 2) (ǫ + 2) 2 (ǫ + 4) + 2-and higher loop corrections, (35) where A = P 0 /(N 0 λ) 2Q .

C. Renormalization Group Analysis of the Model.
The large time asymptotic behavior ofP (t) can be obtained by solving the Callan-Symanzik equation with initial conditions given by Eq. (35) (see Ref. [27] for a review). The coefficients of Callan-Symanzik equation are determined by the law of renormalization of all the relevant couplings of the theory Eqs. (18) and (19). Power counting analogous to that carried out in [10], shows that there are only two relevant couplings of the theory in d < 2: the reaction rate λ and the initial mass distribution P 0 . We will derive the one-loop renormalization law of the initial mass distribution by requiring that Eq. (35) is non-singular in the limit ǫ → 0 if expressed in terms of renormalized relevant couplings.
Let t 0 be a reference time and g 0 = λt ǫ/2 0 be the dimensionless reaction rate. We choose t 0 in such a way that g 0 ≪ 1. The mechanism of renormalization of the reaction rate in the theory is identical to that of the reaction A + A → A. Physically, the renormalization of reaction rate is explained by the recurrent property of random walks. The probability of a reaction between particles at time t is proportional to the "bare" reaction rate, multiplied by the probability that the reaction has not occurred before time t. In d ≤ 2, the latter probability explicitly depends on time t. The law of renormalization of the reaction rate has been worked out in Ref. [20]. If g R is the renormalized reaction rate, then it is related to g 0 by the relation where g * = 2πǫ + O(ǫ 2 ) is the nontrivial fixed point of the renormalization group flow in the space of effective coupling constants. The mass distributionP (t 0 ) can now be expressed in terms of the renormalized reaction rate g R to bē .
(37) The order g R term in Eq. (37) is singular at ǫ = 0. To cancel this divergence, we have to introduce a renormalized initial mass distribution P R : where Z(g R , t 0 , ǫ) is chosen such that is non-singular at ǫ = 0. Substituting Eq. (39) into Eq. (37), we obtain The Callan-Symanzik equation is obtained by noting thatP (t, λ, P 0 , ǫ) does not depend on the reference time t 0 . Therefore, It follows from dimensional analysis that the most general form of the average mass distribution is where Φ is a dimensionless function. Using the scaling function Eq. (42) in Eq. (41), we obtain the Callan-Symanzik equation forP (t): are the beta and gamma functions of the theory.
At large times, the solutions of Eq. (43) are governed by the nontrivial fixed point g * of the beta function. It then follows from the Callan-Symanzik Eq. (43) that Thus we obtain from Eqs. (45) and (46) that The knowledge of θ to the first order in ǫ combined with Eqs. (2) and (3) allows one to calculate the exponents δ and ζ with the same precision: The exponent ζ is proportional to the sum of the anomalous dimensions of P and ρ. As a result, the mass dependence ofP (m, t) can be captured neither by mean field theory nor by Smoluchowski approximation (see [10] for a more detailed discussion of this point).
The results of this subsection can be summarized as follows. The mean mass distributionP (m, t) varies as The upper critical dimension of our model is 2. The non-trivial fixed point of the β-function Eq. (44) vanishes at d → 2. We therefore expect the mean field answers Eqs. (27) and (28) to give the correct large time-small mass of average densities in two dimension, modulo logarithmic corrections. In this subsection we calculate these corrections.
When Q = 1 it was shown that in two dimensions, P (m, t) ∼ ln(t) ln(m)t −2 for t → ∞, m ≪ M (t) [10]. To calculate these corrections for arbitrary Q we need to solve the Callan-Symanzik equation (43) with coefficients calculated at d = 2. In two dimensions, Then Eq. (43) reduces to (53) which has to be solved with the initial condition provided by the mean field theory. The solution to Eq.(53) with this initial condition is ) .
(55) When Q = 1, we recover the result of Ref. [10]. When Q = 0, P (t) ceases to depend on time, as it should. When Q = 1/2,P ∼ ln(t)/t, which coincides with the decay law of the concentration of particles in A + A → ∅ reaction [20].
The dimensional arguments that led to Eq. (3) cannot capture the mass dependent logarithmic corrections that are present in 2-dimensions. Hence, we need to generalize these dimensional arguments. This is provided by the Callan-Symanzik equation obeyed byP (m, t) when considered as a function of both m and t.
The full distributionP (m, t) cannot depend on the choice of reference time t 0 , which we introduced to regularize the perturbative expansion ofP (t). Therefore, From dimensional analysis, it follows that The form in Eq. (57) is different from the scaling form used in Eq. (42) because P R has to be now expressed in terms of m. Substituting Eq. (57) into Eq. (56) we obtain are scaling dimensions of fields N and ρ, which can be obtained from the corresponding loop expansions.
We look for solutions of Eq. (58) of the form The time dependent function F 1 obeys the Callan-Symanzik equation (43). Using this fact and substituting Eq. (59) into Eq. (58), we obtain where Γ(g R ) = (2Q − 1)g R /(2πd) + O(g 2 R ). As expected, when d < 2, Q > 1 2 and t → ∞ the solution of Eq. (60) is given by Eq. (50). Let d = 2. Solving Eq. (60) for Q ≥ 1 2 with the initial condition F 2 (m 0 ) = const, provided by mean field theory, we find that is a reference point in mass space. Combining Eqs. (55) and (61), we conclude that in for m 0 ≪ m ≪ M (t) and t → ∞. For Q = 1 we recover the answer for the average mass distribution in the A i + A j → A i+j model obtained in Ref. [10]. This result has also been verified numerically [10]. When Q = 1/2, P (m, t) no longer depends on mass, as expected.

E. Comparison with results from Monte Carlo simulations
In this subsection, we compare the results obtained for the exponents θ and δ as an ǫ-expansion with results from Monte Carlo simulations in 1-dimension. Let θ 1 denote the value of θ obtained by truncating the ǫ-expansion at order ǫ and setting ǫ = 1. Then, In Table I, we compare this analytic expression with results from numerical simulations (see columns 2 and 3). There is good agreement.
To go beyond the expression in Eq. (63) and to make an estimate of the error arising by neglecting terms of order ǫ 2 and higher, we proceed as follows. Let the corrections from order ǫ 2 and higher orders be denoted by R(ǫ, Q), such that R(ǫ, Q) must vanish at Q = 0 and Q = 1/2 (see Eq. (24)). Moreover, R(1, 1) = 0, since θ = 3/2 when Q = 1 and ǫ = 1 [22]. Therefore, where h 1 and h 2 are unknown functions. Setting ǫ = 1, we obtain The value of function h 1 (1, Q) at Q = 0 can be determined. It was shown in Ref. [14] that Therefore, h 1 (1, 0) = 3 √ 3/π − 1. A two loop calculation carried out in Sec. V shows that the function h 1 (1, Q) is slowly varying in the interval Q ∈ [0, 1]. Therefore, we replace the function h 1 (1, Q) by its value at Q = 0 and denote the resulting expression as θ 2 . Thus, we obtain In Table I, we compare θ 2 with results from Monte Carlo simulations in 1-dimension. The error decreases as compared to θ 1 . The error due to dropping terms of order ǫ 2 and higher can be estimated. The function h 1 (1, Q) in Eq. (66) is or order 1. The function |Q(Q − 1/2)(Q − 1)| takes on a maximum value of 0.05 . . . in the interval Q ∈ [0.5, 1]. Hence the absolute error is of order 0.05, which is in agreement with the results presented in Table I.
We do a similar analysis for δ. Let δ 1 be the value of δ obtained by truncating the series Eq. (48) and then putting ǫ = 1. Then To obtain δ 2 , we substitute Q → (1 − Q) in Eq. (68) to obtain In Table II, we compare the results δ 1 and δ 2 with results from numerical simulations. Very good agreement is seen.

V. THE ANALYSIS OF TWO-AND HIGHER LOOP CORRECTIONS.
A. General structure of the loop expansion.
In this subsection, we examine the contributions from diagrams with 2-and more loops. It will be shown that the coefficient of ǫ n in the ǫ-expansion of θ is a polynomial of degree 2n in Q. It is easier to derive the result, not by using the formalism of renormalization group, but by identifying the principal set of diagrams contributing to the large time limit ofP (t) and deriving a simple integral equation satisfied by the sum of these diagrams.
The polarization operator Π(t 2 , t 1 ) is defined as the sum of all one-particle irreducible diagrams with one outgoing and one incoming P -line, with the external propagator lines stripped off. Using the polarization operator, we can write down the Schwinger-Dyson equation obeyed byP (t). LetP (t) and Π(t 2 , t 1 ) be denoted by a thick dashed-dotted line and by a grey circle respectively. Then,P (t) satisfies the equation shown diagrammatically in Fig. 4(a). In equation form, it is In terms of η, Eq. (71) reduces to where η 0 is a constant independent of t. Differentiating with respect to t, we obtain  The expansion of Π(t, t 1 ) in terms of Feynman diagrams is shown in Fig. 4(b). From the Feynman diagrams in Fig. 1, it follows that the number of dotted lines in any given diagram is conserved as one moves from the right to the left. Any diagram contributing to the polarization operator has a single dotted line threading through it from right to left. Also, the only vertex that contributes a factor Q is the PPN vertex. Therefore, the Q-dependent part of a given diagram has the form Q #of PPN vertices i (t i /t i+1 ) 2Q , where the product is over all vertices involving the P -line. For a kloop diagram, we can have utmost 2k vertices of the type PPN. Also, it was shown in Ref. [10] that n-loop diagrams contribute at the order ǫ n only. Hence, after coupling constant renormalization, one finds that the expression in square brackets of Eq. (73) is of the form t −2 ∞ n=1 ǫ n P 2n (Q), where P 2n (Q) is a polynomial of degree 2n in Q, and where the factor t −2 has been pulled out to give the right dimension.
We can now solve Eq. (73) perturbatively, order by order in ǫ. Simple dimension counting shows that t 2Q Π(t, t 1 )t −2Q is a dimensionless function. The previous argument shows that F (τ ) = ∞ n=1 ǫ n P 2n (Q). Assume a power law solution for η, i.e., η = ct −θp where θ p = n=1 a n ǫ n and c is a constant. Then Eq. (73) simplifies to Expanding F (τ ) and θ p as series in ǫ, we obtain ∞ n=1 a n ǫ n = − (75) Solving Eq. (75) order by order for a n , it is easy to verify that the coefficient of ǫ n is a polynomial of degree 2n in Q, i.e., where C n,p 's are some unknown constants. Given Eq. (76), it is easy to rederive the one loop correction to θ (Eq. (47)) obtained by the renormalization group formalism. The three unknowns in the coefficient of ǫ in Eq. (76) are obtained from the exact results in Eq. (24) giving C 1,0 = 0, C 1,1 = −3/2, and C 1,2 = 1.
As the mean field answer for the exponent ζ is 0, it is desirable to evaluate order ǫ 2 correction to Eq. (49). This requires the knowledge of order ǫ 2 term in θ.
From Sec. V A, we know that the term in ǫ 2 is a polynomial of degree 4 in Q, i.e., Out of the 5 unknown C 2,k 's, two are fixed by the conditions that θ = 0 for Q = 0 and θ = d/2 for Q = 1/2 (see Eq. (24). For Q = 1, it is known that θ = d + ǫ/2 + O(ǫ 2 ) [10]. We assume that the order ǫ 2 term is absent when Q = 1 (see Appendix for a heuristic validation of this assumption). This fixes the third constant and we are left with where A and B are constants.
The constant A is not difficult to calculate. The contribution to θ of order ǫ 2 Q 4 comes from the square of 1-loop polarization operator and from two-loop diagrams with four P P N -vertices. There are only two such diagrams, which are shown in Fig. 5. The numerical computation of corresponding Feynman integrals gives where the first term on the right hand side comes from the order-ǫ 2 term in the one-loop polarization operator. The calculation of constant B seems an almost impossible task, as there are over twenty two-loop diagrams contributing to it. However, the terms proportional to B drop out of the the two loop expression for ζ. Substituting Eq. (78) into Eq. (3) one finds that In Table III, we compare the one-loop expression for ζ (Eq. (49)) and two-loop expression for ζ (Eq. (80)) with numerical results in 1-dimension.

VI. SUMMARY AND CONCLUSIONS
In summary, we develop a systematic method to calculate the persistence exponent θ for a system of coagulating and annihilating random walkers, in arbitrary dimensions. In 1-dimension, this corresponds to persistence probabilities of domain walls in the Potts model evolving via zero temperature Glauber dynamics. We establish an exponent relation by which the number of unknown exponents in the problem is reduced from two to one. The unknown persistence exponent θ is determined perturbatively using the formalism of renormalization group.
The persistence problem studied in this paper can be considered as a special case of a more general problem of the survival probability of a test particle with diffusion constant κ times the diffusion constant of the other particles. In this case it is known that the persistence exponent θ κ (q) depends on κ [3,14]. While simple limiting cases have been studied [14,35,36,37] via numerics, mean-field or perturbative techniques, a general understanding is still lacking. It would be interesting to extend the formalism of perturbative renormalization group to calculate θ κ (q).

APPENDIX A: MASS DISTRIBUTION IN THE
Ai + Aj → Ai+j MODEL In this appendix, we present a heuristic derivation of the distribution of small masses in the A i + A j → A i+j model. This model corresponds to the Q = 1 limit of the model discussed in the paper. For the A i + A j → A i+j model, it is known [10,22] that for m ≪ t d/2 , In this appendix, we give a heuristic argument as to why the terms of order ǫ 2 and higher could be absent, as a result of which the expansion up to order ǫ gives the exact answer.
Putting λ a = 0 in Eq. (8), we obtain that the mass distribution P (m, x, t) evolves according to obeys the equation The function B(s, t) obeys the same equation as the density N [25]. However, the initial conditions at t = 0 are different. If the initial density of particles N (0) = N 0 , then B(s, 0) = N 0 (1 − e −s ). Therefore, if the average particle density isN (N 0 , t), then F (s, t) =N (N 0 , t) −N (N 0 (1 − e −s ), t). (A7) The function F (s, t) will have the scaling form where the scaling function g(x) ∼ x −φ for x ≫ 1. On performing the inverse Laplace transform, we obtain P (m, t) ∼ m φ−1 t d/2(1+φ) .
To calculate φ, we look at the behavior of the particle densityN . It is expected to have the scaling form where the scaling function h(x) behaves for large x as where ψ is some exponent greater than zero. Substituting Eqs. (A10) and (A11) into Eq. (A7), it is straightforward to verify that φ = ψ. (A12) In 1-and 2-dimensions, it is easy enough to verify that ψ = 2 and ψ = 1 respectively, consistent with the exact results in Eq. (A1). In other dimensions, we argue as follows. The density of particles in the A i + A j → A i+j model is inversely proportional to the area swept out by a random walker in time t. This area varies as ( √ t + c) d , where c is some constant. Then, the particle density decays asN ∼ t −d/2 (1 − const/t), such that Substituting into Eq. (A9), we obtain which is the one-loop answer in Eq. (A1).