# Random walk through a fertile site

Michel Bauer, P. L. Krapivsky, Kirone Mallick

► To cite this version:

Michel Bauer, P. L. Krapivsky, Kirone Mallick. Random walk through a fertile site. Physical Review E, 2021, 103 (2), pp.022114. 10.1103/PhysRevE.103.022114. hal-04014567

HAL Id: hal-04014567

https://hal.science/hal-04014567

Submitted on 20 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Random walk through a fertile site

Michel Bauer<sup>1,2</sup> P. L. Krapivsky<sup>3</sup> and Kirone Mallick<sup>1</sup>

<sup>1</sup>Institut de Physique Théorique, Université Paris-Saclay, CEA and CNRS, 91191 Gif-sur-Yvette, France

<sup>2</sup>Département de Mathématiques et Applications, ENS-Paris, 75005 Paris, France

<sup>3</sup>Department of Physics, Boston University, Boston, Massachusetts 02215, USA

<figure><img src="image_2.png" /></figure>

(Received 20 October 2020; accepted 15 January 2021; published 9 February 2021)

We study the dynamics of random walks hopping on homogeneous hypercubic lattices and multiplying at a fertile site. In one and two dimensions, the total number $N(t)$ of walkers grows exponentially at a Malthusian rate depending on the dimensionality and the multiplication rate $\mu$ at the fertile site. When $d > d_c = 2$, the number of walkers may remain finite forever for any $\mu$; it surely remains finite when $\mu \le \mu_d$. We determine $\mu_d$ and show that $\langle N(t) \rangle$ grows exponentially if $\mu > \mu_d$. The distribution of the total number of walkers remains broad when $d \le 2$, and also when $d > 2$ and $\mu > \mu_d$. We compute $\langle N^m \rangle$ explicitly for small $m$, and show how to determine higher moments. In the critical regime, $\langle N \rangle$ grows as $\sqrt{t}$ for $d = 3$, $t/\ln t$ for $d = 4$, and $t$ for $d > 4$. Higher moments grow anomalously, $\langle N^m \rangle \sim \langle N \rangle^{2m-1}$, in the critical regime; the growth is normal, $\langle N^m \rangle \sim \langle N \rangle^m$, in the exponential phase. The distribution of the number of walkers in the critical regime is asymptotically stationary and universal, viz., it is independent of the spatial dimension. Interactions between walkers may drastically change the behavior. For random walks with exclusion, if $d > 2$, there is again a critical multiplication rate, above which $\langle N(t) \rangle$ grows linearly (not exponentially) in time; when $d \le d_c = 2$, the leading behavior is independent on $\mu$ and $\langle N(t) \rangle$ exhibits a sublinear growth.

DOI: 10.1103/PhysRevE.103.022114

I. INTRODUCTION

We study noninteracting random walks (RWs) on the hypercubic lattice with a single fertile site where a RW may give birth to another RW. More precisely, we assume that each RW hops with the same unit rate to any neighboring site, so the overall hopping rate on the $d$-dimensional hypercubic lattice is $2d$. When a RW occupies the fertile site, the multiplication occurs at a rate $\mu$; the newborn RW appears at the same fertile site, it is fully matured, i.e., capable of reproducing right from the moment it was born. We assume that the process begins with a single RW at the fertile site; the extension to the general case when the initial number of RWs and their initial locations are arbitrary is rather straightforward.

This deceptively simple problem exhibits several counterintuitive behaviors, e.g., when the spatial dimension exceeds the lower critical dimension, $d > d_c = 2$, a phase transition occurs at a certain critical multiplication rate $\mu_d$. Many properties of critical behavior are universal (independent of the spatial dimension). Another important feature is the lack of self-averaging. Mathematically, it means that the distribution $P_N(t)$ of the total number $N(t)$ of RWs remains broad. This effect is particularly pronounced at the critical multiplication rate $\mu = \mu_d$ where the probability distribution $P_N(t)$ becomes asymptotically stationary in the $t \to \infty$ limit. We will see that this limiting distribution has a remarkably universal form, valid in any spatial dimension:

$$P_N(\infty) = \lim_{t \to \infty} P_N(t) = \frac{1}{\sqrt{4\pi}} \frac{\Gamma(N - \frac{1}{2})}{\Gamma(N + 1)}. \quad (1)$$

We now give a glimpse of our findings concerning average characteristics. The total number of RWs and the density at any site grow exponentially with time when $d = 1, 2$, and also in the supercritical regime, $\mu > \mu_d$, when $d > 2$. For instance,

$$N(t) \equiv \langle N(t) \rangle \sim e^{C_d t}. \quad (2)$$

The growth rate $C_d$ plays the role of the Malthusian parameter; we computed $C_d$ for hypercubic lattices $\mathbb{Z}^d$.

In the following we consider only hypercubic lattices $\mathbb{Z}^d$ and we choose the simplest initial state with a single RW at the fertile site. In this situation, we show that the threshold multiplication rate is given by

$$\mu_d = \frac{2}{W_d}, \quad (3)$$

where

$$W_d = \int_0^\infty dx [e^{-x} I_0(x)]^d \quad (4)$$

is the Watson integral [1]. For the critical multiplication rate, the average density at the fertile site satisfies

$$n_0(t) \approx \pi^{-2} \times \begin{cases} \mu_d^{-2} A_d^{-1} & d > 4 \\ \mu_4^{-2} [\ln t]^{-1} & d = 4, \\ \mu_3^{-2} (4\pi t)^{-1/2} & d = 3 \end{cases} \quad (5)$$
while the average number of RWs grows as

$$N(t) \approx \pi^{-2} \times \begin{cases} (A_d \mu_d)^{-1} t & d > 4 \\ \mu_4^{-1} \frac{t}{\ln t} & d = 4. \\ \mu_3^{-1} \sqrt{t/\pi} & d = 3 \end{cases} \quad (6)$$

The amplitudes $A_d$ are given in Eq. (66).

Thus above the lower critical dimension, $d > d_c = 2$, the exponential growth is possible only when $\mu > \mu_d$. The behavior at the critical multiplication rate, $\mu = \mu_d$, shows that the upper critical dimension $d^c = 4$ demarcates different growth laws. A qualitative explanation for the value $d^c = 4$ and the scalings in (6) is given in Appendix C 2.

The exponential growth above the lower critical dimension occurs only on average, the number of RWs may remain finite forever. More precisely

$$\mathrm{Prob}[\mathcal{N}(\infty) = \text{finite}] = \begin{cases} 0 & d \le 2 \\ \frac{\mu_d}{\mu} & d > 2, \quad \mu > \mu_d \\ 1 & d > 2, \quad \mu \le \mu_d \end{cases} \quad (7)$$

Therefore when $\mu > \mu_d$, the unlimited growth occurs with probability $1 - \mu_d/\mu$. When $\mu < \mu_d$, the number of RWs remains finite. For instance, the average eternal number of RWs is

$$\langle \mathcal{N}(\infty) \rangle = \frac{\mu_d}{\mu_d - \mu} \quad (8)$$

In physics literature, the problem has been examined in Refs. [2,3]. Our results are much more detailed; e.g., only average characteristics have been probed in those references. Several generalizations including biased RWs, systems with a few fertile sites, etc., have been additionally studied there. We do not treat these systems and merely remark that some extensions are rather straightforward. For instance, since linear equations govern the evolution of averages, the average characteristics in systems with many fertile sites can be deduced from the corresponding results with a single fertile site.

In mathematical literature, RWs on the lattice with branching at a single point have been also studied; see Refs. [4–8]. Death was included in most studies and in such situations the extinction is always feasible. A particular attention has been paid to the critical branching [9–12]. Random walks performing more complicated hopping have been also investigated (our RWs perform nearest-neighbor hopping); systems with a few fertile sites have been studied in Refs. [2,3,7,8,12]. Our results agree with previous findings whenever the models coincide. Initial conditions do not affect qualitative behaviors, so we consider the most natural initial condition with a single RW starting on the fertile site. In this setting, we obtain several explicit asymptotic behaviors, e.g., the long time behaviors of the probability distribution $P_N(t)$ in one and two dimensions are expressed through Catalan numbers. We briefly discuss a general situation applicable to arbitrary graphs and birth rates.

Our analysis as well as all previous studies [2–12] rely on the absence of interactions between random walkers. The extension to interacting many-particle systems is an important challenge. As a first step into this domain, we analyze the influence of the multiplication on a fertile site on the behavior of two interacting particle systems. One is the symmetric exclusion process in which RWs are subjected to the exclusion

constraint. In another example, there are no direct interactions but the birth is allowed only when the fertile site is occupied by a single particle.

The outline of this work is as follows. In Sec. II we study the one-dimensional model. Exact results in two dimensions are established in Sec. III. Explicit calculations become challenging in higher dimensions (Sec. IV), but we still derive a number of exact results like (1) and (8) and asymptotically exact results like (6)–(8). In Sec. V we study fluctuations, e.g., we compute the moments $\langle N^2 \rangle$ and $\langle N^3 \rangle$ for any $d$. The distribution of the number of RWs is studied in Sec. VI; we show that when $d > 2$, this distribution is asymptotically stationary in the critical and subcritical regimes, $\mu \le \mu_d$, and we determine it. In Sec. VII we show that when $d \le 2$, the region occupied by RWs grows ballistically with time and, apart from a few holes, this region is a segment in one dimension and a disk in two dimensions. In Sec. VIII we consider two interacting particle systems. In one system, particles interact through exclusion; in another, the birth is possible only when the fertile site hosts a single particle. In both examples, the growth is greatly suppressed compared to noninteracting RWs. A few technical calculations are relegated to Appendices A and B. In Appendix C we show how to adapt our approach to more general situations (arbitrary graphs, general birth rates, etc.), and we outline more mathematical techniques helpful in studying these generalizations.

## II. AVERAGE GROWTH IN ONE DIMENSION

In this section, we analyze the average growth in the one-dimensional lattice model. The governing equations for the densities are

$$\frac{dn_j}{dt} = n_{j-1} - 2n_j + n_{j+1} \quad (9a)$$

when $j \ne 0$. The density at the fertile site obeys

$$\frac{dn_0}{dt} = n_{-1} - 2n_0 + n_1 + \mu n_0. \quad (9b)$$

Making the Laplace transform with respect to time

$$\hat{n}_k(s) = \int_0^\infty dt e^{-st} n_k(t) \quad (10)$$

and the Fourier transform with respect to lattice sites

$$N(s, q) = \sum_{k=-\infty}^{\infty} \hat{n}_k(s) e^{-iqk}, \quad (11)$$

we recast (9a) and (9b) into

$$N(s, q) = \frac{1 + \mu \hat{n}_0(s)}{s + 2(1 - \cos q)} \quad (12)$$

Using

$$\hat{n}_0(s) = \int_0^{2\pi} \frac{dq}{2\pi} N(s, q)$$

and the identity

$$\int_0^{2\pi} \frac{dq}{2\pi} \frac{1}{s + 2(1 - \cos q)} = \frac{1}{\sqrt{s^2 + 4s}}, \quad (13)$$
we extract from (12) the Laplace transform of the density at the fertile site

$$ \hat{n}_0(s) = \frac{1}{\sqrt{s^2 + 4s - \mu}}. \qquad (14) $$

Thus

$$ N(s, q) = \frac{\sqrt{s^2 + 4s}}{\sqrt{s^2 + 4s - \mu}} \frac{1}{s + 2(1 - \cos q)}. \qquad (15) $$

The Laplace transform of the average number of RWs

$$ \hat{N}(s) = \sum_{k=-\infty}^{\infty} \hat{n}_k(s) = N(s, q = 0) \qquad (16) $$

is therefore given by

$$ \hat{N}(s) = \frac{1}{s} \frac{\sqrt{s^2 + 4s}}{\sqrt{s^2 + 4s - \mu}}. \qquad (17) $$

Inverting (14) we find the density at the fertile site

$$ n_0(t) = \int_{s_* - i\infty}^{s_* + i\infty} \frac{ds}{2\pi i} \frac{e^{ts}}{\sqrt{s^2 + 4s - \mu}}. \qquad (18) $$

An integration contour can go along any vertical line in the complex plane satisfying the requirement that $s_* = \text{Re}(s)$ is greater than the real part of singularities of the integrand. One can also deform a contour simplifying the extraction of asymptotic behavior. Instead, we rely on a useful general identity for inverse Laplace transforms. Suppose we know the inverse Laplace transform $f(t)$ of $\hat{f}(s)$. We actually want to determine the inverse Laplace transform of $\hat{f}(\sqrt{s^2 - a^2})$, and there is an expression through $f(t)$ which is valid for arbitrary $f(t)$. It reads [13]

$$ f(t) + a \int_{0}^{t} d\tau I_{1}(a\tau) f(\sqrt{t^{2} - \tau^{2}}), \qquad (19) $$

where $I_1$ is the Bessel function. Turning to (18) we notice that $\sqrt{s^2 + 4s} = \sqrt{(s + 2)^2 - 4}$ which coincides with $\sqrt{s^2 - a^2}$ if we choose $a = 2$ and make the shift $s \to s + 2$. Equation (18) implies $\hat{f} = 1/(s - \mu)$, the corresponding inverse Laplace transform is $f = e^{\mu t}$. Using these relations together with (19) we obtain

$$ e^{2t} n_0(t) = e^{\mu t} + 2 \int_0^t d\tau I_1(2\tau) e^{\mu \sqrt{t^2 - \tau^2}}. \qquad (20) $$

The integral in Eq. (20) looks simple, but apparently, it does not admit an expression in terms of standard special functions.

The second term on the right-hand side of (20) dominates in the long time limit. Thus

$$ e^{2t} n_0(t) \approx 2 \int_0^t d\tau I_1(2\tau) e^{\mu \sqrt{t^2 - \tau^2}}. \qquad (21) $$

Rescaling the time variable, $\tau = \eta t$, and using the well-known asymptotic

$$ I_1(x) \approx \frac{e^x}{\sqrt{2\pi x}} \quad \text{when } x \gg 1, \qquad (22) $$

we simplify (21) to

$$ e^{2t} n_0(t) \approx \sqrt{\frac{t}{\pi}} \int_0^1 \frac{d\eta}{\sqrt{\eta}} e^{\eta f(\eta)}, \qquad (23) $$

where $f(\eta) = \mu\sqrt{1 - \eta^2} + 2\eta$. The maximum of $f(\eta)$ is reached at $\eta_* = (1 + \mu^2/4)^{-1/2}$. Expanding $f(\eta)$ near $\eta_*$ and computing the Gaussian integral we arrive at a simple exponential asymptotic

$$ n_0(t) \approx A_1 e^{tC_1} \qquad (24) $$

with parameters

$$ C_1 = \sqrt{\mu^2 + 4} - 2, \quad A_1 = (1 + 4/\mu^2)^{-1/2}. \qquad (25) $$

Another way to establish (24) is to argue that the dominant contribution to the integral (18) is provided by an integral over a small contour surrounding the rightmost pole $s_+ = C_1$ of $1/[\sqrt{s^2 + 4s - \mu}]$; another pole is at $s_- = -\sqrt{\mu^2 + 4} - 2$. Near the pole $s_+ = C_1$ the singular part of the integrand in (18) is $A_1(s - s_+)^{-1}$; this leads to (24).

The average number of RWs is found from (17) to give

$$ N(t) = 1 + \mu \int_{s_* - i\infty}^{s_* + i\infty} \frac{ds}{2\pi i} \frac{e^{ts}}{s(\sqrt{s^2 + 4s - \mu})}. \qquad (26) $$

It is simpler to use the identity (valid for arbitrary lattice)

$$ N(t) = 1 + \mu \int_{0}^{t} d\tau n_{0}(\tau). \qquad (27) $$

Combining Eq. (27) with the asymptotic (24) we obtain

$$ N(t) \approx \frac{\mu^2}{\mu^2 + 4 - 2\sqrt{\mu^2 + 4}} e^{tC_1}. \qquad (28) $$

Specializing this result to small and large $\mu$ yields

$$ N(t) \approx \begin{cases} 2 e^{\mu^2 t/4} & \mu \ll 1 \\ e^{(\mu-2)t} & \mu \gg 1. \end{cases} \qquad (29) $$

When $\mu \ll 1$, the asymptotic behavior is the same as in the model with particles undergoing independent Brownian motions instead of random walks, and with birth happening at the origin and mathematically represented by $\mu\delta(x)$. When $\mu \gg 1$, the leading behavior is the same as in the zero-dimensional situation.

The density normalized by the density at the fertile site is asymptotically stationary. This remarkable property allows one to derive the asymptotic density profile in a simple manner. Plugging the ansatz $n_j(t) = n_0(t)m_j$ with time-independent $m_j$ into Eq. (9a) and using the asymptotic formula $\frac{dn_0}{dt} = C_1n_0$ we obtain the recurrence $m_{j-1} + m_{j+1} = (2 + C_1)m_j$ which has an exponential solution $m_j = \lambda^j$ with $\lambda$ satisfying $\lambda^{-1} - 2 + \lambda = C_1$. One root of this quadratic equation corresponds to $j > 0$, another to $j < 0$. Overall,

$$ \frac{n_j(t)}{n_0(t)} = \lambda^{|j|}, \quad \lambda = \frac{\sqrt{\mu^2 + 4} - \mu}{2}. \qquad (30) $$

Inserting $n_{\pm 1}(t) = \lambda n_0(t)$ into (9b) leads to the same result for $\lambda$; this provides a consistency check.

A more rigorous derivation of (30) that does not rely on factorization [i.e., on the ansatz $n_j(t) = n_0(t)m_j$] is given in Appendix A.
### III. AVERAGE GROWTH IN TWO DIMENSIONS

On the square lattice, the governing equations read

$$ \frac{dn_{i,j}}{dt} = \nabla^2 n_{i,j} + \mu n_{i,j} \delta_{i,0} \delta_{0,j}. \quad (31) $$

Here $\nabla^2$ is the discrete Laplacian defined via

$$ \nabla^2 n_{i,j} = n_{i-1,j} + n_{i,j-1} + n_{i+1,j} + n_{i,j+1} - 4n_{i,j} \quad (32) $$

for the square lattice.

Applying the Laplace-Fourier transform

$$ \hat{n}_{a,b}(s) = \int_0^\infty ds e^{-st} n_{a,b}(t), \quad (33a) $$

$$ N(s; p, q) = \sum_{a=-\infty}^{\infty} \sum_{b=-\infty}^{\infty} \hat{n}_{a,b}(s) e^{-i(pa+qb)} \quad (33b) $$

to (31) we obtain

$$ N(s; p, q) = \frac{1 + \mu \hat{n}_0(s)}{s + 4 - 2(\cos p + \cos q)}, \quad (34) $$

where $\mathbf{0} = (0, 0)$ is the fertile site. The definition (33b) allows us to express $\hat{n}_0(s)$ through the double integral

$$ \hat{n}_0(s) = \int_0^{2\pi} \frac{dp}{2\pi} \int_0^{2\pi} \frac{dq}{2\pi} N(s; p, q). $$

To compute the integral we use the identity

$$ \int_0^{2\pi} \frac{dp}{2\pi} \int_0^{2\pi} \frac{dq}{2\pi} \frac{1}{1 - z(\cos p + \cos q)/2} = \frac{2}{\pi} K(z), \quad (35) $$

where

$$ K(z) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1 - z^2 \sin^2 \theta}} \quad (36) $$

is the complete elliptic integral of the first kind. This allows us to fix $\hat{n}_0(s)$ and we arrive at

$$ \hat{n}_0(s) = \frac{1}{\Phi_2(s) - \mu}, \quad (37a) $$

$$ N(s; p, q) = \frac{\Phi_2(s) \hat{n}_0(s)}{s + 4 - 2(\cos p + \cos q)}, \quad (37b) $$

$$ \hat{N}(s) = \frac{1}{s} \frac{\Phi_2(s)}{\Phi_2(s) - \mu}, \quad (37c) $$

where we use the shorthand notation

$$ \Phi_2(s) = \frac{2\pi}{z K(z)} \quad \text{and} \quad z = \frac{4}{s + 4}. \quad (38) $$

Inverting (37a) we find the density at the fertile site

$$ n_0(t) = \int_{s_* - i\infty}^{s_* + i\infty} \frac{ds}{2\pi i} \frac{e^{is}}{\Phi_2(s) - \mu}. \quad (39) $$

In the long time limit, the leading contribution is again provided by an integral over a small circle surrounding the pole of $1/[\Phi_2(s) - \mu]$. Hence

$$ n_0(t) \approx \frac{1}{\Phi'_2(C_2)} e^{C_2 t} \quad (40) $$

<figure><img src="image_3.png" /><figcaption>FIG. 1. The Malthusian growth $C_2$ versus the multiplication rate $\mu$. Bottom curve: An implicit exact result given by Eq. (41). Top curve: An explicit approximate result, Eq. (44), formally valid when $\mu \ll 1$, but providing an excellent approximation up to $\mu < 1.8$.</figcaption></figure>

with $C_2$ determined from $\Phi(C_2) = \mu$. Using (38) we get

$$ kK(k) = \frac{2\pi}{\mu}, \quad k = \frac{4}{4 + C_2}. \quad (41) $$

The average number of RWs is asymptotically

$$ N(t) \approx \frac{\mu}{C_2 \Phi'_2(C_2)} e^{tC_2}. \quad (42) $$

This asymptotic behavior follows from (27) and (40).

When $\mu \gg 1$, the leading behavior is again essentially the same as in the zero-dimensional situation; an accurate analysis of (41) yields $N(t) \simeq e^{(\mu-4)t}$. In the opposite limit of the small multiplication rate, $\mu \ll 1$, the growth is still exponential, although the growth rate $C_2$ is extremely small. Using the asymptotic formula for the complete elliptic integral of the first kind

$$ K(k) = \ln \frac{4}{\sqrt{1-k^2}} + O\left[(1-k^2) \ln \frac{4}{\sqrt{1-k^2}}\right] \quad (43) $$

valid when $k \to 1 - 0$, one finds

$$ C_2 \approx 32 e^{-4\pi/\mu} \quad (44) $$

for $\mu \ll 1$. When $\mu < 1.8$, the explicit formula (44) provides an excellent approximation of the Malthusian growth rate $C_2$ given by the implicit relation (41); see Fig. 1.

The normalized density is asymptotically stationary. Inserting

$$ \frac{n_{i,j}(t)}{n_0(t)} = m_{i,j} \quad (45) $$

into (31) we obtain

$$ C_2 m_{i,j} = \nabla^2 m_{i,j} + \mu \delta_{i,0} \delta_{0,j} \quad (46) $$

with $\nabla^2$ defined in Eq. (32). Combining (46) with the generating function

$$ \mathcal{M}(u, v) = \sum_{i=-\infty}^{\infty} \sum_{j=-\infty}^{\infty} m_{i,j} u^i v^j, \quad (47) $$
we deduce

$$ \mathcal{M}(u, v) = \frac{\mu}{C_2 + 4 - u - u^{-1} - v - v^{-1}}. \quad (48) $$

We can write $m_{a,b}$ as a double contour integral, each over a unit circle in the complex plane:

$$ m_{a,b} = \frac{1}{(2\pi i)^2} \int_{|u|=1} \frac{du}{u^{1+a}} \int_{|v|=1} \frac{dv}{v^{1+b}} \mathcal{M}(u, v). \quad (49) $$

One can express the integrals through generalized hypergeometric functions; see Ref. [14]. We do not present those cumbersome results and just remark that $m_{0,\pm 1} = m_{\pm 1,0}$ can be deduced without computations. Indeed, using (46) and recalling that $m_{0,0} = 1$ we obtain

$$ m_{0,\pm 1} = m_{\pm 1,0} = 1 + \frac{C_2 - \mu}{4}. \quad (50) $$

## IV. AVERAGE GROWTH WHEN $d > 2$

### A. Green function approach

If $\mu = 0$, the initial value problem

$$ \frac{dn_i}{dt} = \nabla^2 n_i, \quad n_i(0) = \delta_{i,0} $$

is readily solved (see, e.g., Ref. [18]):

$$ n_i(t; \mu = 0) = \prod_{a=1}^{d} e^{-2t} I_{i_a}(2t). \quad (51) $$

Generally when $\mu > 0$, we take into account that the governing equations are linear and that additional particles are born with rate $\mu n_0(\tau)$ and express the density as a sum of (51) and a linear combination of similar solutions caused by particles that are born during the time interval $(0, t)$. In particular, at the origin

$$ n_0(t) = [J_0(t)]^d + \mu \int_0^t d\tau n_0(t-\tau)[J_0(\tau)]^d, \quad (52) $$

where we have used the shorthand notation

$$ J_0(t) \equiv e^{-2t} I_0(2t). \quad (53) $$

Therefore instead of the Laplace-Fourier transform, one can rely on the integral equation (52). Some results are easier to derive from the integral equation (52). For instance, the exponential growth $n_0(t) \sim e^{C_d t}$ is consistent with (52) when the growth rate $C_d$ satisfies

$$ 1 = \mu \int_{0}^{\infty} d\tau e^{-C_d \tau} [J_0(\tau)]^d. \quad (54) $$

Recalling the definition (4) of the Watson integral and taking into account that $C_d \ge 0$ we see that the right-hand side of (54) cannot exceed $\mu W_d/2$. Thus

$$ \mu W_d \ge 2. \quad (55) $$

The long time asymptotic $e^{-2t} I_0(2t) \simeq (4\pi t)^{-1/2}$ implies that the Watson integral (4) diverges when $d \le 2$. This together with the obvious fact that the right-hand side of (54) is a decreasing function of $C_d$ shows that (54) admits a single solution which is positive: $C_d > 0$. On the other hand, the Watson integral (4) converges when $d > 2$ and hence if (55)

is not obeyed, there is no solution to (54). This completes the derivation of the announced expression (3) for the critical multiplication strength $\mu_d$.

In three dimensions, the Watson integral (4) can be expressed [1] via Euler's gamma functions:

$$ W_3 = \frac{\sqrt{6}}{96\pi^3} \Gamma\left(\frac{1}{24}\right) \Gamma\left(\frac{5}{24}\right) \Gamma\left(\frac{7}{24}\right) \Gamma\left(\frac{11}{24}\right). \quad (56) $$

Numerically

$$ \begin{aligned} W_3 &= 0.505462020374\dots, \\ W_4 &= 0.309866780462\dots, \\ W_5 &= 0.231261630449\dots, \end{aligned} $$

etc. The asymptotic behavior $\lim_{d \to \infty} dW_d = 1$ follows from the definition (4) of the Watson integral. Hence the threshold value grows according to

$$ \mu_d \approx 2d \quad (57) $$

when $d \to \infty$.

### B. Laplace-Fourier transform

To probe the behavior in the $\mu \le \mu_d$ range, we apply the Laplace-Fourier transform:

$$ \hat{n}_a(s) = \int_0^\infty ds e^{-st} n_a(t), \quad (58a) $$

$$ N(s; \mathbf{q}) = \sum_{\mathbf{a}} \hat{n}_{\mathbf{a}}(s) e^{-i\mathbf{q} \cdot \mathbf{a}}, \quad (58b) $$

where $\mathbf{a} = (a_1, \dots, a_d)$, $\mathbf{q} = (q_1, \dots, q_d)$, and

$$ \begin{aligned} \mathbf{q} \cdot \mathbf{a} &= q_1 a_1 + \dots + q_d a_d, \\ \sum_{\mathbf{a}} &= \sum_{a_1=-\infty}^{\infty} \dots \sum_{a_d=-\infty}^{\infty}. \end{aligned} $$

We find

$$ N(s; \mathbf{q}) = \frac{1 + \mu \hat{n}_0(s)}{s + 2C(\mathbf{q})}, \quad C(\mathbf{q}) = d - \sum_{a=1}^{d} \cos q_a. \quad (59) $$

To avoid cluttered notation we write

$$ \frac{1}{\Phi_d(s)} = \int \frac{d\mathbf{q}}{s + 2C(\mathbf{q})}, \quad (60) $$

where

$$ \int d\mathbf{q} = \int_0^{2\pi} \frac{dq_1}{2\pi} \dots \int_0^{2\pi} \frac{dq_d}{2\pi}. $$

Fixing $\hat{n}_0(s)$ as before we arrive at

$$ \hat{n}_0(s) = \frac{1}{\Phi_d(s) - \mu}, \quad (61a) $$

$$ N(s; \mathbf{q}) = \frac{\Phi_d(s) \hat{n}_0(s)}{s + 2d - 2C(\mathbf{q})}, \quad (61b) $$

$$ \hat{N}(s) = \frac{1}{s} \frac{\Phi_d(s)}{\Phi_d(s) - \mu}. \quad (61c) $$

Note that

$$ \Phi_d(0) = \mu_d \quad (62) $$

remains positive when $d > 2$.
1. Supercritical regime: $\mu > \mu_d$

The growth is exponential, $n_0(t) \sim e^{s^d}$, with $s_d$ following from

$$ \Phi_d(s_d) = \mu. \qquad (63) $$

It is straightforward to verify that $s_d = C_d$, which was determined by (54).

2. Subcritical regime: $\mu < \mu_d$

In this range, the average number of RWs is finite:

$$ \lim_{t \to \infty} N(t) = N_d > 0. \qquad (64) $$

Hence $\lim_{s \to 0} s \hat{N}(s) = N_d$ which is consistent with (61c) when $N_d = \frac{\Phi_d(0)}{\Phi_d(0)-\mu}$. Using (62) we arrive at the general formula (8) for the average number of RWs.

3. Critical regime: $\mu = \mu_d$

Using the definition (60) one finds that

$$ \hat{n}_0(s) \approx \pi^{-2} \times \begin{cases} \mu_d^{-2} A_d^{-1} s^{-1} & d > 4 \\ \mu_4^{-2} [s \ln(1/s)]^{-1} & d = 4 \\ \mu_3^{-2} (4s)^{-1/2} & d = 3 \end{cases} \qquad (65) $$

when $s \to +0$ with

$$ A_d = (2\pi)^{-2} \int \frac{d\mathbf{q}}{[C(\mathbf{q})]^2}. \qquad (66) $$

Inverting (65) we arrive at the announced expressions (5) and (6) for the average density at the fertile site and the average number of RWs. We also establish the values of the amplitudes (66). The integral in (66) converges only when $d > 4$, and this explains why $d^c = 4$ plays the role of the upper critical dimension.

C. Density

The Laplace-Fourier transform of the density is exactly known; Eqs. (61a) and (61b) give

$$ N(s; \mathbf{q}) = \frac{\Phi_d(s)}{\Phi_d(s) - \mu} \frac{1}{s + 2d - 2C(\mathbf{q})}. \qquad (67) $$

Inverting this expression is tedious, and since we are mostly interested in the long time behavior it is convenient to rely on the already established asymptotic behavior of the density at the fertile site.

1. Supercritical regime: $\mu > \mu_d$

The normalized density is asymptotically stationary in this regime

$$ \frac{n_a(t)}{n_0(t)} = m_a. \qquad (68) $$

The generating function

$$ \mathcal{M}(\mathbf{u}) = \sum_{\mathbf{a} \in \mathbb{Z}^d} m_{\mathbf{a}} \mathbf{u}^{\mathbf{a}}, \quad \mathbf{u}^{\mathbf{a}} = \prod_{p=1}^{d} u_p^{a_p} \qquad (69) $$

is found as in two dimensions, and it is an obvious generalization of (48):

$$ \mathcal{M}(\mathbf{u}) = \frac{\mu}{C_d + 2d - \sum_{p=1}^{d} (u_p + 1/u_p)}. \qquad (70) $$

We give again the normalized density at sites neighboring the fertile site:

$$ m_{\pm 1, 0, \dots, 0} = 1 + \frac{C_d - \mu}{2d}. \qquad (71) $$

2. Critical regime: $\mu = \mu_d$

In the critical regime we use the same ansatz (68) and determine the generating function

$$ \mathcal{M}(\mathbf{u}) = \frac{\mu_d}{2d - \sum_{p=1}^{d} (u_p + 1/u_p)}. \qquad (72) $$

There are no simple general expressions for $m_a$ valid for all $\mathbf{a} \in \mathbb{Z}^d$. The normalized density at sites neighboring the fertile site admits a simple expression through the Watson integral

$$ m_{\pm 1, 0, \dots, 0} = 1 - \frac{\mu_d}{2d} = 1 - \frac{1}{dW_d}. \qquad (73) $$

Noting that $m_a$ satisfies a discrete Poisson equation

$$ \nabla^2 m_a + \mu_d \delta_0 = 0, \qquad (74) $$

we replace the discrete Laplacian by the continuous Laplacian far from the fertile site, $r = |\mathbf{a}| \gg 1$, and conclude that the solution approaches the Coulomb solution far away from the fertile site:

$$ m(r) \sim \frac{\mu_d}{r^{d-2}}. \qquad (75) $$

The average number of RWs is therefore

$$ N(t) \sim n_0(t) \int_0^{R(t)} dr r^{d-1} m(r) \sim n_0(t) R(t)^2. $$

Random walkers spread diffusively and the cutoff length is also expected to grow diffusively with time: $R(t) \sim \sqrt{t}$. Hence $N(t) \sim t n_0(t)$. This is consistent with (5) and (6).

V. FLUCTUATIONS

We are mostly interested in $N(t)$, the total number of RWs. Focusing on this global quantity allows us to suppress spatial aspects. The evolution of $N(t)$ can be interpreted as a branching process [15–17]. This change of view greatly helps in calculations.

A. Effective branching process

The mapping onto the branching process is simple. The primordial RW starting at the fertile site at $t = 0$ reproduces at a certain “branching” time $T_1$. The two RWs become the seeds of two independent branching processes. The branching times depend on the multiplication rate and on the first return probability to the fertile site and thus on the geometry of the lattice, but the overall procedure is universal (that is, valid for any lattice).

Let $P_N(t) = \text{Prob}[N(t) = N]$ be the probability distribution of the total number of RWs. The moment-generating
function

$$Z(\lambda, t) = \sum_{N=1}^{\infty} P_N(t) e^{\lambda N} \qquad (76)$$

satisfies an integral equation

$$Z(\lambda, t) = e^{\lambda} \Psi(t) + \int_{0}^{t} d\tau \, \psi(\tau) Z^{2}(\lambda, t - \tau). \qquad (77)$$

Here we shortly write

$$\Psi(t) = \mathrm{Prob}(T_1 = t), \qquad (78a)$$

$$\Psi(t) = \mathrm{Prob}(T_1 \ge t) = 1 - \int_0^t d\tau \, \psi(\tau). \qquad (78b)$$

Indeed, if there were no branching up to time $t$, we have $N = 1$ and $Z(\lambda, t) = e^{\lambda}$. This happens with probability $\Psi(t) = \mathrm{Prob}(T_1 \ge t)$ and results in the first term on the right-hand side of (77). The first branching may also occur at time $\tau$ in the range $\tau \in (0, t)$; this happens with probability density $\psi(\tau)$. There are then two independent processes with moment-generating functions $Z^{(1)}(\lambda, t - \tau)$ and $Z^{(2)}(\lambda, t - \tau)$. The total number of RWs is the sum $\mathcal{N} = \mathcal{N}^{(1)}(t - \tau) + \mathcal{N}^{(2)}(t - \tau)$. This leads to the product of the corresponding generating functions and results in the integral on the right-hand side of (77).

Thus, the problem reduces to solving a nonlinear integral equation (77). The probability density $\mathrm{Prob}(T_1 = \tau)$ encodes all the geometric data of the problem (the structure of the lattice and the spatial dimension).

## B. Zero-dimensional case

As a warm-up, we start with zero-dimensional case. In this situation

$$\Psi(t) = \mu e^{-\mu t}, \quad \Psi(t) = e^{-\mu t}, \qquad (79)$$

so the integral equation (77) becomes

$$Z(\lambda, t) = e^{\lambda - \mu t} + \mu \int_{0}^{t} d\tau \, e^{-\mu \tau} Z^{2}(\lambda, t - \tau). \qquad (80)$$

It is not immediately clear how to directly solve this integral equation. Fortunately, we can determine the moment-generating function since we know the distribution $P_N(t)$ in the zero-dimensional case. Indeed, the probabilities $P_N(t)$ satisfy exact rate equations

$$\frac{1}{\mu} \frac{dP_N}{dt} = (N-1)P_{N-1} - NP_N. \qquad (81)$$

Solving (81) subject to the initial condition $P_N(0) = \delta_{N,1}$ is straightforward (see, e.g., Ref. [18]). The solution reads

$$P_N(t) = e^{-\mu t} (1 - e^{-\mu t})^{N-1}. \qquad (82)$$

Hence the moment-generating function is

$$Z(\lambda, t) = \frac{1}{1 + e^{\mu t}(e^{-\lambda} - 1)} \qquad (83)$$

in the zero-dimensional case. One can verify that (83) is indeed the solution of (80) satisfying the initial condition

$$Z(\lambda, 0) = e^{\lambda}. \qquad (84)$$

## C. Perturbative expansion

Treating $\lambda$ as a small parameter we write

$$Z(\lambda, t) = 1 + \lambda \langle \mathcal{N} \rangle + \frac{\lambda^2}{2!} \langle \mathcal{N}^2 \rangle + \frac{\lambda^3}{3!} \langle \mathcal{N}^3 \rangle + \dots \qquad (85)$$

and plug this expansion into the governing equation (77).

Since $Z(0, t) = 1$ we find that (77) is satisfied at the zeroth order zero we recall (78b). Equation (77) is satisfied at the first order if the average number of particles $N(t) = \langle \mathcal{N}(t) \rangle$ obeys

$$N(t) = \Psi(t) + 2 \int_{0}^{t} d\tau \, \psi(\tau) N(t - \tau). \qquad (86)$$

This linear integral equation can be solved by the Laplace transform if we know the value of the branching probability. To ensure that (77) is satisfied at order 2 we must require that the second moment $M(t) = \langle \mathcal{N}^2 \rangle$ satisfies the same linear integral equation as (86), but with an extra source term

$$M(t) = \Psi(t) + 2 \int_{0}^{t} d\tau \, \psi(\tau) M(t - \tau) \\ + 2 \int_{0}^{t} d\tau \, \psi(\tau) N^{2}(t - \tau). \qquad (87)$$

Similarly the third moment $M_3(t) = \langle \mathcal{N}^3 \rangle$ satisfies

$$M_3(t) = \Psi(t) + 2 \int_0^t d\tau \, \psi(\tau) M_3(t - \tau) \\ + 6 \int_0^t d\tau \, \psi(\tau) N(t - \tau) M(t - \tau). \qquad (88)$$

Performing the Laplace transform of (86) we find

$$\hat{N}(s) = \frac{\hat{\Psi}(s)}{1 - 2\hat{\Psi}(s)}. \qquad (89)$$

The Laplace transform of (78b) gives

$$\hat{\Psi}(s) = \frac{1 - \hat{\psi}(s)}{s} \qquad (90)$$

and hence (89) simplifies to

$$\hat{N}(s) = \frac{1}{s} \frac{1 - \hat{\psi}(s)}{1 - 2\hat{\psi}(s)}. \qquad (91)$$

The consistency with (61c) allows us to fix

$$\hat{\psi}(s) = \frac{\mu}{\mu + \Phi_d(s)}. \qquad (92)$$

## D. Moments in one dimension

The asymptotic behavior of the second moment can be extracted directly from (87). First, we recall the already known asymptotic (28), which we rewrite as

$$N(t) \approx v_1 e^{C_1 t} \qquad (93)$$

with $v_1 = 1 + 2/\sqrt{\mu^2 + 4}$ and $C_1 = \sqrt{\mu^2 + 4} - 2$; see (25). The exponential growth (93) is consistent with (86) if

$$2 \int_{0}^{\infty} d\tau \, \psi(\tau) e^{-C_1 \tau} = 1. \qquad (94)$$

Equivalently, we rewrite (94) as $2\hat{\psi}(C_1) = 1$, and using (92) and $\Phi_1(s) = \sqrt{s^2 + 4s}$ we recover (25).
<figure><img src="image_4.png" /><figcaption>FIG. 2. The ratios $v_2/v_1^2$ (the bottom curve) and $v_3/v_1 v_2$ (the top curve) versus multiplication rate $\mu$. Both ratios are maximal, $\sqrt{8} + 2$ and $\sqrt{27} + 3$, respectively, when $\mu \to 0$; both ratios monotonically decrease and approach to 2 and 3, respectively, when $\mu \to \infty$.</figcaption></figure>

Inserting (93) into (87) we find $M \sim e^{2C_1 t}$ suggesting that we seek the solution in the form

$$M(t) \approx v_2 e^{2C_1 t}. \qquad (95)$$

Inserting (95) into (87) we find

$$v_2 = 2(v_2 + v_1^2) \int_0^\infty d\tau \, \psi(\tau) e^{-2C_1 \tau}. \qquad (96)$$

The integral in the above equation is $\hat{\psi}(2C_1)$, and using (92) we can express the ratio $v_2/v_1^2$ as

$$\frac{v_2}{v_1^2} = \frac{2\mu}{\Phi_1(2C_1) - \mu}. \qquad (97)$$

Using $\Phi_1(s) = \sqrt{s^2 + 4s}$ and (25) we obtain

$$\lim_{t \to \infty} \frac{M(t)}{[N(t)]^2} = \frac{v_2}{v_1^2} = \frac{2\mu}{2\sqrt{\mu^2 + 4 - 2\sqrt{\mu^2 + 4} - \mu}}. \qquad (98)$$

This ratio decreases from $\sqrt{8} + 2 = 4.828427...$ to 2 as $\mu$ increases from 0 to $\infty$; see Fig. 2. If the number of RWs were an asymptotically self-averaging quantity, the ratio would be equal to unity. Therefore $N(t)$ is a non-self-averaging quantity for all $\mu$.

The third moment grows according to

$$M_3(t) \approx v_3 e^{3C_1 t}, \qquad (99)$$

and after straightforward calculations one gets

$$\frac{v_3}{v_1 v_2} = \frac{6\mu}{\Phi_1(3C_1) - \mu}, \qquad (100)$$

which can be rewritten similarly to (98):

$$\frac{v_3}{v_1 v_2} = \frac{6\mu}{\sqrt{9\mu^2 + 48 - 24\sqrt{\mu^2 + 4} - \mu}}. \qquad (101)$$

The qualitative behavior of this ratio is similar to the behavior of the ratio (98), namely, it monotonically decreases from $\sqrt{27} + 3 = 8.19615242...$ to 3; see Fig. 2.

Generally the $n$th moment grows according to

$$M_n(t) \approx v_n e^{nC_1 t}, \qquad (102)$$

and the same calculations as above yield

$$v_n = \frac{\mu}{\Phi_1(nC_1) - \mu} \sum_{a=1}^{n-1} \binom{n}{a} v_a v_{n-a}, \qquad (103)$$

which should be solved for $n \ge 2$ with $v_1 = 1 + 2/\sqrt{\mu^2 + 4}$ playing the role of the initial condition. One can recursively determine any $v_n$. We haven't succeeded in solving (103) analytically, but some asymptotic behaviors can be deduced, see Appendix B.

### E. Moments in higher dimensions

The governing equations (86)–(88) are the same in any spatial dimension; the functions $\Psi(t)$ and $\psi(t)$ appearing in (86)–(88) depend on the dimensionality.

#### 1. Supercritical regime: $\mu > \mu_d$

In two dimensions, and also when $d > 2$ in the supercritical regime, the moments exhibit formally the same asymptotic behaviors as in one dimension:

$$N \approx v_1 e^{C_1 t}, \quad M \approx v_2 e^{2C_1 t}, \quad M_3 \approx v_3 e^{3C_1 t}, \qquad (104)$$

where

$$\Phi_d(C_d) = \mu, \quad v_1 = \frac{\mu}{C_d \Phi'_d(C_d)}. \qquad (105)$$

Equations (97), (100), and (103) remain applicable after the obvious replacement $C_1 \to C_d$ and $\Phi_1 \to \Phi_d$. For instance, the recurrence (103) becomes

$$v_n = \frac{\mu}{\Phi_d(nC_d) - \mu} \sum_{a=1}^{n-1} \binom{n}{a} v_a v_{n-a}.$$

These results are valid in one and two dimensions, and in the supercritical regime $\mu > \mu_d$ when $d > 2$.

#### 2. Critical regime: $\mu = \mu_d$

We perform the Laplace transform of (87) and find

$$\hat{M} = \frac{\hat{\Psi} + 2\hat{\psi}\hat{N}^2}{1 - 2\hat{\psi}}. \qquad (106)$$

The $s \to +0$ asymptotic behavior determines the large time asymptotic. One finds $\hat{\Psi} \ll \hat{N}^2$ when $s \to +0$, so

$$\hat{M} \approx \frac{2\hat{\psi}}{1 - 2\hat{\psi}} \hat{N}^2. \qquad (107)$$

Using additionally $\hat{\psi}(s) = \mu_d/[\mu_d + \Phi_d(s)]$, $\hat{\psi}(0) = 1/2$ and $\hat{n}_0(s) = 1/[\Phi_d(s) - \mu_d]$ we simplify (107) to

$$\hat{M} \approx 2\mu_d \hat{n}_0 \hat{N}^2 \qquad (108)$$

when $s \to +0$. Using (6) we compute

$$\hat{N}^2 \approx \begin{cases} 2\pi^{-4}(A_d \mu_d)^{-2} s^{-3} & d > 4 \\ 2\pi^{-4}\mu_4^{-2} s^{-3} [\ln(1/s)]^{-2} & d = 4. \\ \pi^{-5}\mu_3^{-2} s^{-2} & d = 3 \end{cases} \qquad (109)$$
We insert (65) and (109) into (108) to yield

$$ \hat{M} \approx \begin{cases} 4\pi^{-6}(A_d \mu_d)^{-3} s^{-4} & d > 4 \\ 4\pi^{-6}\mu_4^{-3} s^{-4} [\ln(1/s)]^{-3} & d = 4 \\ \pi^{-7}\mu_3^{-3} s^{-5/2} & d = 3 \end{cases} $$

from which

$$ M \approx \frac{2}{3} \times \begin{cases} \pi^{-6}(A_d \mu_d)^{-3} t^3 & d > 4 \\ \pi^{-6}\mu_4^{-3} t^3 [\ln t]^{-3} & d = 4. \\ 2\pi^{-15/2}\mu_3^{-3} t^{3/2} & d = 3 \end{cases} \quad (110) $$

In contrast to the behavior in the supercritical regime where $M \sim N^2$, we have $M \sim N^3$ in the critical regime. More precisely,

$$ \lim_{t \to \infty} \frac{M(t)}{[N(t)]^3} = \begin{cases} \frac{2}{3} & d \ge 4 \\ \frac{4}{3} & d = 3 \end{cases} \quad (111) $$

Therefore the behavior is strongly non-self-averaging in the critical regime.

Although the moments diverge as $t \to \infty$, the probability distribution $P_N(t)$ is asymptotically stationary:

$$ \Pi(N) = P_N(\infty) = \mathrm{Prob}[\mathcal{N}(\infty) = N]. \quad (112) $$

The divergence of the average, $\sum_{N \ge 1} N \Pi(N) = \infty$, is compatible with stationarity due to an algebraic tail of the distribution $\Pi(N)$. We derive the entire distribution later. Here we show how to establish the most interesting large N behavior relying only on consistency. We postulate $\Pi(N) \sim N^{-a}$ when $N \gg 1$ and note that $a < 2$ to agree with the divergence of the average. The lower bound, $a > 1$, ensures the normalization

$$ \sum_{N \ge 1} \Pi(N) = 1. \quad (113) $$

To match with the actual growth of the moments, we anticipate that the distribution $P_N(t)$ is stationary up to some growing crossover. Thus

$$ P_N(t) \approx \Pi(N) \sim N^{-a}, \quad 1 < a < 2 $$

when $1 \ll N \lesssim N_*$, and we additionally assume that the crossover grows algebraically with time, $N_* \sim t^\xi$; when $N \gg N_*$, the distribution $P_N(t)$ is nonstationary and it quickly vanishes. Using these assumptions one can determine the exponents $a$ and $\xi$. Indeed, we estimate two moments

$$ \langle \mathcal{N} \rangle \sim \sum_{N \ge 1} t^{\xi} \frac{N}{N^a} \sim t^{\xi(2-a)}, \quad \langle \mathcal{N}^2 \rangle \sim \sum_{N \ge 1} t^{\xi} \frac{N^2}{N^a} \sim t^{\xi(3-a)} $$

and use $M(t) \sim N(t)^3$ to get $\xi(3-a) = 3\xi(2-a)$ thereby fixing the exponent $a = 3/2$. Using asymptotic (6) we then fix the second exponent, viz., $\xi = 2$ when $d \ge 4$ and $\xi = 1$ when $d = 3$.

There is actually a logarithmic correction at the upper critical dimension $d^c = 4$ and the more precise expression for the crossover number of RWs is

$$ N_* \sim \begin{cases} t^2 & d > 4 \\ (t/\ln t)^2 & d = 4. \\ t & d = 3 \end{cases} \quad (114) $$

Thus we provided heuristic evidence for the tail

$$ \Pi(N) \sim N^{-3/2}. \quad (115) $$

One extra check of (115) is based on computing higher moments. Using (114) and (115) we find

$$ \langle \mathcal{N}^m \rangle \sim \begin{cases} t^{2m-1} & d > 4 \\ (t/\ln t)^{2m-1} & d = 4. \\ t^{m-1/2} & d = 3 \end{cases} \quad (116) $$

The same time dependence characterizes the critical behavior of the moments in the model of RWs with branching at the origin [4].

The calculation of $M_3 = \langle \mathcal{N}^3 \rangle$ can be done along the same lines as the calculation of $M = \langle \mathcal{N}^2 \rangle$ described above. Instead of (108) one finds

$$ \hat{M}_3 \approx 6\mu_d \hat{n}_0 \hat{N}M. \quad (117) $$

A long but straightforward calculation gives

$$ M_3 \approx \frac{4}{5} \times \begin{cases} \pi^{-10} (A_d \mu_d)^{-5} t^5 & d > 4 \\ \pi^{-10} \mu_4^{-5} (t/\ln t)^5 & d = 4 \\ \frac{16}{3} \pi^{-25/2} \mu_3^{-5} t^{5/2} & d = 3 \end{cases} \quad (118) $$

in agreement with (116). Similarly to (111) we have

$$ \lim_{t \to \infty} \frac{M_3(t)}{[N(t)]^5} = \begin{cases} \frac{4}{5} & d \ge 4 \\ \frac{64}{15} & d = 3 \end{cases} \quad (119) $$

This result and Eq. (111) quantify strongly non-self-averaging behavior in the critical regime.

Below we derive the critical stationary distribution (1), which rigorously confirms the tail (115).

### 3. Subcritical regime: $\mu < \mu_d$

In the subcritical regime, the probability distribution is stationary. We shall derive this stationary distribution below; see (139). Using this distribution one can compute any moment; e.g., the variance is given by

$$ \langle \mathcal{N}^2 \rangle - \langle \mathcal{N} \rangle^2 = \frac{2\mu_d\mu^2}{(\mu_d - \mu)^3}. \quad (120) $$

## VI. THE PROBABILITY DISTRIBUTION $P_N(t)$

### A. One dimension

In Sec. V D we have shown that in one dimension the moments $M_a(t) = \langle \mathcal{N}^a \rangle = \sum_{N \ge 1} N^a P_N(t)$ satisfy

$$ M_a(t) \sim [N(t)]^a, \quad (121) $$

where $N(t) = M_1(t)$. The growth laws (121) suggest that in the long time limit the probability distribution $P_N(t)$ acquires the scaling form

$$ P_N(t) = [N(t)]^{-1} \mathcal{P}(z), \quad z = \frac{N}{N(t)} \quad (122) $$

with $N(t)$ given by (28). More precisely, the scaling form (122) is expected to be valid in the limit

$$ N \to \infty, \quad t \to \infty, \quad z = \frac{N}{N(t)} = \text{finite}. \quad (123) $$
1. Small N behavior

In many problems, the scaling form remains applicable even when $N = O(1)$ and $t \to \infty$, but there are counterexamples, e.g., in submonolayer epitaxial growth [19]. In the present case, $P_N(t)$ also exhibits an unusual behavior for small $N$. Inserting (76) into (77) we obtain

$$P_1(t) = \Psi(t), \qquad (124a)$$

$$P_2(t) = \int_0^t d\tau \, \psi(t-\tau) P_1^2(\tau), \qquad (124b)$$

$$P_3(t) = 2 \int_0^t d\tau \, \psi(t-\tau) P_1(\tau) P_2(\tau), \qquad (124c)$$

$$P_4(t) = \int_0^t d\tau \, \psi(t-\tau) [2P_1(\tau)P_3(\tau) + P_2^2(\tau)], \qquad (124d)$$

etc. Using (90) and (92) with $\Phi_1 = \sqrt{s^2 + 4s}$ we get

$$\hat{P}_1(s) = \hat{\Psi}(s) = \frac{1}{s} \frac{\sqrt{s^2 + 4s}}{\mu + \sqrt{s^2 + 4s}} \qquad (125)$$

from which $\hat{P} \to (2/\mu)s^{-1/2}$ as $s \to +0$, implying that

$$P_1(t) \approx \sqrt{\frac{4}{\pi \mu^2 t}} \quad \text{as } t \to \infty. \qquad (126a)$$

Combining (124b) and (137) we deduce the asymptotic

$$P_2(t) \approx \frac{4}{\pi \mu^2 t} \quad \text{as } t \to \infty. \qquad (126b)$$

Continuing one deduces the asymptotic behavior

$$P_N(t) \approx \frac{1}{N} \binom{2N-2}{N-1} \left(\frac{4}{\pi \mu^2 t}\right)^{N/2} \qquad (126c)$$

with amplitudes being Catalan numbers.

The asymptotically exact results (126a)–(126c) seemingly disagree with the scaling form (122). However, when $N = O(1)$, the scaling variable $z$ in (122) vanishes as $t \to \infty$, while $z$ must be finite; see (123). Thus the scaling form is inapplicable when $N = O(1)$.

Let us estimate $N_*$ where the crossover to scaling form may occur. Catalan numbers grow as $4^N$, so

$$P_N \propto \left( \frac{64}{\pi \mu^2 t} \right)^{N/2} \qquad (127)$$

from (126c). The crossover from (127) to (122) apparently occurs when $t^{-N_*/2} \propto 1/N(t) \propto e^{-C_1t}$. Thus

$$N_* \sim \frac{C_1 t}{\ln t}. \qquad (128)$$

In the boundary layer, $N \ll N_*$, the distribution $P_N(t)$ varies according to (126c); the scaling apparently emerges when $N \gg N_*$. Similar behaviors with a boundary layer structure at small masses were found in models mimicking submonolayer epitaxial growth [19].

2. Large N behavior

To probe the large $z$ tail of the scaled distribution $\mathcal{P}(z)$ we use the identity

$$\int_0^\infty dz \, z^n \mathcal{P}(z) = \frac{v_n}{v_1^n}. \qquad (129)$$

The large $n$ behavior of the amplitudes $v_n$ is established in Appendix B. It gives

$$\int_0^\infty dz \, z^n \mathcal{P}(z) \approx \frac{\sqrt{\mu^2 + 4} - 2}{\mu} n! \left(\frac{\beta}{v_1}\right)^n \qquad (130)$$

and implies an exponential tail

$$\mathcal{P}(z) \approx \frac{\sqrt{\mu^2 + 4} - 2}{\mu} e^{-v_1 z / \beta} \quad \text{when } z \gg 1. \qquad (131)$$

We know $v_1 = 1 + 2/\sqrt{\mu^2 + 4}$, but $\beta = \beta(\mu)$ is unknown.

B. Two dimensions

In two dimensions, the scaling laws (121) hold and the probability distribution $P_N(t)$ is also expected to acquire the scaling form (122).

For small $N$ we again rely on Eqs. (124a)–(124d). Using (90) and (92) with $\Phi_2$ given by (38) we obtain

$$\hat{P}_1(s) = \hat{\Psi}(s) \approx \frac{4\pi}{\mu} \frac{1}{s \ln(32/s)} \qquad (132)$$

as $s \to +0$, implying that the probability for the primordial random walker still being alone at time $t \gg 1$ vanishes very slowly, viz., as the inverse logarithm:

$$P_1(t) \approx \frac{4\pi}{\mu} \frac{1}{\ln t}. \qquad (133a)$$

Combining (124b) and (146) we deduce the asymptotic

$$P_2(t) \approx \left( \frac{4\pi}{\mu} \frac{1}{\ln t} \right)^2 . \qquad (133b)$$

Similarly, using (146), (147), and (124c) we deduce

$$P_3(t) \approx 2 \left( \frac{4\pi}{\mu} \frac{1}{\ln t} \right)^3, \qquad (133c)$$

while from (146)–(148) and (124d) we obtain

$$P_4(t) \approx 5 \left( \frac{4\pi}{\mu} \frac{1}{\ln t} \right)^4 . \qquad (133d)$$

Computing the following asymptotic:

$$P_5(t) \approx 14 \left( \frac{4\pi}{\mu} \frac{1}{\ln t} \right)^5, \qquad (133e)$$

we recognize the pattern, and the amplitudes 1,1,2,5,14 remind us of the Catalan numbers. The general formula is

$$P_N(t) \approx \frac{1}{N} \binom{2N-2}{N-1} \left( \frac{4\pi}{\mu} \frac{1}{\ln t} \right)^N . \qquad (134)$$

The same argument as in the previous subsection shows that (134) is valid when $N \ll N_*$ with

$$N_* \sim \frac{C_2 t}{\ln[\ln t]} . \qquad (135)$$

The scaling form (122) emerges when $N \gg N_*$.
### C. Dimensions $d > 2$

When $\mu \le \mu_d$, the probability distribution $P_N(t)$ becomes asymptotically stationary, $P_N(\infty) = \Pi(N)$, in the long time limit. The moment-generating function is also asymptotically stationary, and $Y(\lambda) = Z(\lambda, t = \infty)$ satisfies a simple quadratic equation

$$Y = [1 - \hat{\psi}(0)]e^{\lambda} + \hat{\psi}(0)Y^2. \quad (136)$$

Recalling that $\hat{\psi}(0) = \mu/(\mu + \mu_d)$ and solving (136) we obtain

$$Y = \frac{\mu + \mu_d}{2\mu} \left\{ 1 - \sqrt{1 - \frac{4\mu\mu_d}{(\mu + \mu_d)^2} e^{\lambda}} \right\}, \quad (137)$$

which is expanded to yield

$$Y = \frac{\mu + \mu_d}{4\mu\sqrt{\pi}} \sum_{N \ge 1} \frac{\Gamma(N - \frac{1}{2}) (4\mu\mu_d)^N e^{\lambda N}}{\Gamma(N + 1) (\mu + \mu_d)^{2N}}. \quad (138)$$

Therefore

$$\Pi(N) = \frac{\mu + \mu_d}{4\mu\sqrt{\pi}} \frac{\Gamma(N - \frac{1}{2}) (4\mu\mu_d)^N}{\Gamma(N + 1) (\mu + \mu_d)^{2N}}. \quad (139)$$

This distribution has an exponentially decaying tail and an algebraically decaying $N^{-3/2}$ prefactor.

In the critical regime, $\mu = \mu_d$, Eq. (139) reduces to the announced formula (1). This remarkably universal result does not depend on the spatial dimension; the growth of the moments does depend on the dimensionality and also on more subtle properties of the lattice (we have considered only hypercubic lattices).

In the supercritical regime, $\mu > \mu_d$, the number of RWs may remain finite forever, although on average it grows exponentially. This suggests that in the long time limit the probability distribution $P_N(t)$ has a stationary part $\Pi(N)$ and an evolving part of the form (122). For $\lambda < 0$, the evolving part (122) does not contribute to the moment-generating function when $t \to \infty$ [since $N(t)$ becomes boundless], and the following asymptotically stationary expression is obtained:

$$Y(\lambda) = Z(\lambda, \infty) = \sum_{N=1}^{\infty} \Pi(N) e^{\lambda N}. \quad (140)$$

Thus (137)–(139) continue to hold in the supercritical regime. Equation (137) shows that the number of RWs remain finite forever with probability

$$\sum_{N=1}^{\infty} \Pi(N) = Y(0) = \frac{\mu_d}{\mu}. \quad (141)$$

With probability $1 - \mu_d/\mu$, the number of RWs diverges when $t \to \infty$. Hence we write $P_N(t)$ as a sum of the stationary distribution and an evolving scaling distribution

$$P_N(t) = \Pi(N) + [N(t)]^{-1} \mathcal{P}(z), \quad z = \frac{N}{N(t)}. \quad (142)$$

The scaled density satisfies

$$\int_0^\infty dz \mathcal{P}(z) = 1 - \frac{\mu_d}{\mu}, \quad (143a)$$

$$\int_0^\infty dz z^n \mathcal{P}(z) = \frac{\nu_n}{\nu_1^n}, \quad n \ge 1. \quad (143b)$$

<figure><img src="image_5.png" /><figcaption>FIG. 3. The rescaled front velocity $v/\mu$ versus the multiplication rate $\mu$.</figcaption></figure>

## VII. SPATIAL CHARACTERISTICS

The total number of RWs grows exponentially when $d=1$ and $d=2$. The region containing occupied sites,

$$\mathcal{D}(t) = \{j | n_j(t) > 0\}, \quad (144)$$

also tends to grow. The question is how. It is intuitively obvious that this region has a few holes, so it is essentially a droplet, that is effectively the region surrounded by the sea of empty sites. Let us disregard holes and determine the size and the shape of the droplet.

### A. One dimension

Denote by $r$ the rightmost occupied site: $n_r(t) > 0$ and $n_j(t) = 0$ for all $j > r$. The front position $r = r(t)$ is a random quantity. The leading behavior of this quantity is deterministic and can be determined using heuristic arguments. Equations (24) and (30) yield

$$n_j(t) = A_1 e^{C_1 t} \lambda^{|j|}, \quad \lambda = \frac{\sqrt{\mu^2 + 4} - \mu}{2} \quad (145)$$

with $A_1$ and $C_1$ given by (25).

The position of the front can be estimated from the criterion $n_r(t) \sim 1$, or the criterion

$$\sum_{j \ge r} n_j(t) \sim 1, \quad (146)$$

asserting that the total average number of RWs to the right of the front is of order one. Using (145) and any criterion we find that the front spreads ballistically

$$r(t) = vt \quad (147)$$

with velocity

$$v = \frac{C_1}{\ln(1/\lambda)} = \frac{\sqrt{\mu^2 + 4} - 2}{\ln[(\sqrt{\mu^2 + 4 + \mu})/2]}. \quad (148)$$

The ratio $v/\mu$ of the velocity to the multiplication rate exhibits the following limiting behaviors (see also Fig. 3)

$$\frac{v}{\mu} = \begin{cases} \frac{1}{2} - \frac{\mu^2}{96} + \cdots & \mu \to 0 \\ \frac{1}{\ln \mu} - \frac{2}{\mu \ln \mu} + \cdots & \mu \to \infty \end{cases}. \quad (149)$$

This ratio vanishes very slowly in the $\mu \to \infty$ limit.

We have used (145) for $j \sim t$, i.e., on distances greatly exceeding the diffusion scale, $j \sim \sqrt{t}$. The derivation of (145)
given at the end of Sec. II assumes the factorization property $n_j(t) = n_0(t)m_j$; a rigorous derivation is given in Appendix A. We have also ignored fluctuations which are substantial—the average value (25) of the amplitude $A_1$ in (145) is known, but different $A_1$ arise in different realizations. Fluctuations do not affect the leading behavior, however. Indeed, (146) gives

$$e^{C_1 t - r \ln(1/\lambda)} = \text{const}$$

with constant fluctuating from realization to realization. Thus a more accurate form of (147) is probably

$$r(t) = vt + \text{const} \qquad (150)$$

with constant fluctuating from realization to realization.

The droplet $D(t) = [\ell(t), r(t)]$, where $\ell(t)$ is the position of the leftmost particle, has a certain number of holes $H(t)$. It would be interesting to understand the statistics of this random quantity. It is not even clear whether it becomes stationary in the long time limit. Even if it does and the probability distribution $Q(h)$ is well defined, the moments may diverge.

## B. Two dimensions

Conjecturally, the droplet has a deterministic limiting shape as $t \to \infty$. More formally, this means that

$$\lim_{t \to \infty} t^{-1} D(t) = D_{\infty}. \qquad (151)$$

The normalized droplet is a disk, that is, the growth is asymptotically isotropic. The growth proceeds with a certain velocity $v$ which we determine below. The triviality of the limit shape is a nontrivial statement. Indeed, limit shapes often depend on the lattice and just a few are known even in two dimensions. As an example of the known limit shape different from the disk we mention an Ising droplet. This droplet is formed in the Ising ferromagnet on the square lattice endowed with zero-temperature spin-flip dynamics. More precisely, when a large domain of one phase is inside the sea of the opposite phase, the minority domain shrinks and approaches to the limit shape [20,21] different from the disk. The Eden-Richardson growth model [22–24] on the square lattice is among the known unknowns—the unknown limit shape is known to be different from the disk. The general rule is that if the growth is driven by the boundary like in the Eden-Richardson model, the lack of local isotropy results in a nontrivial limit shape. In our model, in contrast, most of the RWs are near the origin and the diffusion process is known to be asymptotically isotropic (see also Fig. 4).

The average normalized density is asymptotically stationary [see (45)], and it satisfies (46). Far away from the fertile site the governing equation (46) for the normalized density can be written in a continuous form

$$C_2 m = \nabla^2 m. \qquad (152)$$

The solution of this rotationally isotropic equation also enjoys rotational symmetry (far away from the fertile site). Thus we can rewrite (152) as

$$m'' + r^{-1} m' = C_2 m, \qquad (153)$$

<figure><img src="image_6.png" /><figcaption>FIG. 4. The evolution of the 2d droplet sampled at total population $10^n$ for $n = 4, \dots, 9$. The diffusion constant and multiplication rate are of the same order. A few faraway isolated walkers whose position is determined by diffusion alone are not shown. The picture is consistent with the linear in time growth of a disk.</figcaption></figure>

where prime denotes a derivative with respect to the radial coordinate $r = \sqrt{i^2 + j^2}$. The solution to (153) is

$$m(r) = CC_2 K_0(\sqrt{C_2} r). \qquad (154)$$

The numerical factor $C = O(1)$ remains undetermined in the realm of continuum framework. A linearly independent solution of Eq. (153) involving another modified Bessel function, $I_0(\sqrt{C_2} r)$, is absent in (154) since this solution diverges when $r \to \infty$. The criterion (146) gives

$$n_0(t) \int_{\sqrt{C_2} R}^{\infty} dx \, x K_0(x) \sim 1, \qquad (155)$$

where $R$ is the boundary of the droplet. By inserting the large time asymptotic $n_0(t) \sim e^{C_2 t}$ and

$$K_0(x) \approx \sqrt{\frac{\pi}{2x}} e^{-x} \quad \text{when } x \gg 1 \qquad (156)$$

into (155) we obtain

$$R(t) = vt, \quad v = \sqrt{C_2} \qquad (157)$$

in the leading order. The asymptotic behaviors are

$$v = \begin{cases} \sqrt{32} e^{-2\pi/\mu} & \mu \to 0 \\ \mu & \mu \to \infty \end{cases}. \qquad (158)$$

The top formula is asymptotically exact when $\mu \to 0$ but actually works very well up to $\mu < 1.8$.

In deriving (157) we used only the dominant exponential factor from (156). Taking into account an algebraic $x^{-1/2}$ prefactor and more carefully computing the integral in (155) we obtain

$$R(t) = vt + \frac{1}{2v} \ln(C_2 t), \quad v = \sqrt{C_2}. \qquad (159)$$
A logarithmic correction to the front position is known to occur (see Refs. [25–30] and references therein) in many traveling wave phenomena. In the present case, a logarithmic correction apparently arises only in two dimensions.

VIII. INTERACTING RANDOM WALKS

So far, we have investigated noninteracting particles performing identical RWs and multiplying at the fertile site. There are numerous interesting deformations of these simple dynamical rules. In this section, we discuss two simple deformations of the original model that include interactions.

A. Symmetric exclusion process

Here we consider a deformation of the original model based on including exclusion interaction. We assume that particles undergo identical nearest-neighbor symmetric hopping on $\mathbb{Z}^d$ and satisfy the constraint that each lattice site is occupied by at most one particle so that hopping to an occupied site is forbidden. This interacting particle system, known as the symmetric exclusion process (SEP), has achieved the status of a paradigm in statistical physics (see books and reviews [31–36]).

We should modify the birth rule as the newborn particle must be in a different site than the parent particle. One can postulate that the particle at the fertile site gives the birth, the newborn particle is put into a randomly chosen neighboring site of the fertile site, and the birth event is successful only if the chosen site is empty. Another birth rule is defined as follows: Whenever a particle at the fertile site hops to an (empty) neighboring site, it leaves the daughter particle at the fertile site with probability $p$. These two birth rules are essentially equivalent. Below we use the latter slightly simpler birth rule.

The major simplifying property of the SEP is that the density satisfies the diffusion equation [31–36], exactly as in the case of noninteracting RWs. This property is easy to appreciate. In one dimension, for instance, one writes an exact equation

$$ \frac{d\langle\tau_j\rangle}{dt} = \langle\tau_{j-1}(1-\tau_j)\rangle + \langle\tau_{j+1}(1-\tau_j)\rangle - \langle\tau_j(1-\tau_{j-1})\rangle - \langle\tau_j(1-\tau_{j+1})\rangle, \quad (160) $$

following from the rules of the SEP for all $j \neq 0$. (Hereinafter we use occupation numbers: $\tau_j = 1$ if site $j$ is occupied and $\tau_j = 0$ otherwise.) Massaging Eq. (160) one notices that second-order correlation functions like $\langle\tau_{j-1}\tau_j\rangle$ cancel. Therefore Eq. (160) simplifies indeed to the lattice diffusion equation

$$ \frac{dn_j}{dt} = \nabla^2 n_j \quad (161) $$

for the densities $n_j = \langle\tau_j\rangle$ when $j \neq 0$. The same equation describes the evolution in arbitrary dimension. The diffusion coefficient is the same as for random walkers, our convention being that the hopping rates to neighboring sites are equal to unity.

At the fertile site we have (again for concreteness in one dimension)

$$ \frac{d\langle\tau_0\rangle}{dt} = \langle\tau_{-1}(1-\tau_0)\rangle + \langle\tau_1(1-\tau_0)\rangle - (1-p)\langle\tau_0(1-\tau_{-1})\rangle - (1-p)\langle\tau_0(1-\tau_1)\rangle, \quad (162) $$

which becomes

$$ \frac{dn_0}{dt} = n_1 + n_{-1} - 2(1-p)n_0 - p\langle\tau_0(\tau_1 + \tau_{-1})\rangle. \quad (163) $$

Thus the evolution of the density $n_0$ at the fertile is coupled to the second-order correlation functions $\langle\tau_0\tau_1\rangle$ and $\langle\tau_0\tau_{-1}\rangle$. Exact equations for these correlation functions involve third-order correlation functions. This attempt to get a closed system of equations never ends leading to an infinite hierarchy.

Let us first consider the extreme case of $p = 1$. This case is tractable because the fertile site is always occupied. Therefore we do not need (163), we merely have the boundary condition

$$ n_0(t) = 1 \quad (164) $$

at the fertile site for all $t > 0$. The initial condition is

$$ n_j(0) = 0. \quad (165) $$

Thus in the extreme case we need to solve Eq. (161) subject to (164) and (165). A mathematically identical problem arises in various contexts, e.g., it governs the evolution of the two-body correlation function for the voter model and one simple catalysis problem [37–39]; it also obviously describes the SEP with an infinitely strong localized source [40–42]. Several exact and asymptotically exact behaviors are known. For instance, the average $N(t) \equiv \langle\mathcal{N}(t)\rangle$ exhibits an asymptotic growth [41]

$$ N(t) \approx \begin{cases} 4\sqrt{t/\pi} & d=1 \\ 4\pi t / \ln t & d=2. \\ \frac{d}{W_d} t & d>2 \end{cases} \quad (166) $$

The growth becomes linear in time above the critical dimension, $d > d_c = 2$. The amplitude of this linear growth involves the Watson integral $W_d$ which appeared in some previous formulas, e.g., in Eqs. (4) and (56).

Even in the extreme case the fluctuations of the random quantity $\mathcal{N}$ are essentially unknown, the only result known so far is the variance of $\mathcal{N}$ in one dimension [40,41]: The ratio of the variance to the average is asymptotically

$$ \lim_{t \to \infty} \frac{\langle\mathcal{N}^2(t)\rangle - \langle\mathcal{N}(t)\rangle^2}{\langle\mathcal{N}(t)\rangle} = 3 - \sqrt{8}. \quad (167) $$

Let us look at nonextreme versions of the model parameterized by $p \in (0, 1)$. The crucial feature of the SEP is the absence of correlations in equilibrium. Hence when $i \neq j$, we have $\langle\tau_i\tau_j\rangle = \langle\tau_i\rangle\langle\tau_j\rangle = n_i n_j$ in equilibrium. This is inapplicable in systems with flux, so we cannot, e.g., replace $\langle\tau_0\tau_1\rangle$ in (163) by $n_0 n_1$.

In one dimension, the flux vanishes. More precisely, the flux decays as $t^{-1/2}$ in the long time limit. This follows from Eq. (166) in the extreme case, $p = 1$, and clearly occurs for
all $0 < p \le 1$. Since the flux asymptotically vanishes, the behavior approaches to the behavior of the SEP at equilibrium when the correlators factorize [31–35]. Thus $\langle \tau_0 \tau_1 \rangle = n_0 n_1$ is asymptotically exact, so Eq. (163) simplifies to

$$ \frac{dn_0}{dt} = n_1 + n_{-1} - 2(1-p)n_0 - pn_0(n_1 + n_{-1}). \quad (168) $$

Summing (168) and all Eqs. (161) for $j \ne 0$, and taking into account the $n_j = n_{-j}$ symmetry, we obtain

$$ \frac{dN}{dt} = 2pn_0(1-n_1) \quad (169) $$

for $t \gg 1$.

We proceed on the “physical” level of rigor by making plausible guesses and checking consistency. The starting point is the asymptotic behavior

$$ 1 - n_j \approx \frac{A_j}{\sqrt{\pi t}} \quad (170) $$

valid when $t \gg 1$ and $j \ll \sqrt{t}$. In the extreme case, an exact expression for the density profile valid for all $j \ge 0$ and $t \ge 0$ is known [41]:

$$ n_j = e^{-2t} I_j(2t) + 2e^{-2t} \sum_{k>j} I_k(2t). \quad (171) $$

Using (171) one confirms the asymptotic (170) and gets $A_j = j$ in the extreme case.

Generally for arbitrary $p$ we insert (170) into Eqs. (161) and find $A_{j+1} - 2A_j + A_{j-1} = 0$, from which

$$ A_j = A_0 + j. \quad (172) $$

The amplitude in front of the linear term is fixed by the known asymptotic, $A_j \approx j$ when $j \gg 1$, which can be established by using a continuum approach valid when $j \gg 1$. Substituting (170) into (168) we deduce a relation $A_0 = (1-p)A_1$. Combining this result with (172) specialized to $j = 1$ we obtain $A_0 = (1-p)/p$ and $A_1 = 1/p$.

We can now derive the leading behavior of $N(t)$. Substituting $1 - n_1 \approx A_1/\sqrt{\pi t}$ into Eq. (169) and integrating we find

$$ N \approx 4pA_1 \sqrt{\frac{t}{\pi}} = 4\sqrt{\frac{t}{\pi}}. \quad (173) $$

The leading asymptotic growth is therefore independent on the birth probability $p \in (0, 1]$. This intriguing feature seems present only in the leading behavior. To confirm, or disprove, this assertion one would like to compute subleading terms. In the extreme case, we know the exact answer [41]

$$ N = e^{-2t}[I_0(2t) + 4tI_0(2t) + 4tI_1(2t)] \quad (174) $$

from which one can obtain the entire expansion

$$ \langle N \rangle = 4\sqrt{\frac{t}{\pi}} + \frac{1}{4}\frac{1}{\sqrt{\pi t}} + O(t^{-3/2}). $$

Generally when $p < 1$ we anticipate

$$ \langle N \rangle = 4\sqrt{\frac{t}{\pi}} + C_1(p) + O(t^{-1/2}) \quad (175) $$

with $C_1(p) < 0$ when $p < 1$. A similar expansion has been derived in Ref. [41] for the SEP with a source of finite strength, and it is probably valid in the present model when $0 < p < 1$.

In two dimensions, the flux also vanishes in the long time limit. Generalizing the above arguments one finds that the leading asymptotic remains the same as in the extreme model. The subleading term, however, is only logarithmically smaller than the leading term, and it probably depends on $p$. In other words, we anticipate

$$ \langle N \rangle = \frac{4\pi t}{\ln t} + \frac{C_2(p)t}{(\ln t)^2} + \dots \quad (176) $$

with large subleading correction, so the convergence to the leading asymptotic is extremely slow. The subleading correction in (176) is conjectural in the general case of $0 < p < 1$, but for $p = 1$ such correction and the exact expression for $C_2(1)$ was established in Ref. [42].

Thus when $d \le d_c = 2$, the density at the fertile site approaches to unity and the flux vanishes in the long time limit. This happens for all $p \in (0, 1]$. In contrast, the birth probability $p$ affects the leading behavior when $d > 2$. Similarly to noninteracting RWs, we anticipate different behaviors depending on whether the birth probability $p$ is smaller, equal, or larger than the critical birth probability $p_c(d)$. For any $p < 1$, the total number of particles may remain finite when $d > 2$. Furthermore, the total number of particles will remain surely finite for sufficiently small $p$. This feature makes plausible the existence of the critical value such that the total number of particles is surely finite when $p \le p_c(d)$. Thus for $p \le p_c(d)$ the particle number distribution is expected to be asymptotically stationary. The form of this stationary distribution is unknown.

In the supercritical regime, $p > p_d$, we anticipate the same linear in time growth as in the extreme case:

$$ N(t) \approx S_d(p)t \quad (177) $$

when $d > 2$ and $p > p_c(d)$. We know that the amplitude $S_d(p)$ is a strictly increasing function of the birth probability on the interval $p_c(d) < p < 1$. We also know that $S_d(p) = 0$ when $p = p_c(d)$ and $S_d(1) = d/W_d$.

Recall that for RWs in the critical regime, the density at the fertile site vanishes as $t \to \infty$ when $d = 3$ and $d = 4$, and remains finite when $d > 4$; see (5). The SEP is essentially identical to noninteracting RWs at small density, and hence at least when $d = 3$ and $d = 4$ we anticipate the same qualitative behaviors in the critical regime as for RWs. Thus when $p = p_c(d)$, the density at the fertile site is expected to decay as

$$ n_0(t) \sim \begin{cases} t^{-1/2} & d = 3 \\ [\ln t]^{-1} & d = 4 \end{cases} \quad (178) $$

while the average number of particles is expected to grow according to

$$ N(t) \sim \begin{cases} \sqrt{t} & d = 3 \\ t/\ln t & d = 4 \end{cases} \quad (179) $$

These asymptotic behaviors are consistent with the tail $\Pi(N) \sim N^{-3/2}$ in the critical regime, the same tail as for RWs in the critical regime; see (115). Hence the higher moments
<figure><img src="image_7.png" /><figcaption>FIG. 5. The model with quiet birth in the situation when multiplication may occur, i.e., the fertile site is occupied by a single RW. The birth event may occur at the moment when the particle leaves the fertile site. The random walkers are shown as squares, the particle at the fertile site, j = 0, has a disk indicating that it may multiply at the moment it leaves the fertile site. (a) In the illustration, the occupation numbers are n<sub>-1</sub> = 3, n<sub>-2</sub> = 4, n<sub>-3</sub> = 4, . . . on the left of the fertile site where n<sub>0</sub> = 1, and n<sub>1</sub> = 5, n<sub>2</sub> = 7, n<sub>3</sub> = 4, . . . on the right of the fertile site. (b) The random walker has jumped from the fertile site to the right, and the multiplication event has occurred; the multiplication certainly occurs in the extreme version, p = 1. Immediately after the jump, the fertile site is occupied by a single particle, and hence the multiplication is feasible. If another random walker jumps at the fertile site, the multiplication would be temporarily impossible.</figcaption></figure>

⟨N<sup>m</sup>⟩ probably exhibit the same dynamical behaviors (116) as in the case of RWs.

The behaviors (177)–(179) are conjectural. The rates S<sub>d</sub>(p) in (177), the amplitudes in (178) and (179), and the critical birth probabilities p<sub>c</sub>(d) are unknown. For RWs in the critical regime, the final particle number distribution RWs is universal (independent of the spatial dimension). The derivation of that property relied on the strict absence of interactions between RWs. Furthermore, this property concerns Π(N), not its asymptotic behavior, so there is no ground for any guess about Π(N) in the case of the SEP in the critical regime.

The region of occupied sites at a given time is not a droplet, there are numerous holes in the case of the SEP. If, however, we consider the domain of sites *visited* during the time interval (0, t), this domain is asymptotically a growing ball. In the extreme case of p = 1, equivalently the SEP with an infinitely strong localized source, the radius R<sub>d</sub>(t) of this ball grows according to Ref. [42]

$$
R_d(t) \sim \begin{cases} \sqrt{t \ln t} & d = 1, 2, 3 \\ t^{2/d} & d \ge 4 \end{cases} . \qquad (180)
$$

B. Quiet birth

The fertile site plays a special role. This observation sug-
gests amending the multiplication at the fertile site, while not
altering the hopping rules. Thus we return to noninteracting
RWs, but assume that the birth may occur only in a non-
crowded environment. The simplest implementation allows
birth only when a single particle occupies the fertile site
(Fig. 5). To make closer contact with the model of Sec. VIII A,
we adopt a slightly different rule. We assume that if the
fertile site is occupied by a single particle and this hops to
a neighboring site, it leaves behind the daughter particle with

probability p. A successful multiplication event is illustrated
in Fig. 5. The birth rule implies an indirect interaction between
RWs: The large concentration of particles in the proximity of
the fertile site suppresses the birth events.

Let us start again with the extreme case of p = 1. The
fertile site is therefore always occupied, but in contrast to the
extreme model studied in Sec. VIII A, the fertile site may host
many particles. This extreme model with noninteracting RWs
was studied in [42]. Intriguingly, this extreme model exhibits
more rich behaviors than the extreme model in the case of
SEP. Different behaviors again emerge depending on whether
d ≦ 2 or d > 2. For instance, the average density at the fertile
site grows indefinitely in low dimensions and saturates when
d > 2:

$$
n_0 \approx \begin{cases} \frac{1}{2} \ln t & d = 1 \\ \ln(\ln t) & d = 2. \\ 2W_d[2W_d - 1]^{-1} \ln(2W_d) & d \ge 3 \end{cases} \quad (181)
$$

Combining (181) and (166) one finds

$$
N(t) \approx \begin{cases} 2\sqrt{\frac{t}{\pi}} \ln t & d = 1 \\ 4\pi \frac{t \ln(\ln t)}{\ln t} & d = 2. \\ 2d[2W_d - 1]^{-1} \ln(2W_d) t & d > 2 \end{cases} \quad (182)
$$

The derivations [42] were not rigorous, but the numeri-
cal support was convincing; even at the critical dimension
$d_c = 2$ where the behaviors are very subtle, involving a re-
peated logarithm, the agreement with numerics [42] was quite
good.

The same analysis as in [42] shows the universality of the
leading behaviors when *d* ≤ 2, that is, the predictions (181)
and (182) for *d* ≤ 2 remain the same for all *p* ∈ (0, 1]. When
*d* > 2, we anticipate qualitatively similar behaviors as in the
model of Sec. VIII A, i.e., the emergence of three regimes
and the validity of Eqs. (177)–(179). Finally, we mention
that in the extreme model the droplet of visited sites grows
asymptotically according to the same law (180) as in the case
of the SEP.

IX. CONCLUSIONS

In Secs. II–VII we have studied noninteracting random
walkers on homogeneous hypercubic lattices with one special
fertile site where RWs can reproduce. In particular, we have
explored the statistics of the total number of RWs. When
d > 2 and µ ≦ µ<sub>d</sub>, the distribution of the total number of
RWs is stationary and given by (139); in the critical case,
the distribution is particularly neat, Eq. (1). When the RW is
recurrent (d ≦ 2), the distribution P<sub>N</sub>(t) approaches a scaling
form (122) that is still unknown.

The region occupied by random walkers is asymptotically
a growing segment in one dimension and a growing disk in
two dimensions. In both cases, we have computed the growth
velocity. It would be interesting to probe the roughness of the
boundary in two dimensions.

Our process is a simple example of a random walk in
a nonhomogeneous environment with a fertile site where
random walkers can reproduce. More pronounced inhomo-
geneities in an environment characterized by spatially varying
quenched growth rates arise in diverse settings ranging from population dynamics to the kinetics of chemical and nuclear reactions [43–45]. These systems tend to exhibit highly non-self-averaging behaviors [46–53]. It would be interesting to apply large deviation techniques to such models and search for universal features in high dimensions, similar to one displayed by the elementary model studied in the present work.

We have also analyzed (Sec. VIII) the influence of interactions in two models. In the first model, the particles undergo the symmetric exclusion process. In the second model, the particles do not directly interact, but the birth allowed only when the fertile site is occupied by a single particle; this introduces subtle collective interactions. The behaviors in these models drastically differ from the behavior of noninteracting RWs, e.g., the growth of the number of particles cannot be faster than linear. Some qualitative features such as the emergence of the critical birth rate when $d > 2$ are similar to RWs, although our arguments in $d > 2$ dimensions are heuristic. The most intriguing feature of these two particular models is the remarkable universality of the low-dimensional behavior ($d \le d_c = 2$): The reproduction rate does not affect the leading behavior. In this respect, the behaviors are simpler than the behaviors of noninteracting random walkers in $d \le d_c = 2$ dimensions.

Our work has started as an attempt to devise a classical analog of an open quantum system [54,55] driven by a localized source of identical bosons. This quantum system exhibits tricky behaviors. An exponential growth occurs when the strength $\Gamma$ of the source exceeds a critical value $\Gamma_d$, while when $\Gamma \le \Gamma_d$ the growth is quadratic in time when $d > 2$; some subtleties occur when $d = 1$ and $d = 2$. Intriguing behaviors of this open quantum system are not fully captured by the classical analog, so it would be interesting to find a better classical analog.

The symmetric exclusion process with multiplication resembles an open quantum system driven by a localized source of spinless lattice fermions [54,56,57]. The behavior of this open quantum system is somewhat simpler than the behavior of the classical system—in the quantum case, there is only one regime, the average number of fermions always exhibits a linear growth, $N = C_d(\Gamma)t$. The behavior of the amplitude is subtle, e.g., $C_d(\Gamma) \to 0$ as $\Gamma \to \infty$ which a signature of the quantum Zeno effect. The interacting classical systems we studied almost exhibit the Zeno effect in dimensions $d \le 2$ where the leading behaviors are independent of the birth rate. However, the subleading behaviors seem normal, namely increasing with the increase of the birth rate.

## ACKNOWLEDGMENTS

We are grateful to Baruch Meerson for fruitful discussions and collaboration on the earlier stages of this project. We also benefited from discussions with Sid Redner. P.L.K. thanks the Institut de Physique Théorique for hospitality and excellent working conditions.

## APPENDIX A: DENSITY IN ONE DIMENSION

Equation (15) encapsulates the Laplace transforms of the densities:

$$ \hat{n}_j(s) = \frac{\Lambda^{|j|}}{\sqrt{s^2 + 4s - \mu}} \quad (A1) $$

with

$$ \Lambda(s) = \frac{s + 2 - \sqrt{s^2 + 4s}}{2} \quad (A2) $$

The inverse Laplace transform reads

$$ n_j(t) = \int_{s_* - i\infty}^{s_* + i\infty} \frac{ds}{2\pi i} \frac{e^{ts} \Lambda^{|j|}}{\sqrt{s^2 + 4s - \mu}} \quad (A3) $$

An integration contour can go along any vertical line in the complex plane such that $s_* = \text{Re}(s)$ is greater than the real part of singularities of the integrand. We are interested in the asymptotic behavior, so we can employ the saddle point technique. First, we rewrite (A3) as

$$ n_j(t) = \int_{s_* - i\infty}^{s_* + i\infty} \frac{ds}{2\pi i} \frac{e^{tf(s)}}{\sqrt{s^2 + 4s - \mu}} \quad (A4) $$

with $f(s) = s + J \ln \Lambda(s)$, where $J = |j|/t$. The saddle point is found from $f'(s_*) = 0$ to give

$$ s_* = -2 + \sqrt{4 + J^2} $$

and take the vertical contour in (A4) passing through the saddle point. Computing the integral we obtain

$$ n_j(t) = \frac{J}{J - \mu} (2\pi t^3 \sqrt{4 + J^2})^{-1/2} e^{-t \mathcal{D}} \quad (A5) $$

with

$$ \mathcal{D} = 2 - \sqrt{4 + J^2} + J \ln \frac{\sqrt{4 + J^2} + J}{2} \quad (A6) $$

The asymptotic (A5) becomes erroneous when $J \le \mu$. The reason is easy to understand: The above computation tacitly assumed that $s_*$ is greater than the real part of the singularities of the integrand in (A4). These singularities are found from $\sqrt{s^2 + 4s} = \mu$, so the rightmost singularity is located at $C_1 = -2 + \sqrt{4 + \mu^2}$. Since $s_* > C_1$ when $J > \mu$, the asymptotic (A5) is applicable in this region.

When $s_* < C_1$, we still take a contour mostly going through the saddle point, but deform it near the real axis. Namely, we take the contour $(s_* - i\infty, s_* - i0)$, then a contour $(s_*, C_1)$ just below the real axis, then a small circle around $C_1$, then the contour $C_1, s_*$) just above the real axis, and finally $(s_* + i0, s_* + i\infty)$. The leading contribution is provided by the circle integral which is computed (there is a simple pole at $s = C_1$) to yield (145).

To justify the computations in Sec. VII A we notice that near the front $J = t^{-1}r \equiv v < \frac{1}{2}\mu$; see (148) and (149) and Fig. 3. Therefore the inequality $J < \mu$ is obeyed and we can indeed use (145).

One can verify that $\mathcal{D}$ given by (A6) is positive when $J > 0$. Therefore the asymptotic (A5) accounts for exponentially small density, i.e., the range where average quantities like the density are not useful.

## APPENDIX B: RECURRENCE (103)

When $n \gg 1$, the recurrence (103) simplifies to

$$ n\beta_n \approx \frac{\mu}{\sqrt{\mu^2 + 4} - 2} \sum_{a=1}^{n-1} \beta_a \beta_{n-a}, \quad (B1) $$
where $\beta_n = v_n/n!$. Using the generating functions

$$B(z) = \sum_{n \ge 1} \beta_n z^n, \quad z \frac{dB}{dz} = \sum_{n \ge 1} n \beta_n z^n, \quad (B2)$$

we rewrite (B1) as

$$z \frac{dB}{dz} \approx \frac{\mu}{\sqrt{\mu^2 + 4 - 2}} [B(z)]^2. \quad (B3)$$

Making a natural guess

$$\beta_n \approx C n^{-\alpha} \beta^n, \quad (B4)$$

we deduce the leading singular behavior of the generating functions

$$B(z) \approx C \frac{\Gamma(1-\alpha)}{(1-\beta z)^{1-\alpha}}, \quad z \frac{dB}{dz} \approx C \frac{\Gamma(2-\alpha)}{(1-\beta z)^{2-\alpha}} \quad (B5)$$

as $1 - \beta z \to +0$. By inserting (B5) into (B3) we get $\alpha = 0$ and also determine the amplitude $C$ to yield

$$v_n \approx \frac{\sqrt{\mu^2 + 4 - 2}}{\mu} n! \beta^n \quad \text{when } n \gg 1. \quad (B6)$$

We emphasize that $\beta$ is an unknown function of $\mu$.

The only solvable case appears to be the $\mu \to \infty$ limit. In this situation the recurrence (103) becomes

$$v_n = \frac{1}{n-1} \sum_{a=1}^{n-1} \binom{n}{a} v_a v_{n-a}.$$

Recalling $\beta_n = v_n/n!$, one gets $\beta_1 = 1$ and

$$(n-1)\beta_n = \sum_{a=1}^{n-1} \beta_a \beta_{n-a}$$

for $n \ge 2$, from which $\beta_n = 1$ leading to $v_n = n!$. Thus

$$\int_0^\infty dz \, z^n \mathcal{P}(z) = n!$$

from which

$$\mathcal{P}(z) = e^{-z}. \quad (B7)$$

The $\mu \to \infty$ limit corresponds to the zero-dimensional situation where the exact solution is known, Eq. (82), whose scaling form is indeed given by (B7).

## APPENDIX C: MISCELLANIES

In this Appendix we discuss more general variants of the models investigated in the main text and outline some other ways to study them.

### 1. Other discretizations

The analysis in the main text dealt with walkers on the hypercubic lattice $\mathbb{Z}^d$. The study of other lattices would be similar. For $d = 1$, the problem has a well-defined limit when the mesh goes to 0, but not so when $d \ge 2$: the naive continuum space equations are singular and one has no choice but to discretize. This raises the question of universality. It is expected that the existence or not of a threshold for $\mu$

and the exponential growth of the population for instance are universal, while the precise numerical factors are not.

As an illustration, we use the rotation invariance that is present in the continuum with a single fertile site to discretize only the radial part of the problem. One convenient choice is a nearest-neighbor random walk on the semi-infinite line

$$\left\{ x_j := \frac{d-1}{2} + j, \ j = 0, 1, 2, \dots \right\} \quad (C1)$$

with jump rates

$$D_{j,j+1} = \frac{2(d-1+j)}{d-1+2j}, \quad D_{j,j-1} = \frac{2j}{d-1+2j} \quad (C2)$$

in dimension $d$.

The origin of the semi-infinite line on which the walker moves and the jumps rates are determined by the normalization condition $D_{j,j+1} + D_{j,j-1} = 2$, and by the condition that the position of the walker $X(t)$ satisfies

$$\langle X(t)^2 \rangle = X(0)^2 + 2dt \quad (C3)$$

so that $X(t)$ behaves like the distance to the origin for a walker on the lattice $\mathbb{Z}^d$ with diffusion constant $D = 1$.

The resulting time evolution of the average density at $j \ge 1$ is governed by

$$\frac{dn_j}{dt} = \frac{2(d-2+j)}{d-3+2j} n_{j-1} + \frac{2(j+1)}{d+1+2j} n_{j+1} - 2n_j. \quad (C4a)$$

The density at the fertile site obeys

$$\frac{dn_0}{dt} = 2 \left( -n_0 + \frac{1}{d+1} n_1 \right) + \mu n_0. \quad (C4b)$$

This system reduces to (9a) and (9b) when $d = 1$, keeping in mind that on the semi-infinite lattice $n_j$, $j \ge 1$, is the sum of the populations at site $j$ and $-j$.

The continuum space limit of the right-hand side of (C4a) under the substitution $n(x) := n_{x/a}$, where $a$ is the physical mesh of the lattice, is

$$a^2 \frac{\partial}{\partial x} \left( \frac{\partial}{\partial x} - \frac{d-1}{x} \right) n(x), \quad (C5)$$

and the dual of the differential operator $\frac{\partial}{\partial x}(\frac{\partial}{\partial x} - \frac{d-1}{x})$ is indeed $\frac{\partial^2}{\partial x^2} + \frac{d-1}{x} \frac{\partial}{\partial x}$, i.e., the radial part of the Laplace operator in dimension $d$ as should be.

It is easy to solve (C4a) and (C4b) for $d = 3$. Making the Laplace transform with respect to time and taking a generating function

$$N(s, z) := \sum_{j=0}^{\infty} z^j \hat{n}_j(s) \quad (C6)$$

leads to

$$N(s, z) = \frac{d}{dz} \left[ z \frac{1 + (\mu - z^{-1})N(s, 0)}{s - (z - 2 + z^{-1})} \right]. \quad (C7)$$

Imposing that $N(s, z)$ be analytic in the unit disk yields

$$N(s, 0) = \frac{2}{s + 2 + \sqrt{s^2 + 4s - 2\mu}} \quad (C8)$$
and

$$N(s, z) = \frac{N(s, 0)}{\left(1 - z^{\frac{s+2-\sqrt{s^2+4s}}{2}}\right)^2}. \quad (C9)$$

The average occupation numbers exhibit exponential growth if and only if $N(s, 0)$ has a pole at some $s > 0$ which occurs if and only if $\mu > 1$, and then the inverse timescale is $(\mu - 1)^2/\mu$. As expected, there is a threshold for exponential growth just like with the $d = 3$ model on the cubic lattice, but the threshold itself, as well as the inverse timescale and the amplitudes are different.

## 2. An effective model for the critical case

The aim of the forthcoming discussion is to give some intuition for the results in (5) and (6). In particular, we explain the origin of the upper critical dimension $d^c = 4$. This number appears in many guises in random walk theory, but here it is neatly explained, as well as the qualitative features of (5) and (6), in terms of last passage times for a single random walker.

For $d \le 2$ the last passage time at the origin is infinite, meaning that the walker visits the origin infinitely many times. For $d > 2$, the walker visits the origin only finitely many times (and then leaves forever) with probability 1. More precisely, it is a classical result that for $d > 2$ the probability that the last passage at the origin occurs after time $t$ decreases like $t^{1-d/2}$. The fact that the exponent decreases with $d$ is intuitive, and the explanation for the precise exponent is straightforward, at least if universality is assumed. For instance, if $A$ is a bounded domain in $\mathbb{R}^d$ containing the origin, the probability that a Brownian motion is in $A$ at time $t$ is immediately seen to scale like $(4\pi t)^{-d/2}$ Volume($A$) for large $t$ (i.e., for $t$ large compared to the square of the linear size of $A$). If $A$ is a small ball for instance, the Brownian will exit soon afterwards, and then never come back with a positive probability if $d > 2$, so that this scaling is also the scaling for a last passage time. In discrete time and space, this scaling is most easily recovered on the body-centered hypercubic lattice, so that each step of the random walk amounts to a step of independent simple random walks along each of the $d$ axes. The probability to be at the origin at time $t = 2l$ (an even integer) factors as $[2^{-2l} (2l)^l l^{d/2}]^d \sim (\pi l)^{-d/2}$. As expected, for $d \le 2$ the decrease at large $l$ is too slow to sum over $l$ but for $d > 2$, setting $Z_d := \sum_{l \ge 0} [2^{-2l} (2l)^l l^{d/2}]^d$ we infer that $Z_d^{-1}$ is the probability of no return (starting from a neighbor of the origin) and $Z_d^{-1} [2^{-2l} (2l)^l l^{d/2}]^d$ is the probability that the last passage at the origin occurs at time $t = 2l$. Of course the same large $t$ behavior of the law of the last passage time can be computed explicitly in our main model. Using $Z_d^{-1}$ to denote the no return probability in any random walk model, the number of passages at the origin follows a geometric distribution of parameter $Z_d^{-1}$, and in particular the average number of passages at the origin is $Z_d$. The critical fertility rate is fixed by the condition $\mu_d Z_d = 1$ [in particular, for our main model, $Z_d$ is one half of the Watson integral (4)], ensuring that the average offspring size is precisely 1.

The important feature of the law of the last passage time is that it has no first moment for $d \in ]2, 4]$ (i.e., $d = 3$ or 4 in the context of the random walk on the hypercubic lattice, but the

same arguments apply for arbitrary $d \in ]2, 4]$ for the models in Appendix C 1, for instance), but it has a first moment for $d > 4$, yielding a timescale. Thus, the following intuitive picture emerges. For $d > 4$ each walker has (in average) a finite time to reproduce and then leaves forever (so that as long as he is around, he contributes a finite density at the origin). He leaves behind him (in average) one new walker which does the same. So overall the number of walkers grows linearly in time (in average) and the average density ad the origin remains finite forever. However, for $d \in ]2, 4]$, a walker may take an arbitrary long time before it leaves some offspring, explaining at the same time why the density at the origin goes to 0 at large times and why the average growth of the number of walkers is sublinear.

We can elaborate on this qualitative argument by introducing an effective model: instead of a fertility rate ensuring that the average size of the offspring of a walker is 1, the effective model assumes that at his last passage at the origin (and only then) each walker leaves a single descendant. Though the above qualitative discussion applies more directly to this effective model than to the original model, the detailed mathematical analysis is more involved than the simple computations in the fertile random walk model, because the Markov property is replaced by a weaker renewal property. For instance, the study of $n_0(t)$ involves the renewal theorem and the ergodic theorem. An explicit computation does not fit here, but the results of (5) are recovered. The case of $N(t)$ is simpler: the time $t_N$ at which the population reaches $N+1$ is the sum of $N$ independent random variables, each with the law of the last passage time. Thus we can apply an appropriate law of large numbers, the precise result depending only on the large $t$ behavior of the law of the last passage time recalled above. For $d > 4$, there is an average and $t_N/N$ converges with probability 1 to the average of the last passage time, the usual strong law of large numbers. Then by inversion the number of walkers grows like $N(t) \propto t$. On the other hand, for $d \in ]2, 4]$, $t_N$ receives its dominant contribution from some very large summands, changing the scaling. For $d \in ]2, 4[$, $N(t) \propto t^{d/2-1}$ while at $d = 4$ there is a logarithmic correction and $N(t) \propto t/\log t$. The similarity with (6) is thus complete. Heuristic arguments indicate that the relation $n_0(t) \simeq Z_d N(t)/t$ holds at large $t$. The factor $Z_d$ appears as the average number of passages of a random walker at the origin. As $\mu_d Z_d = 1$ the precise relationship between $n_0(t)$ and $N(t)$ exhibited in (5) and (6) is presumably universal.

## 3. Other birth functions

The main text concentrates on the simplest reproduction mechanism, when an individual gives birth to another one, equivalently dies while giving birth to two new individuals. A more general reproduction pattern would be to have a jump rate $\mu(n)$ to die and leave $n$ new individuals for $n = 0, 2, 3, \dots$. It is useful to recast these rates in a generating function $E(z) := \sum_{n \neq 1} \mu(n) z^n$. The rate $\mu(0)$ covers the possibility to die without leaving any offspring. The rate $\mu(2)$ is what was called $\mu$ in the main text; in the general situation we set

$$\mu := E'(1) - E(1) = \sum_{n \neq 1} (n-1)\mu(n). \quad (C10)$$
The generalization of many results to this more general setting is straightforward though cumbersome and less explicit: with the binary reproduction rule, many things can be computed explicitly by solving a quadratic equation, while in the general case one relies on the (implicit) inversion of monotonous functions.

### 4. More general models

We consider a more general Markov model for multiplication and diffusion. Models with several fertile sites were studied for instance in Refs. [7,8,12], for walkers on a lattice, but sometimes in a semi-Markovian context. The lattice structure is crucial for some sharp probabilistic estimates, but for the generalities below, the natural setting is an arbitrary Markov process with countable state space. The sites $j \in A$ (a countable set) each come with their own offspring rate function

$$E_j(z) := \sum_{n \neq 1} \mu_j(n) z^n \quad (C11)$$

with walkers jumping from site $j$ to site $k$ with rates $K_{jk}$. To be consistent with the main text, we set

$$2D_j = -K_{jj} := \sum_{k \in A, k \neq j} K_{jk}. \quad (C12)$$

Thus each walker at site $j$ carries two independent exponential clocks, one for offspring with parameter $E_j(1)$ and one for diffusion with parameter $2D_j$. If the offspring clock rings first (probability $E_j(1)/[E_j(1) + 2D_j]$), the walker dies and leaves $n$ new individuals at site $j$ (each with its new pair of independent clocks) with probability $\mu_j(n)/E_j(1)$, while if the diffusion clock rings first (probability $2D_j/[E_j(1) + 2D_j]$), the walker jumps to site $k \neq j$ (and starts a new pair of independent clocks) with probability $K_{jk}/(2D_j)$.

An observable carrying the one-time information is the generating function

$$\mathfrak{H}(z_{\bullet}, t) := \left\langle \prod_{j \in A} z_j^{N_j(t)} \right\rangle. \quad (C13)$$

Here $z_j$ are independent variables, $N_j(t)$ is the population of the site $j$ at time $t$ and $N(t) := \sum_{j \in A} N_j(t)$ denotes the total population. From the Markov property, one infers the master equation

$$\frac{\partial \mathfrak{H}}{\partial t} = \sum_{j \in A} \left[ E_j(z_j) - E_j(1)z_j + \sum_{k \in A} K_{jk} z_k \right] \frac{\partial \mathfrak{H}}{\partial z_j}. \quad (C14)$$

As usual, such a first-order PDE can be reduced to a family of ODEs by the method of characteristics: if $z_{\bullet}(t)$ solves the system of ordinary differential equations

$$\frac{dz_j(t)}{dt} = E_j[z_j(t)] - E_j(1)z_j(t) + \sum_{k \in A} K_{jk} z_k(t) \quad (C15)$$

with initial conditions $z_{\bullet}(0) = z_{\bullet}$, the generating function is $\mathfrak{H}(z_{\bullet}, t) = \mathfrak{H}(z_{\bullet}(t), 0)$. Solving (C15) is a formidable task in general. An exception is when $A$ is a singleton and the offspring function is simply $E(z) = \mu(0) + \mu(2)z^2$. In the even

simpler case $E(z) = \mu z^2$, one retrieves formula (83) with the substitution $z = e^{\lambda}$.

### 5. Asymptotic number of walkers

If no death is possible, i.e., if $\mu_j(0) = 0$ for every $j \in A$, the population may only increase and it is obvious that $N(t)$ has a (sample by sample) limit at large times $N(\infty)$, which is possibly infinite (this may happen in the supercritical regime). This was used in the main text. Under mild assumptions, $N(\infty)$ remains well defined even if death is possible at some sites: the situation when the random process $N(t)$ oscillates, returning to some minimum $N_*$ at arbitrary large times without ever stabilizing to this value has probability zero. The intuition is that each time the total population returns to the value $N_*$, there is some probability that the next change of population will be a decrease because some walker may diffuse to a site where death is possible. So the fact that the next transition is an increase of population costs some phase space. Intuitively, it is like playing head and tails: even if the probability to toss head is very small, the probability that only tail shows up forever is 0. The difference here is that the different tosses are not independent, and also the bias of the coins may vary from one toss to the next. But if the rates for offspring and diffusion satisfy certain bounds, this annoyance can be controlled.

In particular, this happens when there is a single fertile site, and we concentrate on this situation now. Let $0 \in A$ be the label of the fertile site. Set $D_0 = 1$ for the diffusion constant at 0 to make contact with the notations from the main text. Also set $E(z) = E_0(z)$. If $j \neq 0$, let $R_j$ denote the probability that a walker started at $j$ never returns to the fertile site 0. These probabilities are characteristics of the diffusion on $A$ and do not involve the offspring function. For the site 0, set

$$R = R_0 := \frac{1}{2} \sum_{j \neq 0} R_j K_{0j}. \quad (C16)$$

The computation of the $R_j$ is complicated in general. As an example when the result is simple, the model with jump rates (C2) for $d=3$ leads to

$$R = R_0 = 1/2 \quad R_j = j/(j+1) \quad \text{for } j = 1, 2, \dots \quad (C17)$$

If the process starts with a single walker at 0, the Markov property implies that the generating function $\Pi_0(w) := \sum_n \text{Prob}(N(\infty) = n)w^n$ satisfies

$$[E(1) + 2R]\Pi_0(w) = 2Rw + E[\Pi_0(w)]. \quad (C18)$$

The quantities $E(z)$ and $D$ are input data. The computation of $R$ may be quite involved as already mentioned, but if $R$ is known, (C18) determines $\Pi_0(w)$ either locally via a formal power series expansion or globally via the functional equation itself. For instance, to show the uniqueness of the perturbative expansion, it is enough to do so for the first term, which follows from the fact that $E(z) - [E(1) + 2R]z$ is convex on $[0,1]$, $\ge 0$ at 0 and $\le 0$ at 1 (even $<0$ if $R > 0$). Of course, if $E$ is quadratic (the only offspring is none or twins), $\Pi_0(w)$ is obtained simply by taking the appropriate branch of the solution of a quadratic equation.
If the process starts with single walker at $j \in A$, the analo-
gous function $\Pi_j(w)$ satisfies

$$
\Pi_j(w) = R_j w + (1 - R_j) \Pi_0(w) \text{ for } j \neq 0 \quad (C19)
$$

by the Markov property again. Then

$$
\prod_{j \in A} \Pi_j(w)^{N_j(0)} \qquad (C20)
$$

is the generating function for an asymptotic state with a given
number of individuals for a general initial condition [with
$N(0) < +\infty$]. Even if $\Pi_0(w)$ and the $R_j$s are known explic-
itly, this infinite product is not an elementary function.

If $R = R_0 = 0$, i.e., if a walker leaving 0 returns there with
probability 1, then any site $j$ that has a finite probability to
be visited by the walker has $R_j = 0$ so if $R = 0$ we may as-
sume that $R_j = 0$ for $j \in A$. Then $\Pi_j(w) = \Pi_0(w) = \Pi_0(0)$
is $w$-independent and the study of the asymptotics reduces to
a zero-dimensional analysis: $E(1)\Pi_0(w) = E(\Pi_0(w))$ is the
familiar equation from birth-death processes.

The functional equation (C18) determines the condition for
criticality. Because $N(\infty)$ is well defined,

$$
\mathrm{Prob}(N(\infty) = \infty) = 1 - \prod_{j \in A} \Pi_j(1)^{N_j(0)}. \quad (C21)
$$

The supercritical regime corresponds to $\Pi_0(1) < 1$. If $\Pi_0(1) = 1$ and the derivatives of $\Pi_0(w)$ are finite at $w = 1$, the model is an a subcritical regime. In the generic case, the boundary separating the supercritical and the subcritical regime is $\Pi_0(1) = 1$ and $\Pi'_0(1) = \infty$ (the divergence of a higher derivative while $\Pi'_0(1)$ remains finite would indicate a multicritical point). Taking $w \to 1^-$ in the derivative of (C18),

$$
[E(1) + 2R]\Pi'_0(w) = 2R + E'[\Pi_0(w)]\Pi'_0(w), \quad (C22)
$$

and using the definition of $\mu$ in (C10) leads to the criticality
criterion

$$
\mu_c = E'(1) - E(1) = 2R. \quad (C23)
$$

For criticality conditions when the walkers hop on a lattice;
see Refs. [9–12]. For the models studied in the main text, we
recover the well-known interpretation of the Watson integral
in $d \ge 3$ as the inverse of the return probability to the origin
(starting from the origin, or from any nearest neighbor of the
origin) on the hypercubic lattice. Finally, $\mu_c = 1$ for the model
with jump rates (C2) is also recovered correctly as $R = 1/2$ in
that case.

The fact that $N(\infty)$ is well defined has a number of im-
portant consequences. To mention only one, (C20) can be
rephrased as

$$
\langle w^{N(\infty)} \rangle = \prod_{j \in A} \Pi_j(w)^{N_j(0)}. \quad (C24)
$$

Then the Markov property implies that the process

$$
U(t, w) := \prod_{j \in A} \Pi_j(w)^{N_j(t)} \quad (C25)
$$

is what is called in probability theory a closed martingale (see,
e.g., Refs. [58–60]), i.e., a quantity conserved on average and

converging sample by sample at large times—not only is the
expectation time independent

$$
U(0, w) = \langle U(t, w) \rangle = \langle w^{N(\infty)} \rangle, \quad (C26)
$$

but even

$$
\lim_{t \to \infty} U(t, w) = w^{N(\infty)}. \qquad (C27)
$$

In fact,

$$
\langle U(t, w) \rangle = \Im(z_{\bullet} = \Pi_{\bullet}(w), t), \quad (C28)
$$

and it is instructive (if tedious) to check that the time indepen-
dence of $\langle U(t, w) \rangle$ is also a consequence of (C14).

When $R > 0$, $U(t, w)$ depends on $w$ and is a generating
function for conserved quantities. But as a basic application
of such conserved quantities, we content to compute the law
of the maximal population when $R_j = 0$ for $j \in A$ so that
there is no $w$-dependence. Then $U(t, w) = x^{N(t)}$ where $x$ is
the extinction probability of a walker starting at 0 (or at
any $j \in A$ because $R_j = 0$ for $j \in A$ by assumption). In the
identity

$$
\lim_{t \to \infty} x^{N(t)} = w^{N(\infty)} \quad (C29)
$$

for $w \in [0, 1]$ the left-hand side is $w$-independent, and so
must be the right-hand side. Thus $N(\infty)$ is either 0 or $\infty$ and
$w^{N(\infty)} = \mathbf{1}_{N(\infty)=0}$. Thus (C26) implies

$$
x^{N(0)} = \langle x^{N(t)} \rangle = \langle \mathbf{1}_{N(\infty)=0} \rangle. \quad (C30)
$$

The martingale property is robust: under mild assumption, (C30) holds not only for deterministic times, but also for random times. Thus fix a (large) time horizon $T$ and let $\tau_n$ be the minimum of $T$ and the smallest time at which the number of RWs reaches $n$ (which we take to be infinite is this never occurs). Note that

$$
x^{N(\tau_n)} = x^n \mathbf{1}_{\tau_n < T} + x^{N(T)} \mathbf{1}_{\tau_n = T}. \quad (C31)
$$

Taking the average

$$
x^{N(0)} = x^n \mathrm{Prob}(\tau_n < T) + \langle x^{N(T)} \mathbf{1}_{\tau_n = T} \rangle. \quad (C32)
$$

Now

$$
\lim_{T \to \infty} \mathbf{1}_{\tau_n < T} = \mathbf{1}_{\sup_t N(t) \ge n}, \quad (C33)
$$

while

$$
\lim_{T \to \infty} x^{N(T)} \mathbf{1}_{\tau_n = T} = \mathbf{1}_{N(\infty)=0} \mathbf{1}_{\sup_t N(t) < n}, \quad (C34)
$$

and the right-hand side is simply $\mathbf{1}_{\sup_t N(t) < n}$ because on the
event $\sup_t N(t) < n$, automatically $N(\infty) = 0$. Taking the
$T \to \infty$ limit of (C32) and rearranging gives

$$
\mathrm{Prob}(\sup_t N(t) \ge n) = \frac{1 - x^{N(0)}}{1 - x^n} \quad \text{for } n \ge N(0), \quad (C35)
$$

which has a scale-invariant limit in the critical limit when the
extinction probability goes to 0, namely,

$$
\mathrm{Prob}_{\mathrm{crit}}(\sup_t N(t) \ge n) = \frac{N(0)}{n} \quad \text{for } n \ge N(0). \quad (\mathrm{C36})
$$
[1] G. N. Watson, Q. J. Math. Oxford os-10, 266 (1939).

[2] S. Redner and K. Kang, Phys. Rev. A 30, 3362 (1984).

[3] D. ben-Avraham, S. Redner, and Z. Cheng, J. Stat. Phys. 56, 437 (1989).

[4] S. Albeverio, L. V. Bogachev, and E. B. Yarovaya, C. R. Acad. Sci. Paris 326, 975 (1998).

[5] S. Albeverio and L. V. Bogachev, Positivity 4, 41 (2000).

[6] E. B. Yarovaya, Theory Probab. Appl. 55, 661 (2011).

[7] Ph. Carmona and Y. Hu, Ann. Inst. H. Poincaré, Probab. Stat. 50, 327 (2014).

[8] E. VI. Bulinskaya, Stoch. Proc. Appl. 128, 2325 (2018).

[9] V. A. Topchii, V. A. Vatutin, and E. B. Yarovaya, Theory Probab. Math. Statist. 69, 1 (2004).

[10] V. A. Vatutin and V. A. Topchii, Theory Probab. Appl. 49, 498 (2004).

[11] Y. Hu, V. A. Vatutin, and V. A. Topchii, Theory Probab. Appl. 56, 193 (2012).

[12] E. VI. Bulinskaya, Theory Probab. Appl. 55, 120 (2011).

[13] A. Erdélyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi, Tables of Integral Transforms (McGraw-Hill, New York, 1953), Vol. 1.

[14] K. Ray, arXiv:1409.7806.

[15] T. E. Harris, The Theory of Branching Processes (Dover, New York, 1989).

[16] K. B. Athreya and P. E. Ney, Branching Processes (Dover Publications, Mineola, NY, 2004).

[17] P. Haccou, P. Jagers, and V. A. Vatutin, Branching Processes: Variation, Growth, and Extinction of Populations (Cambridge University Press, New York, 2005).

[18] P. L. Krapivsky, S. Redner, and E. Ben-Naim, A Kinetic View of Statistical Physics (Cambridge University Press, Cambridge, 2010).

[19] P. L. Krapivsky, J. F. F. Mendes, and S. Redner, Eur. Phys. J. B 4, 401 (1998); Phys. Rev. B 59, 15950 (1999).

[20] P. L. Krapivsky, Phys. Rev. E 85, 011152 (2012).

[21] H. Lacoin, F. Simenhaus, and F. L. Toninelli, J. Eur. Math. Soc. 16, 2557 (2014).

[22] M. Eden, in Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability. Vol. IV: Biology and Problems of Health, edited by Jerzy Neyman (University of California Press, Berkeley, Los Angeles, 1961), p. 223.

[23] D. Richardson, Proc. Cambridge Philos. Soc. 74, 515 (1973).

[24] H. Kesten, Aspects of first passage percolation, in École d'Été de Probabilités de Saint Flour XIV - 1984, edited by P. L. Hennequin, Lecture Notes in Mathematics (Springer-Verlag, Berlin, 1986), Vol. 1180, pp. 125-264.

[25] M. Bramson, Convergence of Solutions of the Kolmogorov Equation to Traveling Waves (American Mathematical Society, Providence, RI, 1983).

[26] E. Brunet and B. Derrida, Phys. Rev. E 56, 2597 (1997).

[27] U. Ebert and W. van Saarloos, Phys. Rev. Lett. 80, 1650 (1998); Physica D 146, 1 (2000).

[28] W. van Saarloos, Phys. Rep. 386, 29 (2003).

[29] P. L. Krapivsky and S. N. Majumdar, Phys. Rev. Lett. 85, 5492 (2000); S. N. Majumdar and P. L. Krapivsky, Phys. Rev. E 65, 036127 (2002).

[30] S. N. Majumdar and P. L. Krapivsky, Physica A 318, 161 (2003).

[31] H. Spohn, Large Scale Dynamics of Interacting Particles (Springer-Verlag, New York, 1991).

[32] C. Kipnis and C. Landim, Scaling Limits of Interacting Particle Systems (Springer, New York, 1999).

[33] G. M. Schütz, Exactly Solvable Models for Many-Body Systems Far from Equilibrium, Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic Press, London, 2000), Vol. 19.

[34] R. A. Blythe and M. R. Evans, J. Phys. A 40, R333 (2007).

[35] B. Derrida, J. Stat. Mech.: Theory Exp. (2007) P07023.

[36] T. Chou, K. Mallick, and R. K. P. Zia, Rep. Prog. Phys. 74, 116601 (2011).

[37] P. L. Krapivsky, Phys. Rev. A 45, 1067 (1992).

[38] L. Frachebourg and P. L. Krapivsky, Phys. Rev. E 53, R3009 (1996).

[39] M. Mobilia, Phys. Rev. Lett. 91, 028701 (2003).

[40] J. E. Santos and G. M. Schütz, Phys. Rev. E 64, 036107 (2001).

[41] P. L. Krapivsky, Phys. Rev. E 86, 041103 (2012).

[42] P. L. Krapivsky and D. Stefanovic, J. Stat. Mech.: Theory Exp. (2014) P09003.

[43] Y. C. Zhang, Phys. Rev. Lett. 56, 2113 (1986).

[44] Y. B. Zeldovich, S. A. Molchanov, A. A. Ruzmaikin, and D. D. Sokolov, Usp. Fiz. Nauk. 152, 3 (1987) [Sov. Phys. Usp. 30, 353 (1987)].

[45] S. A. Molchanov, Acta Appl. Math. 22, 139 (1991).

[46] W. Ebeling, A. Engel, B. Esser, and R. Feistel, J. Stat. Phys. 37, 369 (1984); A. Engel and W. Ebeling, Phys. Rev. Lett. 59, 1979 (1987).

[47] M. N. Rosenbluth, Phys. Rev. Lett. 63, 467 (1989).

[48] H. Leschke and S. Wonneberger, J. Phys. A 23, 1475 (1990).

[49] R. Tao, Phys. Rev. A 43, 5284 (1991).

[50] K. A. Dahmen, D. R. Nelson, and N. M. Shnerb, J. Math. Biol. 41, 1 (2000).

[51] M. M. Desai and D. R. Nelson, Theor. Popul. Biol. 67, 33 (2005).

[52] P. L. Krapivsky and K. Mallick, J. Stat. Mech.: Theory Exp. (2011) P01015.

[53] T. Gueudré and D. G. Martin, Europhys. Lett. 121, 68005 (2018).

[54] M. Butz and H. Spohn, Ann. Henri Poincaré 10, 1223 (2010).

[55] P. L. Krapivsky, K. Mallick, and D. Sels, J. Stat. Mech.: Theory Exp. (2020) 063101.

[56] H. Fröml, A. Chiocchetta, C. Kollath, and S. Diehl, Phys. Rev. Lett. 122, 040402 (2019).

[57] P. L. Krapivsky, K. Mallick, and D. Sels, J. Stat. Mech.: Theory Exp. (2019) 113108.

[58] J. L. Doob, Stochastic Processes (Wiley, New York, 1953).

[59] J. Jacod and P. Protter, Probability Essentials (Springer-Verlag, Berlin, 2004).

[60] D. Williams, Probability with Martingales (Cambridge University Press, Cambridge, 1991).