Asymptotic behavior of branching diffusion processes in periodic media

We study the asymptotic behavior of branching diffusion processes in periodic media. For a super-critical branching process, we distinguish two types of behavior for the normalized number of particles in a bounded domain, depending on the distance of the domain from the region where the bulk of the particles is located. At distances that grow linearly in time, we observe intermittency (i.e., the $k$-th moment dominates the $k$-th power of the first moment for some $k$), while, at distances that grow sub-linearly in time, we show that all the moments converge. A key ingredient in our analysis is a sharp estimate of the transition kernel for the branching process, valid up to linear in time distances from the location of the initial particle.


Introduction
Consider a collection of particles Y 1 (t), Y 2 (t), . . . in R d that move diffusively and independently according to where W k denote independent Brownian motions in R d . Each particle independently branches into two particles or is annihilated at rates the depend on its location: a particle at x ∈ R d branches into two particles at rate α(x) ≥ 0, and is annihilated at rate β(x) ≥ 0. The newly created particles starting at the location of their parent then repeat this process independently of each other. This process is referred to as a d-dimensional branching diffusion process. We suppose that the drift b(x), the non-degenerate diffusion matrix σ(x), and the rates α(x) and β(x) are all Lipschitz continuous and Z d periodic (and thus bounded). That is, b(x + k) = b(x) for all x ∈ R d and k ∈ Z d , and similarly for σ, α and β. In addition, we assume that all the matrix a(x) is C 1 (R d ).
The main topic of interest here is the limiting behavior of branching diffusion processes in periodic media in the supercritical regime. Our main goal is to study the distribution of the number of particles in regions whose spatial location depends on time. With probability that tends to one, the entire population is confined to a region that grows linearly in time (see Chapter 7.3 in the book of Freidlin [10]). The effective drift of a branching process can be understood heuristically as the speed at which the bulk of the particles is traveling in space. We will give a precise definition of the effective drift later in Section 2. For a bounded region at a fixed location, assuming that the effective drift is zero, the structure of the population is similar to that in the compact setting. See, for example, Engländer, Harris, Kyprianou [8] and references therein. For a time dependent region inside the linearly growing front, the normalized number of particles converges almost surely (see, for example, Uchiyama [27] in the case of constant coefficients). The nature of this convergence, however, depends on how distant the region is from the location of the initial particle (assuming for simplicity that the effective drift is zero). At linear in time distances, we will show that intermittency may occur (i.e., the k-th moment dominates the k-th power of the first moment for some k), while, at distances that grow sub-linearly in time, we will prove that all the moments converge. For the case of homogeneous media and for the case of compactly supported branching term, this question has been studied in the work of Koralov [17] as well as Koralov, Molchanov [18].
Given a single particle initially at x ∈ R d , the transition kernel u(t, x, y) is defined by where f ∈ C b (R) and the sum is over all particles alive at time t ≥ 0. The function (t, y) → u(t, x, y) satisfies (2) ∂ t u = L x u, x, y ∈ R d , t > 0, with initial condition u(0, ·, y) = δ y (·), where L is the operator a(x) = σ(x)σ * (x), and r(x) = α(x) − β(x). The operator L − r(x) is the generator of the process (1). The first step in our analysis is a precise asymptotic description of the transition kernel u(t, x, y), valid up to the large deviation scale, that is, for x−y = O(t).
There are two main parts in the asymptotic analysis of u(t, x, y). First, we transform the operator L in order to alter the effective drift of the process, while simultaneously turning the branching rate into a constant. Thus, the problem reduces to studying the transition kernel of an altered diffusion process near the diagonal, where x−y = O( √ t). The next part is to prove a local limit theorem for the new transformed kernel at this diffusive scale.
The ingredients we use to obtain the asymptotics of the transition kernel -exponential change of measure, homogenization and local limit theorems for the resulting diffusion process are fairly standard. In spite of this, the precise asymptotics of the transition kernel that holds up to linear in time distances has not been published, as far as we know (in 2007, Agmon gave a talk [1] where this result was announced). Here, we provide a simple probabilistic proof that establishes uniform asymptotics of the transition kernel for d-dimensional second-order parabolic operators with periodic coefficients. The precise asymptotics in the 1-dimensional case has been obtained previously by Tsuchida in [26].
Prior results in this direction, in d dimensions, give estimates of the heat kernel, as opposed to precise asymptotics. The seminal work of Aronson [2] gives global estimates on the heat kernel, while in [23] Norris proves a generalization of Aronson's Gaussian bounds in the case of periodic coefficients and identifies an effective drift of the heat flow. The upper and lower bounds of Norris [23] have different constant prefactor in front of the Gaussian term, although the logarithmic asymptotics are sharp. We provide a stronger result that correctly identifies the main term of the asymptotic expansion of the transition kernel, which is precise up to the domain of large deviations (up to distances in space that are linear in time). The asymptotics of Green's function for the corresponding elliptic problem for different values of the spectral parameter has been studied extensively (see, e.g., Murata, Tsuchida [22], Kuchment, Raich [19]).
The asymptotics proved in Section 2 plays a crucial role in analyzing the behavior of the branching diffusion process in periodic media, in Section 3. The bulk of the particles will be seen to be located around ofvt wherev denotes the effective drift of the process (defined later at (16)). Let n y (t, x) denote the number of particles located in a unit ddimensional cube containing y ∈ R d , assuming that, initially, there is one particle located at x ∈ R d . In Section 3.1, for a super-critical branching process, we study the asymptotic behavior of n y (t, x) in the domain of large deviations, that is when y −vt = O(t). We observe the effect of intermittency, that is, for each vector v ∈ R d , v =v, there exists k ≥ 2 such that the k-th moment of n vt (t, x) grows exponentially faster than the k-th power of the first moment. This result was first proved in [18] in the case of a supercritical branching diffusion process in R d with identity diffusion matrix, zero drift, and a positive constant potential. Here, in contrast to [18], we do not have explicit expressions for the transition kernel, but only have asymptotic formulas. This makes the analysis of the higher order moments more involved.
In Section 3.2, we define a sequence of periodic functions f k (x) that serve as limits for the k-th moments of N(t, x)/E(N(t, x)), where N(t, x) denotes the total number of particles in R d , assuming that, initially, there is one particle located at x ∈ R d .
In Section 3.3, we again study n y (t, x), but here we assume that y −vt = o(t). That is, we study the distribution of particles near the region where the bulk of the particles is located (i.e, nearvt). In this region, we show that the k-th moment of n y (t, x)/E(n y (t, x)) converges to the periodic function f k (x) identified in Section 3.2.
There have been several other works on different aspects of branching diffusions in periodic media, and the topic is closely related to reaction-diffusion equations with periodic coefficients. After presenting our results more precisely below, we discuss the relation to some of these other works in Section 3.4.
Acknowledgements: The work of James Nolen partially funded by grant DMS-1351653 from the US National Science Foundation and the work of Leonid Koralov and Pratima Hebbar was partially funded by grant W911NF1710419 from the Army Research Office.

Asymptotics of the transition kernel
Given a positive function h : R d → R that is sufficiently smooth, the h -transform of the operator L (given in (3)) is defined as for each real valued C 2 (R d ) function f . For each t ≥ 0 and x, y ∈ R d , the transition kernel u h (t, x, y) corresponding to L h satisfies: where u(t, x, y), satisfying (2), is the transition kernel corresponding to L (see Theorem 4.1.1 of [24]). We choose h from among a special family of eigenfunctions of L having exponential growth in a given direction. For ζ ∈ R d , let ϕ ζ be the principal positive periodic eigenfunction of the operator e −ζ·x L(e ζ·x ·). That is ϕ ζ satisfies with eigenvalue µ(ζ) ∈ R. Let ϕ * ζ denote the solution of the adjoint problem, that is, where µ * (ζ) is the principal eigenvalue of the adjoint operator, and hence µ * (ζ) = µ(ζ). We normalize ϕ ζ and ϕ * ζ by With this choice of h = h ζ , (4) can be written as Let us define p ζ (t, x, y) := e −tµ(ζ) u h ζ (t, x, y). The function p ζ (t, x, y) is the transition kernel for the operator Compared to L, this operators K ζ has an additional periodic drift a∇ log h ζ = aζ + a∇ log ϕ ζ , but no branching term r(x). Let ψ ζ and ψ * ζ denote the principal eigenfunctions corresponding to the principal eigenvalue (which is equal to zero) of the operator K ζ and K * ζ on the torus, respectively, and suppose that It is easy to see that ψ ζ (x) ≡ 1 and ψ * ζ (x) ≡ ϕ * ζ (x)ϕ ζ (x). Now we choose the direction ζ ∈ R d in an optimal way. Let Φ denote the Legendre transform of µ(ζ): The properties of µ, from Theorem 2.10 in Chapter 8 of the book of Pinsky [24], guarantee that Φ ∈ C 2 is well-defined. In particular Φ is strictly convex. For each c ∈ R d , the supremum in (10) is attained at a unique point which will be denoted byζ =ζ(c), that is Thus, c = ∇µ(ζ). In addition, for each c ∈ R d , we have ∇Φ(c) =ζ(c). Now, given Corresponding to this c, we choose the uniqueζ satisfying: Substituting ζ =ζ(c) in to (7), we obtain the identity Therefore, to obtain the exact asymptotics of u(t, x, y) in the domain of large deviations, we need to chooseζ appropriately, and provide an exact asymptotics of the transition density pζ(t, x, y). The reason for introducing this transformed kernel is that, momentarily assuming y = y(t) = x + ct, the effective drift of the process corresponding to pζ(t, x, y) is c. And therefore, the problem reduces to estimating the density of the transition kernel of the operator Kζ at a diffusive scale. The following proposition, which will be proved later, gives the exact asymptotics of the transition density pζ(t, x, y).
From Proposition 2.1, the following theorem now follows easily, giving the exact asymptotics of u(t, x, y). As we have mentioned, this result was announced in a talk of Agmon [1] in 2007: The following asymptotic relation holds as t → ∞ for all x, y ∈ R d such that y − x ≤ Lt: whereζ =ζ(t, x, y) = ∇Φ( y−x t ).
Proof of Theorem 2.2. Fix L > 0. From Proposition 2.1 and (13), for all (t, x, y) ∈ uniformly for y − x ≤ Lt. This concludes the proof of Theorem 2.2.
Since Φ is strictly convex, we definev ∈ R d to be the unique minimizer of Φ: We call thisv the effective drift of the branching diffusion process. The logarithmic asymptotics in Theorem 2.2 imply that a majority of the particles are located where The bounds (15) are valid at the large deviation scale, where y − x ≤ O(t). The following Aronson-type estimate provides a Gausian bound on the u that holds for all x, y ∈ R d , although it is less precise than (15). It is a consequence of Theorem 1.1 from Norris [23]: Lemma 2.3. Letv be the effective drift. There is a constant c > 0 such that Proof of Lemma 2.3. From (7) with ζ = 0, we have where p 0 (t, x, y) is the transition kernel for the operator K 0 in (9), having periodic coefficients, but without a potential term. The effective drift for p 0 is preciselyv = ℓ(0) = ∇µ(0). By Theorem 1.1 from Norris [23] there exists C > 0 such that for all x, y ∈ R d and t > 0, (See [25] Lemma 5.3 for an outline of the comparison of the setting in [23] to the setting here). Recall that µ(0) = −Φ(v). In terms of u, this implies that

Asymptotic behavior of a super-critical branching process in periodic media
In this section, we study the distribution of the number of particles in regions whose spatial location depends on time. Throughout this section, we will assume that the branching diffusion process is super-critical, meaning that wherev is the effective drift defined at (16). In view of Theorem 2.2, this condition implies that the total mass u(t, x, y) dy grows exponentially fast, as t → ∞.
Recall that n y (t, x) denotes the number of particles located in a unit d-dimensional cube containing y ∈ R d , assuming that, initially, there is one particle located at x ∈ R d . We state three theorems that describe different behaviors of the distribution of n y (t, x)/E(n y (t, x)). The main theorem in this section (Theorem 3.1) shows intermittency (i.e., the k-th moment dominates the k-th power of the first moment for some k), at locations with linear in time distances from the origin (recall that the bulk of the particles is located at the origin).

3.1.
Intermittency in the domain of large deviations. For y = (y 1 , y 2 , · · · , y d ), let Q d y denote the d-dimensional cube: . Recall that n y (t, x) denotes the number of particles located in Q d y , assuming that, initially, there is one particle located at x ∈ R d . (a) For each v ∈ R d , there exists the limit, Jensen's inequality implies that E(n tv (t, x) k ) ≥ (E(n tv (t, x))) k for each k ∈ N and v ∈ R d . Therefore, as long as the limits (19) This is the phenomenon of intermittency. This behavior is markedly different from the behavior in the case when the branching rate c(x) is compactly supported in space. In fact, for super-critical branching processes with compactly supported branching rates, in [18], it is shown that, that n tv (t, x) converges after appropriate scaling as t → ∞, and the quantities E(n tv (t, x) k )/E(n tv (t, x)) k converge to the corresponding quantities for the limiting random variable.
Remark 3.2. Formula (20) essentially provides a criterion for establishing weather intermittency occurs or not, in terms of a variational problem. To see this, we demonstrate the case k = 2. If (w, u) = (v(1 − u), u), and u ↑ 1, then, the term inside the supremum achieves the value 2γ 1 (v). This value of (w, u) lies on the boundary of the domain R d × (0, 1). Thus, intermittency would occur if there exists a different pair (w, u) ∈ R d ×(0, 1) such that the value of the supremum is greater that 2γ 1 (v). Otherwise, intermittency can not occur.

3.2.
Distribution of total number of particles. Following notation introduced in Section 2, recall that ϕ 0 is the principal periodic eigenfunction of the operator L. It satisfies with eigenvalue µ(0) ∈ R. The function ϕ * 0 will denote the solution of the adjoint eigenvalue problem: is the principal eigenvalue of the adjoint operator, and hence µ * (0) = µ(0). We normalize ϕ 0 and ϕ * 0 by In this section, to simplify notation, we will denote ϕ 0 , ϕ * 0 and µ(0) by ϕ, ϕ * and µ. For t > 0, x, y ∈ [0, 1) d , let ̺(t, x, y) denote the fundamental solution of the following PDE on the torus: ∂ t ̺(t, x, y) = L x ̺(t, x, y), ̺(0, x, y) = δ y (x). Observe that C 0 , ε > 0 such that, for every t > 0, (23) [0,1) d ̺(t, x, z)dz ≤ C 0 e tµ Let N(t, x) denote the total number of particles in R d at time t, assuming that, at time t = 0, there is one particle at x ∈ [0, 1) d . In the following theorem, all the moments of the normalized total number of particles are shown to converge. Theorem 3.3. For each k ∈ N, the following limit exists uniformly in x ∈ [0, 1) d : where the functions f k are defined recursively as follows, and, for k ≥ 2, where β k i = k!/(i!(k − i)!). In addition, there exists a real valued random variable ξ x , whose distribution is determined uniquely, such that E(ξ k x ) = f k (x) for each k ∈ N. The functions f k (x) defined recursively by the formulas (25) will be shown to be well defined, that is, the integrals in (25) will be shown to be convergent.
The above theorem implies that the total number of particles N(t, x), normalized by its expected value behaves "regularly". That is, the k-th moment of N(t, x) is commensurate with the k-th power of the first moment. In the next section, we show that n y (t, x) also exhibits the same "regular" behavior when y − tv = o(t). In contrast, in Section 3 we have shown that n tv (t, x) exhibits intermittent behavior when v =v, i.e., the k-th moment of n tv (t, x) grows much faster than the k-th power of the first moment for some k ∈ N.

3.3.
Distribution of the number of particles near the region where the bulk of the particles is located. We show that, at distances that grow sub-linearly in time from the bulk of the particles, all the moments converge. Let n y (t, x) denote the number of particles in Q d y at time t ∈ R + , given that there was one particle at x ∈ [0, 1) d at time t = 0. Define . From this formula for g(t, y), since the minimum of the twice continuously differentiable function Φ(v) is achieved at v =v, for each α ∈ (0, 1), we get (26) g(t, αy) ≥ g(t, y).

3.4.
Discussion. There have been several other works on different aspects of branching diffusions in periodic media, and the topic is closely related to reaction-diffusion equations with periodic coefficients. In particular, many authors have studied the spreading of wave fronts for reaction diffusion equations with periodic coefficents, having of the general form [13,10,28,4,5,6] and references therein. In one space dimension, the distribution of the maximal particle in the branching process, (the particle with largest spatial coordinate) can be expressed in terms of the solution to a reactiondiffusion equation of this KPP-type (see for example [21]), so that the asymptotic behavior of wave fronts as t → ∞ gives information about the behavior of the extremal particle in the branching process. A similar interpretation holds in the higher-dimensional setting.
When f (y, w) = g(y)w(1 − w) and g > 0 is strictly positive, a spreading phenomenon occurs: (See Chapter 7 of [10]). Hence, the set {tv ∈ R d | Φ(v) = 0} is understood as the asymptotic front of the wave as t → ∞. This front matches exactly the set t∂G 1 , where G 1 is defined in Theorem 3.1. The condition that g(y) > 0 is not necessary for such a spreading phenomenon. Berestycki, Hamel and Roques [5], [6] proved a necessary and sufficient condition for the spreading phenomenon (long-time survival of the branching process), which corresponds to the super-critical condition (18). They also analyze the effect of heterogeneity on the principle eigenvalue µ(0) of the associated linearized problem, and provide conditions under which the super-critical condition holds (see Theorem 2.12 of [5]). Refinements of the linear spreading rate have been obtained, even in the case of periodic media. For example, Hamel, Nolen, Roquejoffre, and Ryzhik [14] give sharper asymptotics for such fronts in periodic media in one space dimension, extending to the periodic case a well-known result of Bramson [7] which shows that the front (median of the extremal particle) moves as c 1 t−c 2 log(t)+O(1) as t → ∞. A key part of the analysis in [14] involves an estimate for a heat kernel analogous to p ζ (t, x, y) (for the transformed operator K ζ at (9)), except with Dirichlet boundary condition. This result was extended to fronts in multiple dimensions by Shabani [25]. Lubetzky, Thornett, and Zeitouni [20] have proved related asymptotics for the distribution of the extremal particle of a branching diffusion in periodic media. Unlike these works mentioned above, Theorem 3.1 pertains to the structure of the branching process behind the front, where the population is growing.

Proof of Proposition 2.1
Let X t be the diffusion process with generator K ζ (defined in (8)), From homogenization theory (see Freidlin [11] and the books of Bensoussan, Lions, and Papanicolaou [3] and of Jikov, Kozlov, Oleinik [16]), it is well known that the following result holds for diffusion processes with periodic coefficients: There exists a vector ℓ(ζ) ∈ R d (called the effective drift of X t ) and a positive definite matrix Ξ ζ (called the effective diffusivity of X t ) such that in distribution, where N (0, Ξ ζ ) denotes the normal random vector with mean zero and covariance matrix Ξ ζ . These quantities are given by the formulas: where η ζ (y) is a periodic (vector-valued) solution to which is determined uniquely up to an additive constant. These ℓ(ζ) and Ξ ζ are often called the effective drift and the effective diffusivity of the operator K ζ and hence, of the operator L h ζ since it only differs from K ζ by a constant potential term. For the operator L, notice that effective driftv, as defined at (16), corresponds tov = ℓ(0). We now state the following lemma about properties of the principal eigenvalue µ(ζ). The proof of this lemma can be found in the book of Pinsky [24] (Chapter 8, Theorem 2.10).
Lemma 4.1. The function µ : R d → R is twice continuously differentiable and strictly convex. In addition, for each ζ ∈ R d , and, Remark 4.2. From equation (11) and (12) we observe that corresponding to each (t, x, y) ∈ Since Φ is the Legendre transform of the function µ, we have the relation The proof of Proposition 2.1 is based on estimates of the local averages for z ∈ Z d and for appropriate choice of test functions f . We will choose f ∈ B, where B is the Banach space of Z d periodic continuous functions f : R d → C, equipped with the supremum norm. Observe that (36) has the form we define a family of measures on Z d : Before proving this, let us use this to finish the proof of Proposition 2.1.
Proof of Proposition 2.1. From Lemma 4.3 above, for the function g(k) = 1 0 (k), we have To prove Proposition 2.1, we would like to be able to replace f by a delta function at w ∈ [0, 1) d . This is easily justified if we have an appropriate bound on the derivative of p ζ (t, x, y) in the y variable. In this case, the weighted average of p ζ over a small domain approximates the value of p ζ at any point inside the domain. To get such bounds on the derivative of p ζ , we observe that is the fundamental solution of the PDE with periodic coefficients, with no potential term (see, for example, arguments in the proof of Lemma 2.3). From the Schauder estimate (see, Friedman [12]), it then follows that, . This is enough to conclude from that Note that the exponent in the above formula is slightly different. But the difference is negligible in the limit. Now suppose that L 0 > 0 is fixed, and y − x /t ≤ L 0 , for all x, y ∈ R d and t > 0. Then, recall from (34) that if we choose c = (y − x)/t, we have a correspondingζ such that ℓ(ζ) = c, or equivalently, ∇Φ(c) =ζ. Morevoer, there is L, depending on L 0 , such that |ζ| ≤ L holds if y − x /t ≤ L 0 . Thus, (41) can be applied to those c andζ uniformly to obtain lim We claim that for any L < ∞ fixed, the periodic eigenfunctions normalized by (6) satisfy Finally, we establish the claim (42). If this is not the case, then there must be sequences Since ζ and y are confined to a compact set, we can extract a subsequence of the pairs {(y n , ζ n )} ∞ n=1 that converges to some (ȳ,ζ). Joint continuity of (y, ζ) → ϕ ζ (y)ϕ * ζ (y) implies that ϕζ(ȳ)ϕ * ζ (ȳ) = 0, although the normalization (6) holds for ϕζ and ϕ * ζ . This is a contradiction, since the periodic principal eigenfunctions of elliptic operators L and L * have a strict sign. We conclude that (42) holds.
To complete the proof of Proposition 2.1, we now prove Lemma 4.3. This follows an argument of Hennion and Hervé [15] where a very similar lemma was proved (see Lemma VI.4 of [15]) in the discrete time one dimensional setting; we will explain the technical differences in Remark 4.6 below.
Let us define For g : Z d → R, for θ ∈ (S 1 ) d , z ∈ Z d , we use the following definitions of Fourier Transform and Inverse Fourier Transform: ). Now recalling the definition (37), observe that Using the Fourier inversion formula and Fubini's theorem, (45) can be written as The Fourier kernels {Q ζ (θ, t)} t≥0 are a family of compact operators on B and e iθℓ(ζ)t Q ζ (θ, t) is 2πZ d periodic in the parameter θ. One can show that for a fixed θ ∈ (S 1 ) d , the family {Q ζ (θ, ·)} t≥0 forms a semigroup. That is, for each x ∈ R d , t, s ≥ 0, Observe that, for θ = 0, Q ζ (0, t) is the Markov operator corresponding to the process X t , which is generated by K ζ . Therefore, since zero is the principal simple eigenvalue of the operator K ζ , 1 is the principal simple eigenvalue of the operator Q ζ (0, t). By a perturbation theorem (see, for example, Theorem III.8 in [15]), there exists a small θ 0 > 0 such that, for each θ ∈ (S 1 ) d with θ ≤ θ 0 , the principal eigenvalues of the operators Q ζ (θ, 1) are simple, for each ζ ≤ L. We denote these principle eigenvalues of the operators Q ζ (θ, 1) by λ(ζ, θ) ∈ C, for θ ≤ θ 0 . Thus, from the semigroup property (46) and the time homogeneity of the coefficients of the partial differential operator K ζ , we conclude that the principal eigenvalue of the operator The proof of Lemma 4.3 is based on the following spectral decomposition of the operator Q ζ (θ, t).
Proof of Lemma (4.4). In the discrete time one dimensional setting, Lemma (4.4) is proved in Hennion and Hervé [15] (see Proposition VI.2, therein), but the arguments there also go through in the continuous time d-dimensional setting. The assumptions of that Proposition, denoted by H ′′ [2] in [15] (assumptions on the Banach space being sufficiently big, Q ζ (0, 1) having 1 as its simple eigenvalue corresponding to the eigenfunction f ≡ 1, and the operators Q ζ (θ, 1) being sufficiently regular in the variable θ in a small neighborhood around θ = 0) are all satisfied in our setting, uniformly in ζ ≤ L. The proof of (49) (or, rather, its analog in [15]) relies on the fact that ∇ θ λ(ζ, θ) θ=0 = 0 and D 2 θ λ(ζ, θ) θ=0 = −Ξ ζ , which follows from arguments similar to those used in proving (35).
To apply Lemma 4.4 in the proof of Lemma 4.3, we will need the following fact about the eigenvalues λ(ζ, θ). For a bounded linear operator Q on Banach space B, let r(Q) denote its spectral radius.
Proof of Lemma 4.3. From Lemma 4.4, we know that there exists a θ 0 > 0 such that, for all θ ≤ θ 0 the decomposition (47) holds. Therefore, we can write m χ uniformly over χ ∈ (Z d , R d , B +,r , B 0 (L)). The change of variable θ = s √ t gives On the other hand, we havē where, Hence, We observe from (49) that the sequence {k t } t≥1 converges point-wise to k. Since the function g has bounded support in Z d , g ∞ < ∞. Thus, setting c g := g ∞ , we have By defining we now have β t ≤ α(θ 0 , L) t → 0 exponentially fast, as t → ∞. Now, where It is clear to see that lim t→∞ ǫ 4 t = 0, uniformly over χ ∈ (Z d , R d , B +,r , B 0 (L)). Combining these estimates, we conclude that Another difference is that, in our setting, the operators Q also vary with respect to the additional parameter ζ ∈ R d .

Proof of Theorem 3.1
The main idea of the proof is to look at the higher order correlation functions and the corresponding PDEs they solve and then use the asymptotics of the density function obtained in Theorem 2.2 and techniques developed in [18] to obtain logarithmic asymptotics of the moments E(n tv (t, x) k ).
Recall thatv = ℓ(0) = ∇µ(0) is the effective drift of the branching process defined at (16) (also see Lemma 4.1), and Φ(v) = −µ(0). Without loss of generality, we may assume thatv = 0, which simplifies our notation. Let B δ (y) denote a ball of radius δ > 0 centered at y ∈ R d . For t > 0 and x, y 1 , y 2 , ... ∈ R d with all y i distinct, define the particle density ρ 1 (t, x, y) and the higher order correlation functions ρ n (t, x, y 1 , ...., y n ) as the limits of probabilities of finding n distinct particles in B δ (y 1 ), ...B δ (y n ), respectively, divided by the n-th power of the volume of B δ (0) ⊂ R d . For a fixed y 1 , the density satisfies where L x is the linear operator defined at (3), acting on the variable x. The equations on ρ n , n > 1, are as follows (53) ∂ t ρ n (t, x, y 1 , y 2 , ..., y n ) = L x ρ n (t, x, y 1 , y 2 , ..., y n ) + α(x)H n (t, x, y 1 , y 2 , ..., y n ), ρ n (0, x, y 1 , y 2 , ..., y n ) ≡ 0, where H n (t, x, y 1 , y 2 , ..., y n ) = where Y = (y 1 , ..., y n ), U is a proper non-empty subsequence of Y , and |U| is the number of elements in this subsequence. See Section 2 of [18], for a derivation of these equations. Define m y k (t, x) = Q d y .... Q d y ρ k (t, x, y 1 , y 2 , ..., y k )dy 1 ...dy k . By integrating (53), it follows that The functions m tv i are related to the moments E(n tv (t, x) k ) according to where S(k, i) is the Stirling number of the second kind (the number of ways to partition k elements into i nonempty subsets). As explained in Section 9 of [18], this follows by partitioning Q d tv into small subdomains, and taking a limit as their diameters shrink uniformly to zero. 5.1. Proof of part (a). We first proof part (a) of Theorem 3.1. We will use induction to show the following: (i) For each k ≥ 1, there exists a constant a k > 0 such that (ii) For each k ≥ 1, for each L > 0, the following two limits exist uniformly for v ∈ R d , with v ≤ L and for x ∈ [0, 1) d , and satisfy when v ≤ L, k ≥ 2. In addition, γ k (v) ≥ γ k−1 (v) for k ≥ 2. Starting with k = 1, we estimate By Lemma 2.3, we know that there is c > 0, such that In view of (60) and the inequality − y − holds for some constant a 1 > 0. This proves (i) for k = 1. Now suppose that (i) holds up-to k − 1. From (55) and the Duhamel's Formula, we see that Note that, since (57) holds up-to k − 1, it also holds for k−1 i=1 β k i m y i (s, z)m y k−i (s, z) (with a different constantã k−1 ). Thus there exists a constant a k > 0 such that m y k (t, x) ≤ a k exp a k t − y−x 2 a k (t+1) , since the convolution of two functions satisfying the estimate (57), with two different constants also satisfies (57). That is, (i) holds for k, as well.
We next show that (ii) holds for k = 1 and k = 2, and (iii) holds for k = 2. Here is where we will need the sharp estimate for ρ 1 , provided by Theorem 2.2: for any fixed L > 0 and for all x, y ∈ R d with x − y ≤ Lt, we have where Φ, ϕ 0 and ϕ * 0 are defined before Theorem 2.2. From (64) and (60), we obtain and γ 1 is continuous since Φ is continuous. In addition, from (56), for each t > 0, E(n tv (t, x)) = m tv 1 (t, x). Thus (ii) holds for k = 1.
Next we show that, for k = 2, the first limit on the right hand side of (58) exists and satisfies formula (59). In the arguments below, we treat x and v as fixed, but all the estimates are easily seen to be uniform in v ≤ L and x ∈ [0, 1) d . Let us recall that We will apply Laplace's method to estimate the integral. For 0 < ε < κ < 1 and M > 0, consider the following partition of the domain [0, t] × R d : Then we write Using (62) in the region R 1 , where s < εt and vt − z > ε 1/4 t, we see that which can be made exponentially small (as t → ∞), with an arbitrarily large negative exponent, by choosing ε small enough. Therefore, using the estimate on ρ 1 from (61), we infer that for each r > 0, for all sufficiently small ε > 0, lim sup t→∞ ln I 1 (t, x, tv) t ≤ −r.
Similarly, considering the integral I 5 over the region R 5 , we may exchange the roles of (m tv 1 (s, z)) 2 and ρ 1 (t − s, z, x), to obtain for each r > 0, for all κ ∈ (0, 1) sufficiently close to 1, lim sup t→∞ ln I 5 (t, x, tv) t ≤ −r.
Let us now examine the asymptotics of I 4 = I 4 (t, x, tv), the integral over R 4 : Changing variables s/t = u and z/t = w, this is equivalent to The asymptotic behavior of m 1 and ρ 1 is available in (65) and Theorem 2.2. Observe that α(x) is periodic, non-negative and not identically 0. Therefore, following Laplace's method, we have as t → ∞, Combining this with the estimates on I 1 , I 2 , I 3 , I 5 , I 6 , we obtain that the first limit in (58) exists for k = 2, and is given by formula (59): From the formula above, since γ 1 (v) is continuous, we conclude that γ 2 is also continuous. Next we show that γ 2 (v) ≥ γ 1 (v) for all v ∈ R d . This is complete the proof that (iii) holds for k = 2. In view of (56), this also implies that, for k = 2, the second limit in (58) exists and is equal to γ 2 (v). Recall that m y 1 (t, x) and m y 2 (t, x) solve the following PDEs: x)) 2 , m y 2 (0, x) ≡ 0. We will show that there exists a C L > 0 such that, for each t ≥ 1 and x, y ∈ R d with x − y ≤ Lt, we have m y 2 (t, x) ≥ C L m y 1 (t, x). Fix R > 0 such that [0, 1) d ∈ B R (0). Observe that, since m y 1 (0, x) = χ Q d y (x), there exists a δ 1 > 0 such that m y Also observe that there exists a δ 2 > 0 such that, for all x, y ∈ R d with x − y ≤ 2R and t ∈ [1/4, 1/2], ρ 1 (t, x, y) ≥ δ 2 .
Now, observe that α(x) is periodic, non-negative and not identically 0. Thus, from (72), using Duhamel's Formula, for x ∈ Q d y , . Now, comparing the PDEs (71) and (72), and taking into account (73), we see that for all t ≥ 0, x, y ∈ R d , . For a fixed L > 0, for all x, y ∈ R d with x−y t ≤ L, t ≥ 1/2, from Theorem 2.2, there exists c > 0 such that (75) m y 1 (t, x) ≥ cm y 1 (t + 1/2, x). From (75) and (74), we conclude that there exists a constant C L > 0 such that (76) m y 2 (t, x) ≥ C L m y 1 (t, x), for all x, y ∈ R d with x−y t ≤ L and t ≥ 1. In particular, for each v ∈ R d , we have Because E(n tv (t, x) 2 ) is a linear combination of m tv 1 and m tv 2 (by (56)) this also implies Thus (ii) and (iii) hold for k = 2. This completes the basis for induction.
Next, suppose that (ii) and (iii) hold up to k − 1 with k ≥ 3: we will now show that (ii) and (iii) must also hold for k, completing the induction. From (63), there exists a constant C 1 > 0 such that Since E(n tv (t, x) k ) is a convex function of k, for each 1 ≤ i ≤ k − 1, Thus, using (56), there exists a constant C 2 > 0 such that, In order to prove that the first limit on the right hand side of (58) exists, we need to show that, We claim that, for all sufficiently large M > 0, As before, let 0 < ε < κ < 1, and partition the domain according to (66)-(67), and define the integrals so that I ℓ (t, x, tv) = 6 j=1 I ℓ j (t, x, tv). Using the same arguments as above, it is not difficult to show that, for each r > 0, for each δ > 0, for all sufficiently small ε > 0, for all κ ∈ (0, 1) sufficiently close to 1, for all sufficiently large M, Now consider integral I ℓ 4 (t, x, tv)t. Changing variables s/t = u and z/t = w, as before, this is equivalent to The logarithmic asymptotics of m 1 , m k−1 , and ρ 1 are given by (58) and Theorem 2.2. Therefore, following Laplace's method, as t → ∞, 1 t ln I ℓ 4 (t, x, tv)t is asymptotic to Combining these estimates, we conclude that Now, we justify (79), that is, the logarithmic asymptotics of the integrals I ℓ and I u are equal. The difference between I ℓ and I u is that E(n tv (t, x) i ) in I u replaces m tv i (t, x) in I ℓ . The properties of m tv i (t, x) that were used to derive the asymptotics of I 1 included estimate (57) and the uniform asymptotics of the logarithm (formula (58)). By the inductive assumption, the same uniform asymptotics holds for E(n tv (t, x) i ) for i ≤ k − 1. Moreover, by formula (56), the analogue of (57) holds for E(n tv (t, x) i ). That is, there exist constants d i > 0 such that for all (t, x, y) ∈ R + ×R d ×R d , for all 1 ≤ i ≤ k−1. Therefore, the logarithmic asymptotics of I 2 are the same as that of I 1 i.e., (79) holds. From (56), From the formula (80) which now holds for k and k − 1 and the inductive hypothesis that This, along with the inductive hypothesis that Therefore both the limits in (58) exist and are equal.
Form the inductive assumption that γ i is a continuous function for 1 ≤ i ≤ k − 1, using formula (82), we conclude that γ k is continuous. This concludes the proof of (i)-(iii) through induction.
We have shown that, for all sufficiently large M > 0, Therefore, letting M → ∞, we obtain the formula This completes the proof of (a) in Theorem 3.1.

Proof of part (b).
Using Hölder's inequality, it is easily seen that ln E(n tv (t, x) k ) is a convex function of k for each fixed t ∈ R + , v ∈ R d . In addition, γ 0 ≡ 0 and therefore γ k (v)/k is a non-decreasing function of k, which implies that, if γ k (v) > kγ 1 (v), then γ k+1 (v) > (k + 1)γ 1 (v). Therefore, G k+1 ⊆ G k must hold for each k ∈ N. We will complete the proof of Theorem 3.1 by showing that there exists a sequence of constants To justify this, we use induction. For k = 1, the statement is obvious since γ 1 achieves its maximum at 0. Now suppose the statement holds up to k − 1. Then, from the definition of γ k , We know that −Φ(0) = γ 1 (0) = µ(0) > 0, and Φ is continuous, therefore, the region G 1 is non-empty. Since the function v → γ k (v) is continuous for each k ≥ 1, the sets G k must be closed subsets of R d . Next let us show that each set G k contains a small ball centered at the origin. As a first step, the following lemma establishes an important property of the functions γ k .
Proof. We use induction for this proof. For k = 1, the statement of the lemma holds since γ 1 (v) is a twice differentiable strictly concave function andv = 0 is its maximizer.
Suppose the statement of the lemma holds for each 1 ≤ i ≤ k − 1. To show this for k, we have, for 0 ≤ α ≤ 1, Now, substituting w = αz, we have γ k (αv) = sup z∈R d ,u∈(0,1) Now, in order to prove that each set G k = {v ∈ R d γ k (v) = kγ 1 (v), γ 1 (v) ≥ 0} contains a small ball centered at the origin, we introduce functions f k defined below. For each k ≥ 2, we will first show that there is a small ball centered around the origin on which f k (v) = kγ 1 (v) ≥ 0. Then we will use induction to show that there is a (smaller) ball centered around the origin on which f k (v) = γ k (v).
Observe that f 2 = γ 2 . For k > 2, the formula for function f k is similar to the formula of γ k , but with (γ k−1 + γ 1 ) replaced by kγ 1 .
The analysis of g v k (w, u) is detailed in the following three lemmas. They show that, for each k ≥ 2, there is a small ball centered around the origin B β k (0), such that, for v ∈ B β k (0), the value of the supremum of g v k (w, u) on R d × (0, 1) is kγ 1 (v) which, as shown above, can be nearly achieved when w is close to 0 and u is close to 1.
Step II: Recall that Differentiating with respect to ε we obtain, Using the fact that the maximum of the function γ 1 (v) is achieved at v = 0 and the fact that γ 1 (v) is strictly concave, choose δ 2 ∈ (0, δ 1 ) be such that Let ε 2 ∈ (0, ε 1 ) be such that for each v ∈ B δ 2 2 (0), for all ε < ε 2 , and ℓ ≤ M, the vector v−ℓε 1−ε belongs to the B δ 2 (0). Now choose δ ∈ (0, δ 2 /2) such that, for each v ∈ B δ (0), we have for all ℓ ≤ M. This is possible since γ 1 achieves its maximum at 0, that is ∇γ 1 (0) = 0. Choose ε 0 > 0 with ε 0 < (0, ε 2 ) such that for all v ∈ B δ (0) and for all ℓ ≤ M, we have Thus, for all v ∈ B δ (0), for all ε < ε 0 and for all ℓ ≤ M, Thus, we conclude that, for all v ∈ B δ (0), But we know that, if (w, u) = (v(1 − u), u) and u approaches 1, the value of g v k (w, u) approaches kγ 1 (v). Therefore, The next lemma shows that the supremum of g cannot be achieved if u is close to 1, and w is separated from the origin.
Thus, for all v ∈ B δ (0), we have The last of the three lemmas, Lemma 5.4, shows that there is a small ball centered around the origin, on which the value of the supremum of g v k (w, u) in the region where w ∈ R d and u is away from 1 is strictly less than kγ 1 (v).
Lemma 5.4. For each ε > 0, there exists a δ > 0 such that, for all v ∈ B δ (0), Therefore, for all v ∈ B δ (0), Thus, by the above three lemmas, there exists a sequence of positive constants {β k } k≥1 such that, for all v ∈ B β k (0), Now let us show that there exists a sequence of positive constants {α k } k≥1 such that, for . This will be proved by induction. For k = 2, by the definition of γ 2 , we have that, . We need to show that there exists α k ∈ (0, β k ] such that, for all v ∈ B α k (0), we have To show this, it is enough to show that the supremum in the definition of γ k is achieved in the part of the space where the values of γ k−1 and (k − 1)γ 1 coincide. Let us define is dominated by the supremum of the same expression over the set Γ k (v). We will show that there exists α k > 0 such that, for all v ∈ B α k (0), we have Note that, for each (w, u) ∈ Γ k (v), the expression on the RHS is This expression, as follows from (87), is equal to (88) is justified by the following lemma.
Proof. The lemma will be proved in 2 steps. In Step I, the part of set Γ K (v) c where u is close to 1 is considered. In this part of the set, we make use of the fact that w is bounded from below. In Step II, the part of set Γ k (v) c where u is away from 1 is considered. In this part of the set, the left hand side of (89) can be made strictly smaller than kγ 1 (0), while the right hand side can be made arbitrarily close to kγ 1 (0) by choosing v in a small enough ball around the origin.

Proof of Theorem 3.3
Without loss of generality, we may assume thatv = 0, which simplifies our notation. Observe that the functions f k are clearly positive and continuous on the d-dimensional cube Q d 0 = [0, 1) d , from their recursive definition. As in (56), for each k ≥ 1, where ρ i 's are the particle density and higher order correlation functions, as defined in (52) and (53). Thus, we observe thatm i (t, x) satisfy the following PDEs on T d : . We will prove the following lemma after completing the proof of the theorem.
Using formula (95) in (92), we get Therefore, Now, we use induction to show that there exists a constant A > 0 such that, for every x ∈ [0, 1) d , f k (x) ≤ A k k!. For k = 1, we know that the eigenfunction ϕ(x) corresponding to the principle eigenvalue µ of the operator L on the d-dimensional torus T d is a positive and continuous function. Therefore, there exists a constant A 1 > 1 such that, for every Then, from the definition of the function f k , we get Recall that the operator L − µ has principle eigenvalue zero, while the principle eigenfunction of the adjoint operator (L − µ) * is ϕ * (with [0,1) d ϕ * (z) dz = 1). Therefore, there exists a constant C > 0 such that, for every x ∈ [0, 1) d , t > 0, Therefore, From the convergence of all the moments of N(t, x)/e µt , it follows that, there exists a random variable ξ x with the moments f k (x) (see [9]). The uniqueness of the distribution of ξ x follows from the bound on f k by the Carleman theorem. Except for a proof of Lemma 6.1, this concludes the proof of Theorem 3.3.
Proof of Lemma 6.1. We use induction to prove this lemma. The principle eigenvalue of the operator L is µ > 0, and the corresponding eigenfunction ϕ(x) > 0. Thus, from the theory of elliptic operators, from (93), there exists a function q 1 (t, x) such that where lim t→∞ q 1 (t, x) = 0 uniformly in x ∈ [0, 1) d . This gives (95) for k = 1 with f 1 (x) = ϕ(x). Suppose that the conclusion of the lemma holds up to k − 1, where k ≥ 2. From (94), using Duhamel's formula, we get By the inductive assumption, . After the change of variables u = t − s, we get Thus, we havem k (t, x) = e kµt f k (x) + q k (t, x) . It remains to show that lim t→∞ q k (t, x) = 0 uniformly in x ∈ [0, 1) d . Since the functions c, f i , f k−i are non-negative and continuous on [0, 1) d , there exists a constant C i > 0 such Therefore, from (23), the right hand side of the (96) goes to zero uniformly in x ∈ [0, 1) d .
To deal with the sum in the definition of q k (t, x), we break up the integral in two parts as follows, From (23), the integral in the first term is bounded and from the inductive hypothesis, and therefore, the first term converges to zero uniformly in x ∈ [0, 1) d . Similarly, from the inductive hypothesis, the supremum in the second term is bounded, while, from (23), the integral in the second term converges to zero uniformly in x ∈ [0, 1) d . Thus, we conclude that lim which completes the proof of Lemma 6.1.

Proof of Theorem 3.4
Without loss of generality, we may assume thatv = 0, which simplifies our notation. Recall from (56), We will show the following two statements by induction: (i) For each r(t) = o(t), there exists the limit lim t→∞ m y(t) uniformly in x ∈ [0, 1) d and y(t) ≤ r(t).
The theorem will then immediately follow from (i) since g(t, y) → ∞ as t → ∞ for y ≤ r(t) and therefore, the term with i = k dominates in the sum in formula (97).
Since the operator L is periodic, we first observe that m y k (t, x) = m . Thus, it is enough to show that there exists a C > 0 such that lim t→∞ sup g(s, y(t)) k s ∈ (εt, t)} g(t, y(t)) k t εt z−y(t) ≥r(t) α(z)ρ 1 (t − s, x, z)dzds ≤ C.
Choosing a sufficiently large L > 0, we use the asymptotic formula for ρ 1 (t, x, z) that was given in Theorem 2.2 in the region z − x ≤ Lt and the estimate (61) elsewhere, to obtain that (103) Note that e kt Φ y(t) t − s t Φ y(t) s e µ(t−s) = 1 when s = t. We show that, for sufficiently large t, the supremum in the above expression is achieved when s = t, when t is large enough.
Recall that Φ is continuous and the minimum value of the function Φ is achieved at 0, which is Φ(0) = −µ < 0. In addition, recall that r(t) = o(t). Thus, since k ≥ 2, we conclude that there exists η > 0 such that, for all sufficiently large t, kδΦ y(t) t + µδ < −η.
Thus, it is enough to show that, for all sufficiently large t, Indeed, for large t, the value of y(t)/t is close to 0, while y(t)/(t − δ) = y(t)/s ≤ 1 ε y(t)/t is also close to 0. Thus, using the fact that ∇Φ(0) = 0 and Taylor's formula, we obtain that (105) holds. Thus, we have shown that the supremum in (104) is achieved when s = t, when t is large enough. This completes the proof of (104). Next we show that lim t→∞ C k (t, x, y(t)) g(t, y(t)) k = f k (x). In the region z − y(t) ≤r(t), by the inductive assumption, we can replace Therefore, using a change of variable, it remains to show that, for each 1 ≤ i ≤ k − 1, Note that, g(t − s, y(t) − [z]) g(t, y(t)) = = Observe that ρ 1 is the fundamental solution of the operator L on R d , while ̺ is the fundamental solution of the same operator on T d . Therefore, for each s ≥ 0, x ∈ [0, 1) d , and each continuous Z d − periodic function h : R d → R, we have the relation for all z ≤ R, y(t) ≤ r(t) y(t) − z ≤r(t) and 0 ≤ s ≤ m. Therefore, we can choose δ > 0 small enough, such that, for all sufficiently large t, z−y(t) ≤r(t) z ≤R g(t − s, y(t) − [z])e µs g(t, y(t)) k − 1 α(z)f i (z)f k−i (z)ρ 1 (s, x, z)dz < η.
Following the same arguments that are detailed before (108), it is enough to show that, for all 1 ≤ i ≤ k − 1, uniformly in ȳ(t) ≥r(t), y(t) ≤ r(t). The idea here is that, ȳ as well as y(t) can be bounded from above by 2 z on the domain of integration. Therefore, repeating the arguments from (107), using Taylor's formula, given δ > 0 small. since z −ȳ(t) ≤ p(t) = o(t), y(t) ≤ r(t) = o(t) and ∇Φ(0) = 0, along with the estimate (61), there exists C > 0 such that, for all sufficiently large t and all 1 ≤ i ≤ k − 1, Since ȳ(t) ≥r(t) and z −ȳ(t) ≤p(t) = o(r(t)), we know, for sufficiently large t, z − x ≥r(t)/2. Thus, there exists a > 0 such that, for sufficiently large t, Therefore, there exists a constantC > 0 such that, for all 1 ≤ i ≤ k − 1, if δ is sufficiently small. This concludes the proof of (ii) for k.