Random walk through a fertile site

TL;DR

This paper analyzes the growth dynamics of random walkers multiplying at a fertile site on hyper-cubic lattices, revealing dimension-dependent growth rates, critical thresholds, and universal distribution properties, with implications for interacting particle systems.

Contribution

It provides a detailed characterization of the growth behavior, moments, and distribution of walkers, introducing critical thresholds and universal properties across dimensions.

Findings

Exponential growth in 1D and 2D for certain rates

Finite walker number in high dimensions below critical rate

Universal distribution in the critical regime

Abstract

We study the dynamics of random walks hopping on homogeneous hyper-cubic lattices and multiplying at a fertile site. In one and two dimensions, the total number $\mathcal{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\leq \mu_d$. We determine $\mu_d$ and show that $\langle\mathcal{N}(t)\rangle$ grows exponentially if $\mu>\mu_d$. The distribution of the total number of walkers remains broad when $d\leq 2$, and also when $d>2$ and $\mu>\mu_d$. We compute $\langle \mathcal{N}^m\rangle$ explicitly for small $m$, and show how to determine higher moments. In the critical regime, $\langle \mathcal{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 \mathcal{N}^m\rangle\sim \langle \mathcal{N}\rangle^{2m-1}$, in the critical regime; the growth is normal, $\langle \mathcal{N}^m\rangle\sim \langle \mathcal{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\mathcal{N}(t)\rangle$ grows linearly (not exponentially) in time; when $d\leq d_c=2$, the leading behavior is independent on $\mu$ and $\langle\mathcal{N}(t)\rangle$ exhibits a sub-linear growth.

Equations

$$P_N(\infty )=\underset{t\to \infty }{\lim }{P}_N(t)=\frac{1}{\sqrt{4\pi }}\frac{\Gamma \left(N-\frac{1}{2}\right)}{\Gamma (N+1)}$$

$$N(t)\equiv ⟨\mathcal{N}(t)⟩\sim e^{C_{d}t}$$

$${\mu }_d=\frac{2}{W_d}$$

$$W_d={\int }_0^{\infty }dx{\left[e^{-x}{I}_0(x)\right]}^d$$

$$n_0(t)\simeq {\pi }^{-2}\times \{\begin{array}{ll}{\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{array}$$

$$N(t)\simeq {\pi }^{-2}\times \{\begin{array}{ll}(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{array}$$

$$\text{Prob}[\mathcal{N}(\infty )=\text{ finite}]=\{\begin{array}{ll}0& d\le 2\\ \frac{{\mu }_d}{\mu }& d>2,\mu >{\mu }_d\\ 1& d>2,\mu \le {\mu }_d\end{array}$$

$$⟨\mathcal{N}(\infty )⟩=\frac{{\mu }_d}{{\mu }_d-\mu }$$

$$\frac{d{n}_j}{dt}=n_{j-1}-2n_j+n_{j+1}$$

$$\frac{d{n}_0}{dt}=n_{-1}-2n_0+n_1+\mu n_0$$

$${\hat{n}}_k(s)={\int }_0^{\infty }dt{\mathrm{e}}^{-st}{n}_k(t)$$

$$N(s,q)=\sum _{k=-\infty }^{\infty }{\hat{n}}_k(s)e^{-\mathrm{i}qk}$$

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

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

$${\int }_0^{2\pi }\frac{dq}{2\pi }\frac{1}{s+2(1-\cos q)}=\frac{1}{\sqrt{s^2+4s}}$$

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

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

$$\hat{N}(s)=\sum _{k=-\infty }^{\infty }{\hat{n}}_k(s)=N(s,q=0)$$

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

$$n_0(t)={\int }_{s_{\ast }-\mathrm{i}\infty }^{s_{\ast }+\mathrm{i}\infty }\frac{ds}{2\pi \mathrm{i}}\frac{e^{ts}}{\sqrt{s^2+4s}-\mu }$$

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

$$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}}$$

$$e^{2t}{n}_0(t)\simeq 2{\int }_0^{t}d\tau I_1(2\tau )e^{\mu \sqrt{t^2-{\tau }^2}}$$

$$I_1(x)\simeq \frac{e^x}{\sqrt{2\pi x}}\text{when}x\gg 1$$

$$e^{2t}{n}_0(t)\simeq \sqrt{\frac{t}{\pi }}{\int }_0^1\frac{d\eta }{\sqrt{\eta }}{e}^{tf(\eta )}$$

$$n_0(t)\simeq A_1{e}^{t{C}_1}$$

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

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

$$N(t)=1+\mu {\int }_0^{t}d\tau n_0(\tau )$$

$$N(t)\simeq \frac{{\mu }^2}{{\mu }^2+4-2\sqrt{{\mu }^2+4}}{e}^{t{C}_1}$$

$$N(t)\simeq \{\begin{array}{ll}2e^{{\mu }^{2}t/4}& \mu \ll 1\\ e^{(\mu -2)t}& \mu \gg 1\end{array}$$

Full text

I Introduction

We study non-interacting random walks (RWs) on the hyper-cubic 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 hyper-cubic 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 counter-intuitive 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 $\mathcal{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:

[TABLE]

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

[TABLE]

The growth rate $C_{d}$ plays the role of the Malthusian parameter; we computed $C_{d}$ for hyper-cubic lattices $\mathbb{Z}^{d}$.

In the following we consider only hyper-cubic 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

[TABLE]

where

[TABLE]

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

[TABLE]

while the average number of RWs grows as

[TABLE]

The amplitudes $A_{d}$ are given in Eq. (64).

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.

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

[TABLE]

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

[TABLE]

In physics literature, our problem has been examined in RK84 ; bARZ89 . Our results are much more detailed, e.g., only average characteristics have been probed in RK84 ; bARZ89 . Several generalizations, e.g., biased RWs, systems with a few fertile sites, etc. have been additionally studied in Refs. RK84 ; bARZ89 . 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 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 Sergio98 ; Sergio00 ; Yarovaya ; Carmona ; Bul18 . 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 Vat1 ; Vat2 ; Vat3 ; Bul11 . 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 RK84 ; bARZ89 ; Carmona ; Bul18 ; Bul11 . 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 RK84 ; bARZ89 ; Sergio98 ; Sergio00 ; Yarovaya ; Carmona ; Bul18 ; Vat1 ; Vat2 ; Vat3 ; Bul11 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), (8) and asymptotically exact results like (6)–(8). In Sec. V we study fluctuations, e.g., we compute the moments $\langle\mathcal{N}^{2}\rangle$ and $\langle\mathcal{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\leq\mu_{d}$, and we determine it. In Sec. VII, we show that when $d\leq 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 non-interacting RWs. A few technical calculations are relegated to Appendices A–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

[TABLE]

Making the Laplace transform with respect to time

[TABLE]

and the Fourier transform with respect to lattice sites

[TABLE]

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

[TABLE]

Using

[TABLE]

and the identity

[TABLE]

we extract from (12) the Laplace transform of the density at the fertile site

[TABLE]

Thus

[TABLE]

The Laplace transform of the average number of RWs

[TABLE]

is therefore given by

[TABLE]

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

[TABLE]

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 $\widehat{f}(s)$. We actually want to determine the inverse Laplace transform of \widehat{f}\big{(}\sqrt{s^{2}-a^{2}}\big{)}, and there is an expression through $f(t)$ which is valid for arbitrary $f(t)$. It reads Bateman

[TABLE]

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 $\widehat{f}=1/(s-\mu)$, the corresponding inverse Laplace transform is $f=e^{\mu t}$. Using these relations together with (19) we obtain

[TABLE]

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

[TABLE]

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

[TABLE]

we simplify (21) to

[TABLE]

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

[TABLE]

with parameters

[TABLE]

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 right-most 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

[TABLE]

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

[TABLE]

Combing (27) with the asymptotic (24) we obtain

[TABLE]

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

[TABLE]

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_{1}n_{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,

[TABLE]

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

[TABLE]

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

[TABLE]

for the square lattice.

Applying the Laplace-Fourier transform

[TABLE]

to (31) we obtain

[TABLE]

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

[TABLE]

To compute the integral we use the identity

[TABLE]

where

[TABLE]

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

[TABLE]

where we use the shorthand notation

[TABLE]

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

[TABLE]

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

[TABLE]

with $C_{2}$ determined from $\Phi(C_{2})=\mu$. Using (38) we get

[TABLE]

The average number of RWs is asymptotically

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

into (31) we obtain

[TABLE]

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

[TABLE]

we deduce

[TABLE]

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

[TABLE]

One can express the integrals through generalized hypergeometric functions, see Ray . 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

[TABLE]

IV Average Growth when $d>2$

IV.1 Green function approach

We have an integral equation

[TABLE]

where we have used the shorthand notation

[TABLE]

The exponential growth $n_{\bf 0}(t)\sim e^{C_{d}t}$ is consistent with (51) when the growth rate $C_{d}$ satisfies

[TABLE]

Recalling the definition (4) of the Watson integral and taking into account that $C_{d}\geq 0$ we see that the right-hand side of (52) cannot exceed $\mu W_{d}/2$. Thus

[TABLE]

The asymptotic (22) shows that the Watson integral (4) diverges when $d\leq 2$. This together with the obvious fact that the right-hand side of (52) is a decreasing function of $C_{d}$ shows that (52) 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 (53) is not obeyed, there is no solution to (52). 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 Watson via Euler’s gamma functions:

[TABLE]

Numerically

[TABLE]

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

[TABLE]

when $d\to\infty$.

IV.2 Laplace-Fourier transform

To probe the behavior in the $\mu\leq\mu_{d}$ range, we apply the Laplace-Fourier transform:

[TABLE]

where ${\bf a}=(a_{1},\ldots,a_{d}),~{}{\bf q}=(q_{1},\ldots,q_{d})$ and

[TABLE]

We find

[TABLE]

To avoid cluttered notation we write

[TABLE]

where

[TABLE]

Fixing $\widehat{n}_{\bf 0}(s)$ as before we arrive at

[TABLE]

Note that

[TABLE]

remains positive when $d>2$.

IV.2.1 Supercritical regime: $\mu>\mu_{d}$

The growth is exponential, $n_{\bf 0}(t)\sim e^{s_{d}t}$, with $s_{d}$ following from

[TABLE]

It is straightforward to verify that $s_{d}=C_{d}$ which was determined by (52).

IV.2.2 Subcritical regime: $\mu<\mu_{d}$

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

[TABLE]

Hence $\lim_{s\to 0}s\widehat{N}(s)=N_{d}$ which is consistent with (59c) when $N_{d}=\frac{\Phi_{d}(0)}{\Phi_{d}(0)-\mu}$. Using (60) we arrive at the general formula (8) for the average number of RWs.

IV.2.3 Critical regime: $\mu=\mu_{d}$

Using the definition (58) one finds that

[TABLE]

when $s\to+0$ with

[TABLE]

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

IV.3 Density

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

[TABLE]

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.

IV.3.1 Supercritical regime: $\mu>\mu_{d}$

The normalized density is asymptotically stationary in this regime

[TABLE]

The generating function

[TABLE]

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

[TABLE]

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

[TABLE]

IV.3.2 Critical regime: $\mu=\mu_{d}$

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

[TABLE]

where $\mu_{d}=\mu_{d}/D$. There are no simple general expressions for $m_{\bf a}$ valid for all ${\bf a}\in\mathbb{Z}^{d}$. The normalized density at sites neighboring the fertile site admits a simple expression through the Watson integral

[TABLE]

Noting that $m_{\bf a}$ satisfies a discrete Poisson equation

[TABLE]

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

[TABLE]

The average number of RWs is therefore

[TABLE]

The cutoff length is expected to grow diffusively with time: $R(t)\sim\sqrt{t}$. Hence $N(t)\sim tn_{\bf 0}(t)$. This is consistent with (5)–(6).

V Fluctuations

We are mostly interested in $\mathcal{N}(t)$, the total number of RWs. Focusing on this global quantity allows us to suppress spatial aspects. The evolution of $\mathcal{N}(t)$ can be interpreted as a branching process teh ; athreya04 ; vatutin . This change of view greatly helps in calculations.

V.1 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}[\mathcal{N}(t)=N]$ be the probability distribution of the total number of RWs. The moment generating function

[TABLE]

satisfies an integral equation

[TABLE]

Here we shortly write

[TABLE]

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)={\rm Prob}(T_{1}\geq t)$ and results in the first term on the right-hand side of (75). 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 (75).

Thus, the problem reduces to solving a non-linear integral equation (75). The probability density ${\rm Prob}(T_{1}=\tau)$ encodes all the geometric data of the problem (the structure of the lattice and the spatial dimension).

V.2 Zero-dimensional case

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

[TABLE]

so the integral equation (75) becomes

[TABLE]

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 0-dimensional case. Indeed, the probabilities $P_{N}(t)$ satisfy exact rate equations

[TABLE]

Solving (79) subject to the initial condition $P_{N}(0)=\delta_{N,1}$ is straightforward (see e.g. book ). The solution reads

[TABLE]

Hence the moment generating function is

[TABLE]

in the 0-dimensional case. One can verify that (81) is indeed the solution of (78) satisfying the initial condition

[TABLE]

V.3 Perturbative expansion

Treating $\lambda$ as a small parameter we write

[TABLE]

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

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

[TABLE]

This linear integral equation can be solved by the Laplace transform if we know the value of the branching probability. To ensure that (75) 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 (84), but with an extra source term

[TABLE]

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

[TABLE]

Performing the Laplace transform of (84) we find

[TABLE]

The Laplace transform of (76b) gives

[TABLE]

and hence (87) simplifies to

[TABLE]

The consistency with (59c) allows us to fix

[TABLE]

V.4 Moments in one dimension

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

[TABLE]

with $\nu_{1}=1+2/\sqrt{\mu^{2}+4}$ and $C_{1}=\sqrt{\mu^{2}+4}-2$, see (25). The exponential growth (91) is consistent with (84) if

[TABLE]

Equivalently, we re-write (92) as $2\widehat{\psi}(C_{1})=1$, and using (90) and $\Phi_{1}(s)=\sqrt{s^{2}+4s}$ we recover (25).

Inserting (91) into (85) we find $M\sim e^{2C_{1}t}$ suggesting us to seek the solution in the form

[TABLE]

Inserting (93) into (85) we find

[TABLE]

The integral in the above equation is $\widehat{\psi}(2C_{1})$, and using (90) we can express the ratio $\nu_{2}/\nu_{1}^{2}$ as

[TABLE]

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

[TABLE]

This ratio decreases from $\sqrt{8}+2=4.828427\ldots$ 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 $\mathcal{N}(t)$ is a non-self-averaging quantity for all $\mu$ and $D$.

The third moment grows according to

[TABLE]

and after straightforward calculations one gets

[TABLE]

which can be re-written similarly to (96):

[TABLE]

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

Generally the $n^{\text{th}}$ moment grows according to

[TABLE]

and the same calculations as above yield

[TABLE]

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

V.5 Moments in higher dimensions

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

V.5.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:

[TABLE]

where

[TABLE]

Equations (95), (98), (101) remain applicable after the obvious replacement $C_{1}\to C_{d}$ and $\Phi_{1}\to\Phi_{d}$. For instance, the recurrence (101) becomes

[TABLE]

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

V.5.2 Critical regime: $\mu=\mu_{d}$

We perform the Laplace transform of (85) and find

[TABLE]

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

[TABLE]

Using additionally $\widehat{\psi}(s)=\mu_{d}/[\mu_{d}+\Phi_{d}(s)],~{}\widehat{\psi}(0)=1/2$ and $\widehat{n}_{\bf 0}(s)=1/[\Phi_{d}(s)-\mu_{d}]$ we simplify (105) to

[TABLE]

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

[TABLE]

We insert (63) and (107) into (106) to yield

[TABLE]

from which

[TABLE]

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,

[TABLE]

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:

[TABLE]

The divergence of the average, $\sum_{N\geq 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

[TABLE]

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

[TABLE]

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 non-stationary and it quickly vanishes. Using these assumptions one can determine the exponents $a$ and $\xi$. Indeed, we estimate two moments

[TABLE]

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\geq 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

[TABLE]

Thus we provided heuristic evidence for the tail

[TABLE]

One extra check of (113) is based on computing higher moments. Using (112)–(113) we find

[TABLE]

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

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 (106) one finds

[TABLE]

A long but straightforward calculation gives

[TABLE]

in agreement with (114). Similarly to (109) we have

[TABLE]

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

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

V.5.3 Subcritical regime: $\mu<\mu_{d}$

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

[TABLE]

VI The probability distribution $P_{N}(t)$

VI.1 One dimension

In Sec. V.4 we have shown that in one dimension the moments $M_{a}(t)=\langle\mathcal{N}^{a}\rangle=\sum_{N\geq 1}N^{a}P_{N}(t)$ satisfy

[TABLE]

where $N(t)=M_{1}(t)$. The growth laws (119) suggest that in the long time limit the probability distribution $P_{N}(t)$ acquires the scaling form

[TABLE]

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

[TABLE]

VI.1.1 Small $N$ behavior

In many problems, the scaling form remains applicable even when $N=O(1)$ and $t\to\infty$, but there are counter-examples, e.g., in sub-monolayer epitaxial growth JM98 . In the present case, $P_{N}(t)$ also exhibits an unusual behavior for small $N$. Inserting (74) into (75) we obtain

[TABLE]

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

[TABLE]

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

[TABLE]

with amplitudes being Catalan numbers.

At first sight, the asymptotically exact results (124a)–(124c) disagree with the scaling form (120). Of course, when $N=O(1)$, the scaling variable $z$ in (120) vanishes when $t\to\infty$, while $z$ must be finite, see (121). 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

[TABLE]

from (124c). The crossover from (125) to (120) apparently occurs when $(Dt)^{-N_{*}/2}\propto 1/N(t)\propto e^{-C_{1}t}$. Thus

[TABLE]

In the boundary layer, $N\ll N_{*}$, the distribution $P_{N}(t)$ varies according to (124c); the scaling apparently emerges when $N\gg N_{*}$. Similar behaviors with a boundary layer structure at small masses were found in models mimicking sub-monolayer epitaxial growth JM98 .

VI.1.2 Large $N$ behavior

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

[TABLE]

The large $n$ behavior of the amplitudes $\nu_{n}$ is established in Appendix B. It gives

[TABLE]

and implies an exponential tail

[TABLE]

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

VI.2 Two dimensions

In two dimensions, the scaling laws (119) hold and the probability distribution $P_{N}(t)$ is also expected to acquire the scaling form (120).

For small $N$ we again rely on Eqs. (122a)–(122d). Using (88) and (90) with $\Phi_{2}$ given by (38) we obtain

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

The scaling form (120) emerges when $N\gg N_{*}$.

VI.3 Dimensions $d>2$

When $\mu\leq\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

[TABLE]

Recalling that $\widehat{\psi}(0)=\mu/(\mu+\mu_{d})$ and solving (134) we obtain

[TABLE]

which is expanded to yield

[TABLE]

Therefore

[TABLE]

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

In the critical regime, $\mu=\mu_{d}$, equation (137) 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 hyper-cubic 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 (120). The moment generating function becomes asymptotically stationary when $\lambda<0$:

[TABLE]

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

[TABLE]

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

[TABLE]

The scaled density satisfies

[TABLE]

VII Spatial characteristics

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

[TABLE]

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.

VII.1 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

[TABLE]

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

[TABLE]

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

[TABLE]

with velocity

[TABLE]

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

[TABLE]

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

We have used (143) for $j\sim t$, i.e., on distances greatly exceeding the diffusion scale, $j\sim\sqrt{Dt}$. The derivation of (143) 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 (143) is known, but different $A_{1}$ arise in different realizations. Fluctuations do not affect the leading behavior, however. Indeed, (144) gives

[TABLE]

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

[TABLE]

with constant fluctuating from realization to realization.

The droplet $\mathcal{D}(t)=[\ell(t),r(t)]$ 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.

VII.2 Two dimensions

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

[TABLE]

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 non-trivial 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 PK-Ising ; Fabio different from the disk. The Eden-Richardson growth model Eden ; Richardson ; Kesten 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 non-trivial 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

[TABLE]

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

[TABLE]

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

[TABLE]

The numerical factor $C=O(1)$ remains undetermined in the realm of continuum framework. A linearly independent solution of Eq. (151) involving another modified Bessel function, I_{0}\big{(}\sqrt{C_{2}}\,r\big{)}, is absent in (152) since this solution diverges when $r\to\infty$. The criterion (144) gives

[TABLE]

where $R$ is the boundary of the droplet. By inserting the large time asymptotic $n_{\bf 0}(t)\sim e^{C_{2}t}$ and

[TABLE]

into (153) we obtain

[TABLE]

in the leading order. The asymptotic behaviors are

[TABLE]

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

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

[TABLE]

A logarithmic correction to the front position is known to occurs (see front1 ; bd ; evs ; van ; km ; mk 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 non-interacting 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.

VIII.1 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 Spohn ; KL99 ; S00 ; BE07 ; D07 ; CKZ ).

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 Spohn ; KL99 ; S00 ; BE07 ; D07 ; CKZ , exactly as in the case of non-interacting RWs. This property is easy to appreciate. In one dimension, for instance, one writes an exact equation

[TABLE]

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. (158) one notices that second order correlation functions like $\langle\tau_{j-1}\tau_{j}\rangle$ cancel. Therefore Eq. (158) simplifies indeed to the lattice diffusion equation

[TABLE]

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, $D=1$, due to our convention that the hopping rates to neighboring sites are equal to unity.

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

[TABLE]

which becomes

[TABLE]

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 (161), we merely have the boundary condition

[TABLE]

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

[TABLE]

Thus in the extreme case we need to solve Eq. (159) subject to (162)–(163). 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 PK92 ; FK96 ; Mauro ; it also obviously describes the SEP with an infinitely strong localized source Santos ; PK-SEP-source ; Darko-source . Several exact and asymptotically exact behaviors are known. For instance, the average $N(t)\equiv\langle\mathcal{N}(t)\rangle$ exhibits an asymptotic growth PK-SEP-source

[TABLE]

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, 54).

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 Santos ; PK-SEP-source : The ratio of the variance to the average is asymptotically

[TABLE]

Let us look at non-extreme 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 (161) 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. (164) in the extreme case, $p=1$, and clearly occurs for all $0<p\leq 1$. Since the flux asymptotically vanishes, the behavior approaches to the behavior of the SEP at equilibrium when the correlators factorize Spohn ; KL99 ; S00 ; BE07 ; D07 . Thus $\langle\tau_{0}\tau_{1}\rangle=n_{0}n_{1}$ is asymptotically exact, so Eq. (161) simplifies to

[TABLE]

Summing (166) and all Eqs. (159) for $j\neq 0$, and taking into account the $n_{j}=n_{-j}$ symmetry, we obtain

[TABLE]

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

[TABLE]

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\geq 0$ and $t\geq 0$ is known PK-SEP-source :

[TABLE]

Using (169) one confirms the asymptotic (168) and gets $A_{j}=j$ in the extreme case.

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

[TABLE]

The amplitude in front of the linear term is fixed by the known asymptotic, $A_{j}\simeq j$ when $j\gg 1$, which can be established by using a continuum approach valid when $j\gg 1$. Substituting (168) into (166) we deduce a relation $A_{0}=(1-p)A_{1}$. Combining this result with (170) 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}\simeq A_{1}/\sqrt{\pi t}$ into (167) and integrating we find

[TABLE]

The leading asymptotic growth is therefore independent on the birth probability $p\in(0,1]$. This remarkable anomaly phenomenon occurs in many branches of science ranging from turbulence (where it is known as a dissipative anomaly, see e.g. Frisch ; Khanin ) to anomalies in quantum field theory (see Bilal and references therein). In the present situation, the anomaly seems particularly tractable and it would be interesting to understand its behavior in detail.

The anomaly seems present only in the leading behavior. To confirm, or disprove, this assertion one would like to compute sub-leading terms. In the extreme case, we know the exact answer PK-SEP-source

[TABLE]

from which one can obtain the entire expansion

[TABLE]

Generally when $p<1$ we anticipate

[TABLE]

with $C_{1}(p)<0$ when $p<1$. A similar expansion has been derived in PK-SEP-source 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 sub-leading term, however, is only logarithmically smaller than the leading term, and it probably depends on $p$. In other words, we anticipate

[TABLE]

with large sub-leading correction, so the convergence to the leading asymptotic is extremely slow. The sub-leading correction in (174) 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 Darko-source .

Thus when $d\leq 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 non-interacting 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\leq p_{c}(d)$. Thus for $p\leq 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:

[TABLE]

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 non-interacting 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

[TABLE]

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

[TABLE]

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 (113). Hence the higher moments $\langle\mathcal{N}^{m}\rangle$ probably exhibit the same dynamical behaviors (114) as in the case of RWs.

The behaviors (175)–(177) are conjectural. The rates $S_{d}(p)$ in (175), the amplitudes in (176)–(177), and the critical birth probabilities $p_{c}(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 $\Pi(N)$, not its asymptotic behavior, so there is no ground for any guess about $\Pi(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_{d}(t)$ of this ball grows according to Darko-source

[TABLE]

VIII.2 Quiet birth

The fertile site plays a special role. This observation suggests amending the multiplication at the fertile site, while not altering the hopping rules. Thus we return to non-interacting 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.1, 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.1, the fertile site may host many particles. This extreme model with non-interacting RWs was studied in Darko-source . 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\leq 2$ or $d>2$. For instance, the average density at the fertile site grows indefinitely in low dimensions and saturates when $d>2$:

[TABLE]

Combining (179) and (164) one finds

[TABLE]

The derivations Darko-source were not rigorous, but the numerical support was convincing; even at the critical dimension $d_{c}=2$ where the behaviors are very subtle, involving a repeated logarithm, the agreement with numerics Darko-source was quite good.

The same analysis as in Darko-source shows the universality of the leading behaviors when $d\leq 2$, that is, the predictions (179)–(180) for $d\leq 2$ remain the same for all $p\in(0,1]$. When $d>2$, we anticipate qualitatively similar behaviors as in the model of Sec. VIII.1, i.e., the emergence of three regimes and the validity of Eqs. (175)–(177). Finally, we mention that in the extreme model the droplet of visited sites grows asymptotically according to the same law (178) as in the case of the SEP.

IX Conclusions

In Secs. II–VII, we have studied non-interacting random walkers on homogeneous hyper-cubic 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 $\mu\leq\mu_{d}$, the distribution of the total number of RWs is stationary and given by (137); in the critical case, the distribution is particularly neat, viz. it is purely algebraic (1). When the RW is recurrent ($d\leq 2$), the distribution $P_{N}(t)$ approaches a still unknown scaling form (120).

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 non-homogeneous environment with a fertile site where random walkers can reproduce. More pronounced inhomogeneities 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 Zhang ; Zeldovichreview ; Molchanov . These systems tend to exhibit highly non-self-averaging behaviors Ebeling ; Rosenbluth ; Leschke ; Tao ; Nelson ; Desai ; PaulKM ; Gueudre . 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 non-interacting 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\leq d_{c}=2$): The reproduction rate does not affect the leading behavior. In this respect, the behaviors are simpler than the behaviors of non-interacting random walkers in $d\leq d_{c}=2$ dimensions.

Our work has started as an attempt to devise a classical analog of an open quantum system Spohn10 ; KMS-bosons 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\leq\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 spin-less lattice fermions Spohn10 ; Kollath ; KMS-fermions . 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\leq 2$ where the leading behaviors are independent of the birth rate. However, the sub-leading 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 benefitted from discussions with Sid Redner. PLK 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:

[TABLE]

with

[TABLE]

The inverse Laplace transform reads

[TABLE]

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 re-write (183) as

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

The asymptotic (185) becomes erroneous when $J\leq\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 (184). These singularities are found from $\sqrt{s^{2}+4s}=\mu$, so the right-most singularity is located at $C_{1}=-2+\sqrt{4+\mu^{2}}$. Since $s_{*}>C_{1}$ when $J>\mu$, the asymptotic (185) 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_{*}-{\rm i}\infty,s_{*}-{\rm i}0)$, 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_{*}+{\rm i}0,s_{*}+{\rm 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 (143).

To justify the computations in Sect. VII.1 we notice that near the front $J=t^{-1}r\equiv v<\frac{1}{2}\mu$, see (146)–(147) and Fig. 3. Therefore the inequality $J<\mu$ is obeyed and we can indeed use (143).

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

Appendix B Recurrence (101)

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

[TABLE]

where $\beta_{n}=\nu_{n}/n!$. Using the generating functions

[TABLE]

we re-write (187) as

[TABLE]

Making a natural guess

[TABLE]

we deduce the leading singular behavior of the generating functions

[TABLE]

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

[TABLE]

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 (101) becomes

[TABLE]

Recalling $\beta_{n}=\nu_{n}/n!$, one gets $\beta_{1}=1$ and

[TABLE]

for $n\geq 2$, from which $\beta_{n}=1$ leading to $\nu_{n}=n!$. Thus

[TABLE]

from which

[TABLE]

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

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.

C.1 Other discretizations

The analysis in the main text dealt with walkers on the hyper-cubic 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 [math], but not so when $d\geq 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

[TABLE]

with jump rates

[TABLE]

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

[TABLE]

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\geq 1$ is governed by

[TABLE]

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

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

[TABLE]

and the dual of the differential operator $\frac{\partial}{\partial x}\left(\frac{\partial}{\partial x}-\frac{d-1}{x}\right)$ 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 (197a,197b) for $d=3$. Making the Laplace transform with respect to time and taking a generating function

[TABLE]

leads to

[TABLE]

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

[TABLE]

and

[TABLE]

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 time scale 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 time scale and the amplitudes are different.

C.2 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_{k}$ to die and leave $k$ new individuals for $n=0,2,3,\cdots$. 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

[TABLE]

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.

C.3 More general models

We consider a more general Markov model for multiplication and diffusion. Models with several fertile sites were studied for instance in Carmona ; Bul18 ; Bul11 , 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

[TABLE]

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

[TABLE]

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 $1$-time information is the generating function

[TABLE]

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

[TABLE]

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

[TABLE]

with initial conditions z_{\bullet}(0)=z_{\mathbin{\vbox{\hbox{\scalebox{0.5}{\bullet}}}}}, the generating function is \mathfrak{Z}(z_{\mathbin{\vbox{\hbox{\scalebox{0.5}{\bullet}}}}},t)=\mathfrak{Z}(z_{\mathbin{\vbox{\hbox{\scalebox{0.5}{\bullet}}}}}(t),0). Solving (208) 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 (81) with the substitution $z=e^{\lambda}$.

C.4 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 [math]. 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 [math] 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 [math]. These probabilities are characteristics of the diffusion on $A$ and do not involve the offspring function. For the site [math], set

[TABLE]

The computation of the $R_{j}$s is complicated in general. As an example when the result is simple, the model with jump rates (195) for $d=3$ leads to

[TABLE]

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

[TABLE]

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, (211) 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]$, $\geq 0$ at [math] and $\leq 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 analogous function $\Pi_{j}(w)$ satisfies

[TABLE]

by the Markov property again. Then

[TABLE]

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 explicitly, this infinite product is not an elementary function.

If $R=R_{0}=0$, i.e. if a walker leaving [math] 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 assume 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 [math]-dimensional analysis: $E(1)\Pi_{0}(w)=E(\Pi_{0}(w))$ is the familiar equation from birth-death processes.

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

[TABLE]

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^{\prime}_{0}(1)=\infty$ (the divergence of a higher derivative while $\Pi^{\prime}_{0}(1)$ remains finite would indicate a multi-critical point). Taking $w\to 1^{-}$ in the derivative of (211),

[TABLE]

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

[TABLE]

For criticality conditions when the walkers hop on a lattice, see Vat1 ; Vat2 ; Vat3 ; Bul11 . For the models studied in the main text, we recover the well-known interpretation of the Watson integral in $d\geq 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 hyper-cubic lattice. Finally, $\mu_{c}=1$ for the model with jump rates (195) is also recovered correctly as $R=1/2$ in that case.

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

[TABLE]

Then the Markov property implies that the process

[TABLE]

is what is called in probability theory a closed martingale (see e.g. Doob ; JP:book ; Williams ), i.e., a quantity conserved on average and converging sample by sample at large times—not only is the expectation time independent

[TABLE]

but even

[TABLE]

In fact,

[TABLE]

and it is instructive (if tedious) to check that the time independence of $\left<U(t,w)\right>$ is also a consequence of (207).

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 [math] (or at any $j\in A$ because $R_{j}=0$ for $j\in A$ by assumption). In the identity

[TABLE]

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 [math] or $\infty$ and $w^{N(\infty)}=\mathbf{1}_{N(\infty)=0}$. Thus (219) implies

[TABLE]

The martingale property is robust: under mild assumption, (223) 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

[TABLE]

Taking the average

[TABLE]

Now

[TABLE]

while

[TABLE]

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 (225) and rearranging gives

[TABLE]

which has a scale invariant limit in the critical limit when the extinction probability goes to [math], namely

[TABLE]