# Clusters in the critical branching Brownian motion

Benoît Ferté  Laboratoire de Physique de l’École normale supérieure, ENS,
Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005
Paris, France  Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France
Pierre Le Doussal  Laboratoire de Physique de l’École normale supérieure, ENS,
Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005
Paris, France  Alberto Rosso  Université Paris-Saclay, CNRS, LPTMS, 91405,
Orsay, France  Xiangyu Cao  Laboratoire de Physique de l’École normale
supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris
Cité, F-75005 Paris, France

###### Abstract

Brownian particles that are replicated and annihilated at equal rate have
strongly correlated positions, forming a few compact clusters separated by
large gaps. We characterize the distribution of the particles at a given time,
using a definition of clusters in terms a coarse-graining length recently
introduced by some of us. We show that, in a non-extinct realization, the
average number of clusters grows as $\sim t^{D_{\mathrm{f}}/2}$ where
$D_{\mathrm{f}}\approx 0.22$ is the Haussdoff dimension of the boundary of the
super-Brownian motion, found by Mueller, Mytnik, and Perkins. We also compute
the distribution of gaps between consecutive particles. We find two regimes
separated by the characteristic length scale $\ell=\sqrt{D/\beta}$ where $D$
is the diffusion constant and $\beta$ the branching rate. The average number
of gaps greater than $g$ decays as $\sim g^{D_{\mathrm{f}}-2}$ for $g\ll\ell$
and $\sim g^{-D_{\mathrm{f}}}$ for $g\gg\ell$. Finally, conditioned on the
number of particles $n$, the above distributions are valid for $g\ll\sqrt{n}$;
the average number of gaps greater than $g\gg\sqrt{n}$ is much less than one,
and decays as $\simeq 4(g/\sqrt{n})^{-2}$, in agreement with the universal gap
distribution predicted by Ramola, Majumdar, and Schehr. Our results
interpolate between a dense super-Brownian motion regime and a large-gap
regime, unifying two previously independent approaches.

##  I Introduction

The critical branching Brownian motion (BBM) is a simple diffusion-reaction
model: non-interacting Brownian particles have the same diffusive constant
$D$, and can replicate themselves and be annihilated with an equal rate
$\beta$; initially, there is one single particle at the origin. It is the
critical point separating two dynamical phases: the sub-critical phase where
the annihilation rate is greater and the process rapidly goes extinct; and the
super-critical phase with a greater branching rate and an exponentially
growing number of particles. At criticality, the number of particles is
conserved on average, but displays strong statistical fluctuations that grow
with time. The process survives up to time $t$ with probability $\propto
1/(\beta t)$; yet, if it survives, the typical number of particles is large
$n_{\text{typ}}\propto\beta t$.

Figure 1: A realization of a discrete analogue of the critical branching
Brownian motion (simulated up to $1024$ time steps). In each time step, each
particle replicates itself with probability $1/2$, and is removed otherwise.
Then, each particle displaces by an independent standard Gaussian variable.
The particles form dense clusters that can be separated by large gaps.

The positions of these particles are also strongly correlated at criticality.
As particle creation takes place only where they already exist, they tend to
agglomerate and form clusters. This clustering property (sometimes referred to
as patchiness) makes the critical BBM model an interesting model in various
domains of physics and biology, such as population dynamics [1, 2, 3, 4, 5],
epidemics [6, 7, 8], genetics [9, 10], and neutron scattering in nuclear
reactors [11, 12, 13]. While its applications often concern higher dimensions,
clustering is the strongest in one dimension. Figure 1 shows a realization of
the 1D BBM in spacetime. The herd of particles form a few dense clusters. They
ebb and flow collectively, and are separated by large gaps from time to time.
As a first step to quantitatively understand the clustering dynamics, we may
seek to characterize the spatial distribution of the particles at a fixed
time. One would like to define and count the clusters, and characterize the
size distribution of the gaps between them. These are nontrivial questions,
and have been recently pursued by two approaches.

The first one [14, 15] focuses on a natural observable: the gaps between
consecutive particles. More precisely, the authors considered realizations
with a fixed number $n$ of particles alive at time $t$, and computed exactly
the probability distribution of the gap between the $k$-th and $(k+1)$-th
consecutive particles, $k=1,\dots,n-1$, using a mapping to a hierarchy of
Kolmogorov–Petrovsky–Piskunov [16] (KPP) equations. They found that this
distribution has a well-defined limit as $t\to\infty$. Moreover, the limit
probability density function $P_{n,k}(g)$ has a universal tail,

$P_{n,k}(g)\simeq\frac{8D}{\beta g^{3}}\,,\,g\to\infty\,,$ (1)

where the prefactor is exact and the same for any $n$ and $k$. Here and below,
$A\simeq B$ means that $A/B$ tends to one in the regime specified; $A\sim B$
means that $A/B$ tends to an unknown constant.

The second approach, mainly pursued in the probability theory literature,
considers a continuum limit of the critical BBM, called the super-Brownian
motion, where the branching rate is so large that it is more suitable to
consider the system as a random fluid 111Usually, one also modifies the
initial condition and consider a large number of particles at $x=0$ instead of
one. However, if we condition on non-extinction, these two initial conditions
are equivalent. See Appendix B.. It is then described by a time dependent
random density profile $\rho\equiv\rho(x,t)$, which can be quite irregular,
and whose time evolution is governed by a stochastic differential equation
[17, 18]. It was pointed out recently [19] that this continuum limit of the
BBM is equivalent to the Brownian force model, a mean-field theory of the
avalanches which occur near the depinning transition of an elastic interface
in a disordered medium [20, 21, 22]; in that context, the local density
corresponds to the local velocity of the interface. The notion of clusters is
well-defined in this continuum limit (see below). Indeed, $\rho$ is positive
on a number of intervals, separated by gaps where $\rho=0$. We can view these
intervals as the clusters of the SBM. We can count them by measuring the size
of their boundary (since every intervals has two boundary points). It was
shown in [23] that in this continuum limit the boundary has a nontrivial
fractal dimension, whose value $D_{\mathrm{f}}$ is determined in terms of the
leading eigenvalue of a particular Ornstein–Uhlenbeck generator with a killing
term. Bounds were provided in [23] and here we determine
$D_{\mathrm{f}}\approx 0.22$ by solving numerically this Ornstein–Uhlenbeck
problem (see (27) below). We can interpret this result as follows. The extent
spanned by the particles grow diffusively as $\propto\sqrt{t}$, hence the
average number of clusters in a non-extinct realization can be heuristically
estimated as follows:

$\left<N_{c}\right>_{\text{non-ex}}\sim
c\left(\sqrt{t}\right)^{D_{\mathrm{f}}}\sim c\,t^{0.11}$ (2)

where $\left<\dots\right>_{\text{non-ex}}$ denotes an average on non-extinct
realizations, and $c$ depends on a short-distance cutoff necessary to make
$N_{c}$ finite [see (6) below]. So, the number of clusters grows as a power
law in time, albeit very slowly. This result raises a number of outstanding
questions concerning the critical BBM model away from the continuum limit. How
can we define the notion of clusters in from a finite number of points? If so,
can we recover (2)? Moreover, does the exponent $D_{\mathrm{f}}$ inform us
about the gaps? Last but not least, the results of the two approaches (2) and
(1) seem completely unrelated. How can we reconcile them in a more complete
picture?

Figure 2: The definition of clusters, illustrated with a realization with $3$
particles (colored for visibility). We surround each of them with an interval
of size $b$. For the value of $b$ illustrated, there are two clusters (red +
blue, and black), i.e. $N_{c}(b)=2$.

##  II Summary of results

###  II.1 Method

In this work, we address these questions using a simple approach recently
introduced by some of us [24]. The basic idea is to consider a coarse-grain
scale $b$, and assign to each particle an interval of size $b$ centered around
it. If two particles are separated by a distance less than $b$, their
respective intervals will overlap and form clusters. See Figure 2 for an
illustration. We shall take this as a definition of clusters and denote the
number of clusters by $N_{c}(b)$, to emphasize its $b$ dependence. It is a
decreasing function of $b$. Indeed, for a fixed realization of $n$ particles,
$N_{c}(b\to 0)=n$ as the intervals are too small to overlap with each other.
As $b$ increases, and exceeds the size of a gap, two neighboring clusters will
merge and $N_{c}(b)$ will decrease by one. Finally, when $b$ is greater than
all the gaps, $N_{c}(b)=1$, as all the intervals are connected. In other
words, $N_{c}(b)$ is one plus the number of gaps larger than $b$ for non-
extinct realizations ($N_{c}(b)=0$ for extinct ones):

$N_{c}(b)=\begin{dcases}\text{number of gaps $\geq b$}&n>0\text{ (non-
extinct)}\\\ 0&n=0\text{ (extinct)}\end{dcases}$ (3)

Therefore its dependence on $b$ is closely related to the distribution of the
gaps. More precisely, in terms of $P_{n,k}(g)$, the probability density
function of the $k$-th gap conditioned on there being $n$ particles, we have

$\left<N_{c}(b)\right>_{n}=1+\sum_{k=1}^{n-1}\int_{b}^{\infty}P_{n,k}(g)\mathrm{d}g\,,$
(4)

where $\left<[...]\right>_{n}$ denotes an average conditioned on the number of
particles at time $t$.

###  II.2 Typical realizations

A first result of this work is to understand the large time $\beta t\gg 1$
behavior of $N_{c}(b)$ of a typical non-extinct realization. Here, we
heuristically derive it with a scaling argument combined with (2) (a more
refined analysis presented later will confirm it). To do this, we identify the
relevant length scales of the problem. They are the following, in decreasing
order:

  1. 1.

the extent spanned by all the particles $\xi\sim\sqrt{Dt}$, i.e., the typical
displacement of the leftmost and rightmost particle of the critical BBM [14].

  2. 2.

the “mean free path” of a particle between consecutive branching/annihilation
events, $\ell=\sqrt{D/\beta}$. This length scale becomes $1$ in the natural
units $D=\beta=1$.

  3. 3.

the inverse density, $a\sim\xi/n_{\text{typ}}\sim\sqrt{D/(\beta^{2}t)}$.

By dimensional analysis, $N_{c}(b)$ must be a function of the two
dimensionless quantities of the problem, $\beta t$ and $b/\ell$. Intuitively,
we expect $N_{c}(b)\sim n_{\text{typ}}\sim\beta t$ for very small gaps
$b\lesssim a$, $N_{c}(b)\sim 1$ for very large gaps $b\gtrsim\xi$; Thus
$N_{c}(b)$ can have the nontrivial $t$ dependence (2) only at the remaining
relevant length scale $b\sim\ell$: $N_{c}(b\sim\ell)\sim(\beta
t)^{D_{\mathrm{f}}/2}$; the $\beta$ dependence here is fixed by the
dimensional analysis. Finally, it is natural to expect that $N_{c}(b)$ behaves
as a power law in each of the scaling regimes $a\ll b\ll\ell$ and $\ell\ll
b\ll\xi$. Gathering the above information, we surmise that:

$\left<N_{c}(b)\right>_{\text{non-
ex}}\sim\begin{dcases}(b/\ell)^{D_{\mathrm{f}}-2}(\beta
t)^{\frac{D_{\mathrm{f}}}{2}}=\beta^{D_{\mathrm{f}}-1}D^{1-\frac{D_{\mathrm{f}}}{2}}\,t^{\frac{D_{\mathrm{f}}}{2}}\,b^{D_{\mathrm{f}}-2}&a\lesssim
b\lesssim\ell\\\ (b/\ell)^{-D_{\mathrm{f}}}(\beta
t)^{\frac{D_{\mathrm{f}}}{2}}=D^{\frac{D_{\mathrm{f}}}{2}}\,t^{\frac{D_{\mathrm{f}}}{2}}\,b^{-D_{\mathrm{f}}}&\ell\lesssim
b\lesssim\xi\end{dcases}\,.$ (5)

In other words, most gaps of a typical realization fall into two categories,
each governed by a different power law distribution. To understand this
intuitively, we may inspect Figure 1 once more in light of the formula (5).
Many small gaps are not visible in that Figure. Eq. (5) predicts that their
sizes are governed by a rather large exponent $N_{c}(b)\propto
b^{D_{\mathrm{f}}-2}\approx b^{-1.78}$. Most such gaps are comparable to the
inverse density $\sim 1/\sqrt{t}$. Meanwhile, the large gaps that are visible
in Figure 1 have a wide distribution, $N_{c}(b)\propto
b^{D_{\mathrm{f}}-2}\approx b^{-0.22}$. These predictions are verified
numerically, see Figure 4 in Section IV below.

We now make connection with the continuum limit of the super-Brownian motion
(SBM) and of the Brownian force model. For this, we take $\beta$ to infinity,
while keeping $D$ of order one in (5) so that the typical interparticle
distance scale $a\to 0$. Our result (5) becomes the following:

$\left<N_{c}(b)\right>_{\text{non-ex}}\sim
D^{\frac{D_{\mathrm{f}}}{2}}\,t^{\frac{D_{\mathrm{f}}}{2}}\,b^{-D_{\mathrm{f}}}\,,\,b\ll\sqrt{Dt}\,\text{
(SBM, or {\hbox{\beta\gg 1}} limit)}$ (6)

Indeed, the small-gap regime of (5) disappears ($\ell\to 0$), and the large-
gap one extends to all values of $b\lesssim\sqrt{Dt}$ and remains well-
defined, being independent of $\beta$. The $t$ dependence of (6) and (2) are
the same, yet (6) has the correct dimension. Its $b$-dependence should be
interpreted as describing the size distribution of the gaps between
neighboring connected components of the support of the density $\rho$. In the
SBM limit, there are infinitely many such gaps, which can be arbitrarily
small. On the other hand, the average number of gaps $\geq b$ is finite, and
decays as $b^{{-D_{\mathrm{f}}}}$. We can view $b$ as the short-distance
cutoff that is necessary to have a finite number of clusters. This cutoff is
different from that introduced in [23], which is a cutoff on the value of the
local density. In that sense the above result is new, and it is interesting to
note that it also applies, as mentioned above, to the gaps in the support of
the velocity in a mean-field avalanche.

###  II.3 Particle number conditioning

The puzzle remains with regard to the first approach [14, 15]: eq. (5) still
looks unrelated to (1). A second result of this work reconciles the two
approaches by calculating the averaged cluster number, conditioned on a fixed
number of particles $n$ (at time $t$). We shall show that, the conditioned
average number of clusters has a long time limit with the following asymptotic
behaviors for $n\gg 1$:

$\left<N_{c}(b)\right>_{n}-1\stackrel{{\scriptstyle
t\to\infty}}{{\sim}}\begin{dcases}c(b/\xi_{n})^{-D_{\mathrm{f}}}&1\ll
b\ll\xi_{n}\\\
4(b/\xi_{n})^{-2}&b\gg\xi_{n}\end{dcases}\,,\,\xi_{n}=(Dn/\beta)^{\frac{1}{2}}\,.$
(7)

Here, the crossover scale $\xi_{n}$ can be understood as the typical extent of
a realization with $n$ particles. The result in the $b\ll\xi_{n}$ regime
involves an unknown factor $c$ that is independent of $b$. Yet, the prefactor
of the $b\gg\xi_{n}$ regime is exact.

Eq. (7) reconciles the result (1) of [14] and (5) just announced. Indeed, the
formula for $b\gg\xi_{n}$ regime is consistent with (1). To see why, recall
that $\left<N_{c}(b)\right>_{n}-1$ is the average number of gaps greater than
$b$. Then, Eq. (1) implies that each of the $n-1$ gaps is larger than $b$ with
probability $4D/(\beta b^{2})$. So the average total number of gaps greater
than $b$ is $4(n-1)D/(\beta b^{2})$ by (4) 222Note that we do not need to
assume independence between gaps, since we are calculating a sum of
expectation values, not a product., which is asymptotically equivalent —
including the prefactor — to (7) as $n\gg 1$. Also, the formula for
$b\ll\xi_{n}$ is consistent with the our result (5) without particle number
conditioning: by setting $n$ to be the typical number of particles alive at
$t$, $n_{\rm typ}\sim\beta t$, we recover the $b\gg 1$ regime of (5).

Last but not least, (7) provides an additional insight on the distinct nature
of the two approaches. It implies in particular
$\left<N_{c}(b=\xi_{n})\right>_{n}\sim 1$, i.e., in a typical realization of
$n\gg 1$ particles, the largest gap size is comparable to $\xi_{n}$. If we
take a single typical configuration with $n\gg 1$ particles and make a
histogram of its $n-1$ gaps, according to (7), it will display a $P(g)\sim
g^{-D_{\mathrm{f}}-1}$ distribution with a cutoff at $g\sim\xi_{n}$.
Meanwhile, if we consider many realizations with $n$ particles, large gaps
$g\gg\xi_{n}$ will appear from time to time. A histogram of the gaps of all
the realizations combined will show a $P(g)\sim g^{-3}$ power law for
$g\gg\xi_{n}$. Observing this numerically is very challenging, especially for
large $n$. Nevertheless, thanks to the knowledge of the exact prefactor, Ref.
[14] was able to identify the $g^{-3}$ regime for $n\leq 10$.

The rest of the paper is organized as follows. In Section III we analyze the
model using a mapping to a KPP equation, and derive the predictions announced
above. Section IV reports a few numerical tests supporting the predictions. We
close with concluding remarks in Section V. The main text is supplemented with
Appendices containing analytical details and numerical methods.

##  III Analysis of the KPP equation

###  III.1 Setting up the mapping

In this section we set up a mapping to a set of KPP equations that allows to
calculate $\left<N_{c}(b)\right>$ and related quantities in the critical BBM
model. From now on we work with the natural units $D=\beta=1$. This amounts to
measure time in the unit of the inverse branching rate $1/\beta$ and length in
the unit of the mean free path $\sqrt{D/\beta}$. Thus, the dependence on $D$
and $\beta$ can be readily restored.

Let $x_{1}(t),\dots,x_{n}(t)$ be the position of the particles; here $n=n(t)$
is the number of the particle, and is a random variable. Consider the
following observable (here and below, $\left<[\dots]\right>$ denotes an
average on all realizations, including those already extinct by time $t$; we
recall that the initial condition is that of a single particle at the origin)
:

$E(x,t):=\left<\prod_{i=1}^{n}f(x-x_{i}(t))\right>\,,$ (8)

$f$ is any function (if $n=0$, the product is equal to $1$ by convention). A
standard backward recursion argument (see Appendix A) shows that $E(x,t)$
satisfies the following equation

$\displaystyle\partial_{t}E=\partial_{x}^{2}E+(E^{2}-E)\,,\,E|_{t=0}=f(x)\,.$
(9)

It is often convenient to rewrite the KPP equation by considering

$F(x,t):=1-E(x,t)\,,$ (10)

in terms of which (9) becomes particularly simple at criticality:

$\partial_{t}F=\partial_{x}^{2}F-F^{2}\,,\,F|_{t=0}=1-f(x)\,.$ (11)

The above discussion applies to any $f$. As a warm-up example, let us consider
$f(x)=e^{-\mu}$, which is independent of $x$. Then $F(x,t)=1-\left<e^{-\mu
n}\right>$ is a moment generating function of the particle number $n=n(t)$.
The KPP equation reduces to an ODE $\dot{F}=-F^{2},F(0)=1-e^{-\mu}$. the
solution is

$F(x,t)=\frac{1}{(1-e^{-\mu})^{-1}+t}:=F_{0}$ (12)

(We gave the expression a name since it will appear a few times.) One can
check that this corresponds to the following well-known distribution of
particle number:

$\mathbb{P}(n(t)=0)=1-\frac{1}{t+1}\,,\,\mathbb{P}(n(t)=k>0)=\frac{t^{k-1}}{(1+t)^{k+1}}\stackrel{{\scriptstyle
t\gg 1}}{{\sim}}\frac{1}{t^{2}}e^{-k/t}\,.$ (13)

In this work, we shall focus on the following choice of $f$ that allows to
make connection with the quantity $N_{c}(b)$:

$f(x)=e^{-\mu}\theta(|x|-b/2)$ (14)

so that

$F(x,t)=1-\left<\prod_{i=1}^{n}\left(e^{-\mu}\theta(|x-x_{i}(t)|-b/2)\right)\right>\,.$
(15)

In plain words, a configuration $x_{1}<\dots<x_{n}$ contributes $1$ to $F$ if
$x$ belongs to at least one of the intervals $[x_{i}-b/2,x_{i}+b/2]$;
otherwise, if $x$ is outside the union of these intervals, the contribution is
$1-e^{-\mu n}$ where $n=n(t)$ is the number of particles. As a consequence,
integrating $F$ over a large interval $|x|<L$ gives

$\int_{-L}^{L}F(x)\mathrm{d}x=2LF_{0}+\left<e^{-\mu n}\ell(b)\right>\,.$ (16)

where $F_{0}=1-\left<e^{-\mu n}\right>$ is given above (12). We have
introduced another geometric quantity of interest, $\ell(b)$, called the
extension, defined as the total length of the union of the intervals (the
overlaps being counted only once). It is not hard to see that

$F(|x|\to\infty)\to F_{0}=\frac{1}{(1-e^{-\mu})^{-1}+t}\,.$ (17)

To get rid of diverging part $2LF_{0}$, we can subtract by $F$ by $F_{0}$ and
take $L\to+\infty$ :

$\left<\ell(b)e^{-\mu
n}\right>=\int_{-\infty}^{+\infty}(F(x,t)-F_{0})\mathrm{d}x\,.$ (18)

How does this relate the number of clusters? Thanks to a geometric relation:
$N_{c}(b)$ can be obtained by deriving $\ell(b)$ with respect to $b$:

$N_{c}(b)=\partial_{b}\ell(b)\,.$ (19)

To see why this is so, consider the change of $\ell(b)$ as $b$ increases by an
infinitesimal amount $\mathrm{d}b$. The contributions come from the two
extremities of each cluster, which expand by $\mathrm{d}b/2$, therefore,
$\mathrm{d}\ell(b)=2N_{c}(b)\mathrm{d}b/2=N_{c}(b)\mathrm{d}b$. This is
exactly (19).

Combining (19) and (18), we obtain

$\left<N_{c}(b)e^{-\mu
n}\right>=\int_{-\infty}^{+\infty}\partial_{b}F(x,t)\mathrm{d}x\,.$ (20)

The equations (18) and (20) allow us to extract geometric information from the
solution of the KPP equation. The presence of $e^{-\mu n}$ will allow us to
obtain $\left<N_{c}(b)\right>_{n}$ the mean cluster number conditioned on the
number of particles, by an inverse Laplace transform and using (13). More
precisely, for any observable $\mathcal{O}$, if its average conditioned on
particle number $n$ depends on $n$ as a power law,

$\left<\mathcal{O}\right>_{n}\simeq n^{a}\,,\,n\gg 1\,,$ (21)

using (13), we have

$\left<\mathcal{O}e^{-\mu
n}\right>\simeq\frac{1}{t^{2}}\int_{0}^{\infty}e^{-\mu
n-n/t}n^{a}\mathrm{d}n\simeq\frac{\Gamma(a+1)}{t^{2}\mu^{a+1}}\,,\,1\ll\mu^{-1}\ll
t\,.$ (22)

The analysis of the KPP equation will provide a result in the form of (22),
which we can readily translate to (21).

###  III.2 Overview of the analysis

In the rest of this section, we will analyze the following KPP equation:

$\partial_{t}F=\partial_{x}^{2}F-F^{2}\,,\,F|_{t=0}=\begin{cases}1&|x|\leq
b/2\\\ 1-e^{-\mu}&|x|>b/2\end{cases}$ (23)

The initial condition will be often referred to as a “plateau” of width $b$.
Our goal is to establish the asymptotic behavior of the solution in a number
of regimes, in order to derive the results announced in the introduction. We
will proceed in two main steps. First, in section III.3, we will set $\mu=0$,
and consider $1\ll b\ll\sqrt{t}$ (typical large gap regime) as well as $1\ll
b^{-1}\ll\sqrt{t}$ (typical small gap regime). This analysis will lead to the
result (5). Technically, it will crucially involve the Onstein-Uhlenbeck
generator introduced in [23]. Next, in section III.4, we will set $\mu>0$, and
consider various regimes where $b\gg 1$, $t\gg 1$ and $\mu\gg 1$, with the
goal of deriving (7).

###  III.3 Typical regimes ($\mu=0$)

In this section, we shall consider the KPP equation (23) with $\mu=0$. Then,
the initial condition $F|_{t=0}=\theta(b/2-|x|)$ becomes zero outside the
plateau: it resembles a Dirac delta function viewed from a scale $\gg b$.
Indeed, we find numerically that the solution with a plateau initial condition
has a similar long-time behavior as that with a delta initial condition, which
has been well studied in the mathematical literature [25, 23]. We will first
review the latter, and then discuss how to adapt the theory to our initial
condition.

####  III.3.1 Delta initial condition: a review

Consider the KPP equation with a delta function initial condition,

$\partial_{t}F=\partial_{x}^{2}F-F^{2}\,,\,\lim_{t\to 0}F=\lambda\delta(x)\,.$
(24)

This equation has no interpretation in BBM but appears naturally in the SBM
limit [23] and we will follow here the approach of that work. A main result
thereof is that the solution has the following long-time asymptotic Ansatz:

$F(x,t)=\frac{1}{t}\left[f(x/\sqrt{t})+Ct^{-\eta}\varphi(x/\sqrt{t})+\dots\right]\,.$
(25)

Let us explain this equation in detail.

The leading term $t^{-1}f(x/\sqrt{t})$ should be a solution to the KPP
equation (23) itself. This imposes the following equation on $f(y)$:

$f^{\prime\prime}+\frac{y}{2}f^{\prime}-f^{2}+f=0\,.$ (26)

This ODE has a unique nonzero solution that decays to $0$ as $y\to\pm\infty$.
It can be also specified by the initial conditions $f^{\prime}(0)=0$ and
$f(0)=0.6898\dots$. It resembles a Gaussian, and has a tail
$f(y)\sim|y|\exp(-y^{2}/4)$ as $|y|\to\infty$ (see Figure 3 for a plot).

Figure 3: The functions $f$ and $\varphi$ in the Ansatz (25). They are
obtained by solving (26) and (27), respectively, using a standard numerical
shooting method. The interpretation of these functions in the BBM is discussed
below after (38).

The subleading term satisfies a linearisation of the KPP equation around the
leading term. One can express this as a an eigenvalue equation involving an
operator $H$ acting on the function $\varphi(y)$:

$\eta\varphi=H\varphi\,,\,H=-\partial_{y}^{2}-\frac{y}{2}\partial_{y}+2f(y)-1\,.$
(27)

Although $H$ is not Hermitian, it is a generator of an Ornstein–Uhlenbeck
process with a killing term and can be transformed into a Schrödinger
Hamiltonian:

${H}_{S}:=e^{y^{2}/8}He^{-y^{2}/8}=-\partial_{y}^{2}+\frac{y^{2}}{16}+2f(y)-\frac{3}{4}\,.$
(28)

Therefore, $H$ has all the spectral properties of a standard Schrödinger
Hamiltonian. Because the potential is confining, it has a discrete spectrum.
We determine $\eta$ and $\varphi(y)$ as the ground state energy and eigen-
function, respectively 333Note that there is another eigenfunction of $H$,
with negative energy. It tends to a nonzero constant as $y\to\infty$, and maps
to a non-normalizable wave-function of $H_{S}$. That eigenfunction corresponds
to an instability, and will be relevant in the crossover to a $\mu>0$
solution, see Appendix C below. Here, we shall restrict to decaying
eigenfunctions of $H$ that correspond to normalizable wavefunctions of
$H_{S}$.. $\varphi(y)$ has a Gaussian decay $\sim e^{-y^{2}/4}$ as
$|y|\to\infty$ (which is determined by the diffusion term of the KPP
equation). We fix the overall normalization of $\varphi$ by requiring that
$\varphi(0)=1$. The value of $\eta$ is not known in a closed form; we
estimated it numerically as equal to

$\eta\approx 0.3904\dots\,.$ (29)

We note that $\eta$ is related to the exponent $\lambda_{0}$ in [23], and to
the fractal dimension $D_{\mathrm{f}}$ by the following ([23], Theorem 1.3):

$\eta=\lambda_{0}-1/2\,,\,D_{\mathrm{f}}=2-2\lambda_{0}=1-2\eta\approx
0.22\,,$ (30)

as quoted in the introduction. The eigen-function $\varphi(y)$ is plotted in
Figure 3. Its shape does not resemble a Gaussian, and has a pair of maxima
away from $0$.

The linear Schrödinger equation (27) alone cannot fix the prefactor $C$ of the
subleading term in (25). Instead, $C$ must be determined by the initial
condition. With $\delta$-function initial conditions, this can be done by a
scaling argument. Indeed, the rescaling transformation

$t=\alpha^{2}\tilde{t}\,,\,x=\alpha\tilde{x}\,,\,F=\alpha^{-2}\tilde{F}$ (31)

leaves the KPP equation $\partial_{t}F=\partial_{x}^{2}F-F^{2}$ invariant,
i.e., one obtains the same equation with tildes. The rescaling transforms the
initial condition in (24) as follows: $\lambda=\alpha^{-1}\tilde{\lambda}$.
Also, the Ansatz (25) transforms to itself, which implies
$C=\alpha^{2\eta}\tilde{C}$. The above two equations together impose the
following scaling law:

$F_{t=0}=\lambda\delta(x)\implies C\propto-\lambda^{-2\eta}$ (32)

Note that when $\lambda\to\infty$, $C\to 0$. The minus sign of (32) must be
there since $C$ must increase with $\lambda$ (this is because the KPP equation
is monotonous: if $F(x,t=0)\geq F^{\prime}(x,t=0)$ for all $x$, then
$F(x,t)\geq F^{\prime}(x,t)$ for all $t$ and $x$).

####  III.3.2 Plateau initial condition

Having reviewed in some detail the solution to the KPP equation with a delta
initial condition, we come back to our plateau initial condition of width $b$.
We expect the solution to obey the Ansatz (25) in the long time limit, with
however a different prefactor $C$, which now depends on $b$. We now determine
that dependence, as well as the time scale above which the Ansatz (25) is
valid, for small and large $b$.

For small gaps, $b\ll 1$, the plateau is indistinguishable from a delta peak
with $\lambda=b$. Therefore

$C\sim-b^{-2\eta}\,,\,b\ll 1\,.$ (33)

We now argue that the asymptotic form (25) is valid only when $t\gg b^{-2}$
for small $b$. Indeed, the initial condition is very singular, so the
diffusive term will dominate (against the nonlinear term) initially, leading
to a Gaussian solution with $F(0,t)\sim b/\sqrt{t}$. Matching that with the
behavior $F(0,t)\sim 1/t$ in (25) gives us a crossover time $t\sim b^{-2}$,
beyond which (25) is valid.

The case of large gaps, $b\gg 1$, is different. Indeed, the asymptotic Ansatz
(25) holds when and only when $t\gtrsim b^{2}$; at earlier times, the solution
$F(x,t)$ remains almost independent of $x$ away from the boundaries of the
plateau. This changes when $t=O(b^{2})$ the time for diffusion over a distance
of order the size of the plateau.. Having identified this time scale, we can
determine $C$ by a scaling argument: the relative amplitude of the subleading
term, $Ct^{-\eta}$, should be of order unity at the only relevant time scale
$t\sim b^{2}$. This fixes

$C\sim b^{2\eta}\,,\,b\gg 1\,.$ (34)

The sign must be positive here since $C$ must increase with $b$.

In summary, the long time asymptotic solution to the KPP equation is given by
(25) with the prefactor $C$ satisfying (33) and (34):

$F(x,t)=\frac{1}{t}\left[f(x/\sqrt{t})+Ct^{-\eta}\varphi(x/\sqrt{t})+\dots\right]\,,\,C=C(b)\sim\begin{dcases}-b^{-2\eta}&1\gg
b\gg 1/\sqrt{t}\\\ b^{2\eta}&1\ll b\ll\sqrt{t}\end{dcases}$ (35)

Applying (18) (with $\mu=0$ and $F_{0}=0$) to this equation we obtain the
average extension

$\left<\ell(b)\right>\sim
c_{0}t^{-\frac{1}{2}}+C(b)t^{-\frac{1}{2}-\eta}+\dots\,,\quad 1/\sqrt{t}\ll
b\ll\sqrt{t}$ (36)

where $c_{0}=\int_{-\infty}^{+\infty}f(y)\mathrm{d}y=2.913\dots$. Now, the
mean cluster number is given by deriving $\left<\ell(b)\right>$ with respect
to $b$ (19), so only the subleading term $\propto C(b)t^{-\frac{1}{2}-\eta}$
contributes, and we obtain the following:

$\left<N_{c}(b)\right>\sim\begin{dcases}b^{2\eta-1}t^{-\frac{1}{2}-\eta}&1\ll
b\ll\sqrt{t}\\\ b^{-2\eta-1}t^{-\frac{1}{2}-\eta}&1\gg b\gg
1/\sqrt{t}\end{dcases}\,.$ (37)

Since the above averages include extinct realizations, it is helpful to
restrict to the non-extinct ones, recalling that the BBM survives with
probability $\sim 1/t$ (13):

$\boxed{\left<\ell(b)\right>_{\text{non-ex}}\simeq
c_{0}t^{\frac{1}{2}}\,,\,c_{0}=2.913\dots,\frac{1}{\sqrt{t}}\ll
b\ll\sqrt{t}\,,\;\left<N_{c}(b)\right>_{\text{non-
ex}}\sim\begin{dcases}b^{2\eta-1}t^{\frac{1}{2}-\eta}&1\ll b\ll\sqrt{t}\\\
b^{-2\eta-1}t^{\frac{1}{2}-\eta}&\frac{1}{\sqrt{t}}\ll b\ll 1\end{dcases}\,.}$
(38)

This is the main result of this section. The part of it concerning $N_{c}(b)$
is announced in (5) in terms of the fractal dimension $D_{\mathrm{f}}=1-2\eta$
(30), where we also restored the dependence on $\beta$ and $D$ by requiring
that $N_{c}(b)$ be dimensionless.

A striking point of (38) is that the leading asymptotic of
$\left<\ell(b)\right>_{\text{non-ex}}\sim t^{1/2}$ is independent of $b$,
including the prefactor. This means that the particles form a rather compact
bulk within the diffusive length scale $|x|\lesssim\sqrt{t}$, densely
populating a finite portion of that length. Changing the coarse-graining scale
affects that portion only by an infinitesimal amount in the $t\to\infty$
limit. Furthermore, the average number of clusters grows qualitatively more
slowly than the extension, indicating the existence of large clusters. All of
the above is reminiscent of the spatial distribution of all the positions
visited by a critical branching fractional Brownian motion in certain
dimensions [24].

We remark that $tF(x,t)$ is the probability, conditioned on non extinction,
that there is at least one particle in $[x-b/2,x+b/2]$. Therefore, (35)
implies that, when $t$ is large, this probability is given by $f(x/\sqrt{t})$,
plotted in Figure 3. Similarly, $t\partial_{b}F(x,t)$ is the probability
density (conditioned on non extinction) that the closest particle to $x$ is at
distance $b/2$ from $x$. For $t$ large, this probability is proportional to
$\partial_{b}C(b)t^{-\eta}\varphi(x/\sqrt{t})$, where $C(b)$ has the two
limiting behaviors given in (35). The shape of $\varphi(y)$ (Figure 3) with
two maxima away from $y=0$ means that it is more probable to have gaps away
from the origin (which is where the BBM started).

###  III.4 Conditioned regimes ($\mu>0$)

We now turn to study the averages with particle number conditioning, by
setting $\mu>0$ but $\mu\ll 1$. We also focus uniquely on the large gaps $b\gg
1$. There are thus two large time scales:

  1. 1.

$b^{2}$ is the characteristic time for a Brownian particle to displace by a
distance $b$; thus, gaps of size $b$ can probably appear when and only when
$t\gtrsim b^{2}$. In terms of the KPP equation, $b^{2}$ is the time it takes
for the plateau to “melt” under diffusion.

  2. 2.

$\mu^{-1}$ is a time scale set by the conditioning on the particle number.
Indeed, the effect of the chemical potential $\mu$ is to suppress the
contribution of the realizations with $n\gg\mu^{-1}$. Thus, by (13), for
$t\ll\mu^{-1}$, there are virtually no such realizations, so $\mu$ has no
effect; for $t\gg\mu^{-1}$, the typical realizations (with particle number
$\gg\mu^{-1}$) are suppressed, and we are favoring realizations with
atypically few particles.

In terms of the KPP equation, $\mu^{-1}$ is the time scale of the plateau
melting under the effect of the nonlinear term. Indeed, If we solve the KPP
equation (23) approximately by ignoring the diffusive term, we get an ODE
$\dot{F}=-F^{2}$ which has a solution

$F(x,t)=\begin{dcases}\frac{1}{1+t}&|x|<b\\\
F_{0}\approx\frac{1}{\mu^{-1}+t}&|x|>b\end{dcases}\,.$ (39)

Here we also approximated $(1-e^{-\mu})\approx\mu$ as $\mu\ll 1$. The ratio
between the two values tends to $1$ when $t\gg\mu^{-1}$.

From the above discussion, it is reasonable to expect that the asymptotic
behavior of the solution to the KPP equation depends crucially on whether
$\mu^{-1}$ or $b^{2}$ is greater; We will call the regime where
$b^{2}\ll\mu^{-1}$ ($b^{2}\gg\mu^{-1}$) the typical gap regime (rare gap
regime, respectively). The analysis of the typical gap regime (Section
III.4.1) relies on the preceding $\mu=0$ results; meanwhile, the rare gap
regime (Section III.4.2) calls for a distinct approach, somehow reminiscent of
that in [14].

####  III.4.1 The typical gap regime

In the typical gap regime where $b^{2}\ll\mu^{-1}$, the plateau of the initial
condition (23) will first melt by diffusion, when $t\sim b^{2}$. After that,
there are two time regimes. At intermediate time, $b^{2}\ll t\ll\mu^{-1}$, the
solution can be viewed, for all practical purposes, as identical to the one
with $\mu=0$. That has been studied above, and we found the following
asymptotic form:

$F(x,t)=\frac{1}{t}\left[f(x/\sqrt{t})+Ct^{-\eta}\varphi(x/\sqrt{t})+\dots\right]\,,\,C\sim
b^{2\eta}\,,\,b^{2}\ll t\ll\mu^{-1}\,.$ (40)

Now, when $t\gg\mu^{-1}$, the nonzero value of $\mu$ can no longer be ignored,
and dictates the large distance asymptotics of the solution
$F(x\to\infty)=F_{0}\sim 1/t$. Thus, the Ansatz above with $f(|y|\to\infty)\to
0$ is invalid.

Nevertheless, a similar Ansatz with a different leading order scaling function
$f(y)\to f_{1}(y)=1$ becomes valid. Indeed, $f_{1}(y)=1$ is a solution to
(26). Then, the eigenvalue equation (27) becomes equivalent to that of an
oscillator: the ground state energy is $\eta\to\eta_{1}=3/2$, and the ground
state wavefunction $\varphi(y)\to\varphi_{1}(y)=e^{-y^{2}/4}$. However, it is
important to note that $\varphi_{1}(y)$ is the dominant (with the slowest
decay) perturbation that vanishes as $|y|\to\infty$. There is one more
perturbation, which is constant in $y$: $\varphi_{0}(y)=1$, and with
$\eta_{0}=1$. In terms of the Schrödinger Hamiltonian, it corresponds to an
un-normalizable wavefunction. Such a perturbation was not allowed in the
$\mu=0$ above since $F$ has to vanish at infinity there. When $\mu>0$, it must
exist, and its amplitude is determined by the large-$x$ behavior $F$:

$F(x\to\infty)=F_{0}=\frac{1}{\mu^{-1}+t}=\frac{1}{t}(1-\mu^{-1}t^{-1}+\dots)\,,$

see (17) above.

Gathering the above, we have

$F(x,t)=\frac{1}{t}\left[1-\mu^{-1}t^{-1}+C_{1}t^{-\frac{3}{2}}e^{-x^{2}/{4t}}+\dots\right]\,,\,b^{2}\ll\mu^{-1}\ll
t$ (41)

Again, we need to determine the prefactor $C_{1}$. In principle, we need to
characterize the crossover between the two asymptotic Ansätze; this is done in
Appendix C. However, the result of that analysis can be mostly recovered by a
simple argument, which consists in matching the above two Ansätzes with $x=0$
at the crossover time $t=\mu^{-1}$:

$C_{1}\mu^{3/2}\sim f(0)+c_{1}b^{2\eta}\mu^{\eta}$ (42)

This fixes $C_{1}$ and gives

$F(x,t)=\frac{1}{t}\left[1-\mu^{-1}t^{-1}+\left(c_{0}^{\prime}(t\mu)^{-\frac{3}{2}}+c_{1}t^{-\frac{3}{2}}b^{2\eta}\mu^{\eta-\frac{3}{2}}\right)e^{-\frac{x^{2}}{4t}}+\dots\right]$
(43)

where $c_{0}^{\prime}\approx f(0)$ 444This is a rough approximation. A
numerically exact estimate is $c_{0}^{\prime}\approx 1$, see appendix C.  and
$c_{1}$ is a order one constant independent of $b,t$ and $\mu$. Plugging that
into (18) and (20) we have

$\left<\ell(b)e^{-\mu n}\right>\sim
t^{-2}\mu^{-\frac{3}{2}}\,,\,\left<N_{c}(b)e^{-\mu n}\right>\sim
t^{-2}\mu^{\eta-\frac{3}{2}}b^{2\eta-1}\,.$ (44)

By comparing to (21) and (22), we see that this corresponds to the following
average conditioned on particle number:

$\boxed{\left<\ell(b)\right>_{n}\sim
n^{\frac{1}{2}}\,,\,\left<N_{c}(b)\right>_{n}\sim
n^{\frac{1}{2}-\eta}b^{2\eta-1}\,,\,1\ll b\ll\sqrt{n}\ll\sqrt{t}\,.}$ (45)

Eq. (45) is the $b\ll\sqrt{n}$ case of the result (7).

It is interesting to remark that (45) is identical to (38) upon replacing $n$
by $t$. This can be understood as follows: when the particle number $n\ll t$
is atypically small, the particles at $t$ have a recent common ancestor at
time $t^{\prime}=(t-n)$. As far as $\ell(b)$ and $N_{c}(b)$ are concerned, the
effective “age” of the BBM is $n$ instead of ${t}$. Thus, the cluster
structure discussed below (38) prevail as we condition on a large particle
number $n$.

####  III.4.2 The rare gap regime

In the rare gap regime, where $b^{2}\gg\mu^{-1}$, the plateau of the initial
condition in (23) will melt first from the action of the nonlinear term in the
KPP equation. Indeed, (39) is a good approximation at $t\ll b^{2}$ almost
everywhere 555except near the edges of the plateau. Indeed, the edge is
softened by the diffusion term. It has a width $\sim\sqrt{t}$, much smaller
than $b$.. As a consequence, as $t$ exceeds $\mu^{-1}$, $F(x,t)$ becomes
almost a constant. Indeed, the ratio between the values $F$ inside and far
from the plateau is

$\frac{F(|x|\ll b)}{F_{0}}=\frac{t+\mu^{-1}}{t+1}\approx 1+\frac{1}{t\mu}\,.$
(46)

The last approximation is valid for $t\gg 1,\mu^{-1}\gg 1$. Thus, we can fix a
time $t_{0}$ such that $t_{0}\ll b^{2}$ (so the approximation (39) is still
valid) and that $1/(\mu t_{0})$ is small. At later time, the solution can only
become closer and closer to the constant $F_{0}$, since both terms
$\partial_{x}^{2}F$ and $-F^{2}$ make the solution more uniform. Therefore, we
can treat the KPP equation for $t\geq t_{0}$ by a perturbation expansion
around the constant solution $F_{0}=1/(\mu^{-1}+t)$ (12). The expansion is
controlled by the small parameter $1/(\mu t_{0})$. We write

$F=F_{0}+aF_{1}+a^{2}F_{2}+\dots\,,\,\text{where }a=(\mu t_{0})^{-1}$ (47)

is the small parameter. As we explained above, the zero-th order solution is
the constant $F_{0}=1/(\mu^{-1}+t)$ (12). Plugging (47) into (23) and
comparing order by order in $a$, we find

$\displaystyle(\partial_{t}-\partial_{x}^{2}+2F_{0})F_{1}=0\,,\,F_{1}|_{t=t_{0}}=F_{0}\,\theta(b/2-|x|)=\frac{1}{\mu^{-1}+t_{0}}\theta(b/2-|x|)\,,$
(48)
$\displaystyle(\partial_{t}-\partial_{x}^{2}+2F_{0})F_{2}=-F_{1}^{2}\,,\,F_{2}|_{t=t_{0}}=0\,,$
(49)

and so on. Here, the initial condition of $F_{1}$ comes from (46), which
implies $F(|x|\ll b,t=t_{0})\approx F_{0}(1+a)$ (while outside the plateau,
$F(t=t_{0})\approx F_{0}$). We will then extract the extension and gap
distribution order by order:

$\displaystyle\left<\ell(b)e^{-\mu
n}\right>=\ell^{(1)}+\ell^{(2)}+\dots\,,\quad\ell^{(m)}:=a^{m}\int
F_{m}(x,t)\mathrm{d}x\,,$ (50) $\displaystyle\left<N_{c}(b)e^{-\mu
n}\right>=N_{c}^{(1)}+N_{c}^{(2)}+\dots\,,\quad
N_{c}^{(m)}:=\partial_{b}\ell^{(m)}\,,$ (51)

To carry out the perturbation calculations, it is convenient to introduce a
retarded propagator $G$ defined by the following

$(\partial_{t}-\partial_{xx}+2F_{0})G(x,t|y,s)=\delta(x-y)\delta(t-s)\,,\,G|_{t<s}=0\,.$
(52)

whose solution can be explicitly found

$G(x,t|y,s)=p(x-y,t-s)\frac{(\mu^{-1}+s)^{2}}{(\mu^{-1}+t)^{2}}\,,\,p(z,u)=\frac{1}{\sqrt{4\pi
u}}e^{-\frac{z^{2}}{4u}}\,.$ (53)

Here $p$ is the usual diffusive kernel. Note that the Green function is
translation invariant in space but not in time.

Using the Green function, we can readily solve the perturbative KPPs, (48) and
(49). At first order, we have

$F_{1}(x,t)=\int
G(x,t|y,t_{0})F_{1}(y,t_{0})\mathrm{d}y=\frac{(\mu^{-1}+t_{0})}{(\mu^{-1}+t)^{2}}\int_{-b/2}^{b/2}p(x-y,t-t_{0})\mathrm{d}y\,.$
(54)

Plugging this into (50) and (51), the integral over $p$ being unity, we find
rather trivial geometrical information:

$\ell^{(1)}\sim\frac{b}{\mu t^{2}}\,,\,N_{c}^{(1)}\sim\frac{1}{\mu t^{2}}\,.$
(55)

These can understood as follows. The coarse-grain scale $b$ is so large that a
realization will either have a single cluster of length $b$ or no clusters (if
it is extinct). So a realization with $n>0$ particles contributes $e^{-\mu n}$
(and $e^{-\mu n}b$) to $N_{c}^{(1)}$ (and $\ell^{(1)}$, respectively); extinct
ones do not contribute. We can check that (55) follows from (13) in the limit
$t\gg\mu^{-1}\gg 1$: indeed, for the cluster number, (13) implies
$\left<e^{-\mu n}\mathbf{1}_{n>0}\right>\approx t^{-2}\int_{0}^{\infty}e^{-\mu
n-n/t}\mathrm{d}n=1/(t(1+\mu t))\approx 1/(\mu t^{2})$ (a similar calculation
applies to the extension).

At second order,

$F_{2}(x,t)=-\int_{t_{0}}^{t}\mathrm{d}s\int\mathrm{d}yG(x,t|y,s)F_{1}(y,s)^{2}\,.$
(56)

The computation of $\ell^{(2)}$ is a bit lengthy, see Appendix D.1. The result
is simple when $t\gg b^{2}$:

$\ell^{(2)}\simeq-\frac{1}{\mu^{2}t^{2}}\frac{4}{b}\,,\,N_{c}^{(2)}\simeq\frac{1}{\mu^{2}t^{2}}\frac{4}{b^{2}}\,,\,$
(57)

The higher order contributions $\ell^{(k)}$ and $N_{c}^{(k)}$ are more tedious
to compute, but their $t,\mu$ and $b$ dependence can be obtained by a
diagrammatic power-counting (see Appendix D.2):

$\ell^{(k)}\sim\frac{b^{3-2k}}{t^{2}\mu^{k}}\,,\,N_{c}^{(k)}\sim\frac{b^{2-2k}}{t^{2}\mu^{k}}\,.$
(58)

Since we are interested in the regime $b^{2}\gg\mu^{-1}$, these contributions
are more and more subleading as the order increases. Thus, it is justified to
stop the perturbation expansion at second order (by contrast, in the
$b^{2}\ll\mu$, the present perturbation theory would break down). Combining
both orders, we obtain

$\left<N_{c}(b)e^{-\mu n}\right>=N_{c}^{(1)}+N_{c}^{(2)}+\dots\sim\frac{1}{\mu
t^{2}}+\frac{4}{\mu^{2}t^{2}b^{2}}+\dots\,,\,1\ll\mu^{-1}\ll b^{2}\ll t\,.$
(59)

One might check, by (21) and (22), that (59) is equivalent to the following
gap distribution conditioned on the number of particles:

$\boxed{\left<N_{c}(b)\right>_{n}-1\simeq\frac{4n}{b^{2}}\,,\,1\ll\sqrt{n}\ll
b\ll\sqrt{t}\,.}$ (60)

Sending $t\to\infty$, we obtain the $b\gg\sqrt{n}$ case of the result (7)
announced in the introduction. As we discussed there, eq. (60) is in full
agreement with the universal gap distribution prediction (1) of [14, 15],
including the exact prefactor, as we discussed below (7). In fact, the
perturbative approach here is closely related to with the hierarchy of
traveling wave equations considered in these works, which allows to obtain
detailed information on the $k$-th gap. The order of expansion here
corresponds to the number of particles. It is thus not surprising that we need
to go to second order to access information about gaps.

##  IV Numerical study

We studied the model with two numerical methods: direct simulation, and
numerical integration of the KPP equations.

###  IV.1 Direct simulation

We measured directly the cluster number $\left<N_{c}(b)\right>$ in a discrete-
time analogue of the critical BBM. The time evolution takes place by
stroboscopic steps. In each step, each particle replicates itself (into 2
particles at the same position) with probability $1/2$, and is removed with
probability $1/2$. Then, each of the remaining particles displaces by an
independent Gaussian variable of zero mean and unity variance. This model is
expected to be in the same universality class as the continuous-time critical
BBM, while being simpler to simulate. We also applied an importance sampling
technique to generate more realizations with large number of particles, see
Appendix E.1 for details.

Figure 4: Results of the direct numerical simulation of a discrete-time
analogue of the critical BBM. Main: the average cluster number
$\left<N_{c}(b)\right>$ as a function of $b$, for various values of
$t=2^{8},2^{9},\dots,2^{17}$ (see inset for color code). The data are divided
by $\left<N_{c}(b_{0})\right>$, $b_{0}=5$ to highlight the $b$ dependence. The
dashed lines indicate the predicted power laws $b^{-2\eta-1}\approx b^{-1.78}$
(for small $b$) and $b^{2\eta-1}\approx b^{-0.22}$ (for large $b$), with
arbitrarily adjusted prefactors. Inset: The average cluster number with
$b=b_{0}$ restricted to non-extinct realizations. The solid line indicates the
predicted power law $t^{1/2-\eta}\approx t^{0.11}$, with arbitrarily adjusted
prefactors. See Appendix E.1 for numerical methods. We average over $10^{6}$
gaps for each value of $t$.

The main result of the simulation is presented in Figure 4. First, we plot the
average cluster number, conditioned on non-extinction, with a coarse-graining
scale of order unity ($b=b_{0}=5$), as a function of time. It shows a slow but
visible increase, in nice agreement with the prediction $t^{1/2-\eta}\approx
t^{0.11}$ (38) with $\eta\approx 0.39$ (29), see inset of Figure 4. Then, we
focus on the $b$ dependence. As $t$ increases, we observe the emergence of two
distinct power laws, with a crossover at $b\sim b_{0}$. The data are
compatible with the predictions $b^{-2\eta-1}\approx b^{-1.78}$ for small $b$
and $b^{2\eta-1}\approx b^{-0.22}$ for large $b$. We note however that the
relatively small value of the large-gap exponent makes it challenging to
confirm without ambiguity; indeed, in Figure 4, the large-gap power law
appears essentially flat compared to the small-gap one. We shall find more
convincing evidence by numerical integrating the KPP equation.

###  IV.2 Integrating the KPP equation

We have seen that precisely observing the large gap exponent in (38) requires
many orders of magnitude in space and time. In fact, it remains laborious even
if we resort to the semi-numerical approach of integrating the KPP equation.
To overcome this difficulty we used a integrate-and-coarse-grain scheme that
allows us to access large space-time scales with moderate computation
resource, see Appendix E.2 for details.

The main results are displayed in Figure 5. In the left panel, we plotted the
average cluster number, for large $b$, in the same way as in Figure 4 above,
i.e., factoring away the dependence on $t$. By going up to $t\sim 10^{8}$, we
observe a well-established power law $\left<N_{c}(b)\right>\sim b^{2\eta-1}$,
further corroborating the result (38).

We also integrated the KPP equation with $\mu\neq 0$, in order to observe the
crossover between the typical and rare gap regimes. More precisely, we test
the predictions (44) and (59), which can be reformulated as follows: in the
long time limit,

$-b\partial_{b}\left<N_{c}(b)e^{-\mu
n}\right>t^{2}\mu\sim\begin{dcases}(b\mu^{\frac{1}{2}})^{2\eta-1}&1\ll
b\ll\mu^{-\frac{1}{2}}\\\
8(b\mu^{\frac{1}{2}})^{-2}&b\gg\mu^{-\frac{1}{2}}\end{dcases}\,.$ (61)

where the prefactor in the $b\gg\mu^{-\frac{1}{2}}$ is exact, and the
$b\partial_{b}$ derivative is introduced to get rid of the $b$-independent
term in (59). This prediction is nicely confirmed by numerics, see Figure 5,
right panel. This result complements the numerical work in [14, 15] performed
on fixed particle number sectors.

Figure 5: Results obtained from numerical integration of the KPP equation.
Left panel: Average cluster number, with the predicted $t$ dependence removed,
$t^{\eta+1/2}\left<N_{c}(b)\right>$, for $b\in[5,5000]$ and
$t=20,\dots,8\times 10^{7}$ (represented by color, from red to blue). The
dashed line indicates the $b^{2\eta-1}$ ($\eta=0.39$) prediction, with
arbitrarily adjusted prefactor. Right panel: Testing the prediction (61) on
the cluster number with particle number conditioning, for a few values of
$\mu$. The value of $t$ alternates between $3.3\times 10^{8}$ and $1.3\times
10^{9}$. The blue dashed line plots the power laws $(b\mu^{1/2})^{2\eta-1}$
with an arbitrarily adjusted prefactor. The black solid lines depicts
$8(b\mu^{1/2})^{-2}$ with the predicted exact prefactor. See Appendix E.2 for
numerical methods.

##  V Discussion

We characterized the clustering structure of the 1D critical branching
Brownian motion, adopting an approach introduced in [24]. This allowed us
reconcile the result of [14, 15] focusing on gap statistics in fixed particle
number sectors, and that of [23] focusing on the continuum limit. Our approach
is based on a mapping to the KPP equation. Its asymptotic solution is always
given by a perturbative expansion, yet around different backgrounds depending
on the regime. Our analysis of the KPP equation is not rigorous, but is well
supported by numerics. We note that the analysis of the KPP equation —in
particle in the regime $t\ll\mu^{-1}$, corresponding to no conditioning on the
particle number — is reminiscent of that of the instanton equation in our
previous work [24]. There, we analyzed the clusters of all the positions
(integrated over time) visited by a critical branching fractional Brownian
motion. The number of clusters is also governed by a nontrivial exponent
governing the subleading asymptotics of the solution, analogous to the
exponent $\eta$ (and to the fractal dimension of the boundary
$D_{\mathrm{f}}=1-2\eta$) crucial to the present work.

Critical branching Brownian motion displays clustering (patchiness) in higher
dimensions as well. The results of this work apply to the projection of a
$d$-dimension BBM onto one direction. We can also consider generalizing the
present method to $d>1$, by assigning a patch (a disk in 2D, a ball in 3D, …)
of radius $b$ to each point, so that nearby balls overlap and form clusters. A
similar mapping to a KPP equation allows us to compute directly the avarage
$d$-dimensional volume $V$ occupied by the patches. However, differentiating
with respect to $b$ no longer yields the cluster number $N_{c}(b)$ for $d>1$
but other topological numbers. For example, in 2D, the second derivative of
the area $\partial_{b}^{2}\mathcal{A}$ is proportional to the difference
between the number of clusters and holes. Therefore, counting clusters in
higher dimensions remains an interesting problem. On the other hand, we expect
it to be rather straightforward to extend the present approach to 1D branching
fractional Brownian motion, i.e., random walks with long-range hopping;
results on this topic will be reported elsewhere.

###### Acknowledgements.

We thank Satya Majumdar for useful discussions, in particular for pointing out
the Reference [14]. P. Le Doussal acknowledges support from ANR grant
ANR-17-CE30-0027-01 RaMaTraF.

##  Appendix A Derivation of the KPP equation

We derive of the KPP equation (9) by a backward recursion argument. The
initial condition follows from definition. Now, the backward recursion
consists in considering what can happen during the initial elapse
$t\in[0,\mathrm{d}t]$:

  1. 1.

The initial particle branches into two with probability $\mathrm{d}t$. In that
case $E(x)\to E(x)^{2}$ (because from that moment, the two individuals act
independently from now on with the same law), thus giving a contribution
$(E^{2}-E)\mathrm{d}t$ to $\mathrm{d}E$.

  2. 2.

The initial particle is annihilated with probability $\mathrm{d}t$. Then
$E(x)\to 1$, giving a contribution $(1-E)\mathrm{d}t$ to $\mathrm{d}E$.

  3. 3.

The initial particle displaces to $\mathrm{d}x_{0}$, then $E(x)\to
E(x-\mathrm{d}x_{0})$, giving a contribution

$E(x-\mathrm{d}x_{0})-E(x)=-\partial_{x}E\left<\mathrm{d}x_{0}\right>+\frac{1}{2}\partial_{x}^{2}E\left<(\mathrm{d}x_{0})^{2}\right>=\partial_{x}^{2}E\mathrm{dt}$

to $\mathrm{d}E$: here we used Ito calculus (expanding to second order) and
that $\left<\mathrm{d}x_{0}\right>=0,\left<(\mathrm{d}x_{0})^{2}\right>=1$.

Combining all the contributions, we obtain (9).

##  Appendix B About the dense (super-Brownian motion) limit

The super-Brownian motion (SBM) can be defined as a dense limit of the
critical BBM. A common way to proceed (see e.g. [17, 18, 19]) is to consider a
large parameter $M\gg 1$, and perform the following steps:

  1. 1.

Set the branching/annihilation rate to be large $\beta=M$ (while keeping the
diffusion constant $D$ of order one).

  2. 2.

Modify the initial condition, and start with $n(0)\propto M$ particles at the
origin.

  3. 3.

Define the local density as

$\rho(x,t):=\frac{1}{M}\sum_{i=1}^{n(t)}\delta(x-x_{i}(t))\,,$ (62)

where $n(t)$ is the number of particles at time $t$ and
$x_{i}(t),i=1,\dots,n(t)$ their positions.

Then, it is known that the cumulant generation function of $\rho(x,t)$
satisfies a KPP equation with a delta-function initial condition. This is the
basic framework of [23]. The goal of this Appendix is to discuss that only
step 1 above is essential for the observables studied here.

To see why the initial condition is not essential, consider the critical BBM
with $\beta=M$ (step 1) but only a single particle initially (no step 2).
Then, consider what happens after a small lapse $M^{-1}\ll t\ll 1$. Since
$Mt\gg 1$, if the process has not gone extinct, it will have $\sim Mt$
particles. Since $t\ll 1$ and $D$ is of order one, these particles have no
time to diffuse and thus are still infinitesimally close to the origin.
Therefore, we have shown that, conditioned on non-extinction after an
infinitesimal time, starting with one particle is equivalent to starting with
$O(M)$ particles at the same position. We note that the procedure of
restricting to non-extinct realizations is often achieved in another way, by
rescaling the probability measure by $M$. That gives rise to a un-normalized
measure, called the “canonical measure” [17, 18]. In the context of mean-field
theory of avalanches (Brownian force model), the canonical measure corresponds
to studying the avalanche response with respect to an infinitesimal kick, see
[19], around eq. (25) therein.

Step 3 is a choice of observable. The ones considered in this work are
different; fundamentally, they are all related to the probability that
$[x-b/2,x+b/2]$ contains at least one particle at time $t$, denoted $F(x,t)$
[see (15) above with $\mu=0$]. If we consider a single initial particle at
$x=0$, $F(x,t)$ satisfies a KPP equation, which is the following with
dimensions restored:

$\partial_{t}F=D\partial_{x}^{2}F-\beta F^{2}\,,\,F(x,t=0)=\theta(b/2-|x|)\,.$
(63)

Now, in the SBM limit $\beta=M\gg 1$, we can consider $F=\tilde{F}/M$ so that

$\partial_{t}\tilde{F}=D\partial_{x}^{2}\tilde{F}-\tilde{F}^{2}\,,\,\tilde{F}(x,t=0)=M\theta(b/2-|x|)\,.$
(64)

Initially, the nonlinear term dominates with respect to the diffusion one, so
that

$\tilde{F}\approx\frac{1}{M^{-1}+t}\theta(b/2-|x|)$ (65)

This approximation is valid until $\tilde{F}$ becomes of order unity. Observe
that the RHS of the above equation has a limit as $M\to\infty$ (uniformly in
$t\in[t_{0},\infty)$ for any fixed $t_{0}>0$). Therefore, we may conclude that
in the SBM limit, the solution to (63) has the following form:

$F(x,t)=\frac{1}{M}\tilde{F}(x,t)$ (66)

where $\tilde{F}(x,t)$ satisfies a KPP equation with a diverging initial
condition

$\partial_{t}\tilde{F}=D\partial_{x}^{2}\tilde{F}-\tilde{F}^{2}\,,\,\tilde{F}(t\to
0)\simeq\frac{1}{t}\theta(b/2-|x|)\,.$ (67)

To interpret these formulas, recall that $1/M$ is the probability that the
process is not extinct at $t_{0}=1$. Hence, by (66), for $t>t_{0}$,
$\tilde{F}(x,t)$ is the probability that $[x-b/2,x+b/2]$ is not empty at $t$,
conditioned on non-extinction at an early time. (As per the discussion above,
$\tilde{F}(x,t)$ is also the probability of $[x-b/2,x+b/2]$ not being empty at
$t$ with an initial condition of $O(M)$ particles near the origin.) Therefore,
the observables we study in this work have a well-defined behavior in the SBM
limit. This limit can be directly studied via the KPP equation (67). In this
work, we took the alternative approach: we focus on the BBM, and then discuss
the SBM limit of the results, see (6) in Section II.2.

##  Appendix C Crossover between the two asymptotic Ansätze

Figure 6: Crossover between the two asymptotic Ansätze. The left panel plots
the solution to (69) that interpolates between the $t\ll\mu^{-1}$ Ansatz (40)
and the $t\gg\mu^{-1}$ one (41), as a function of $y=x/\sqrt{t}$, for
$s=\ln(t\mu^{-1})\in[-7,7]$ (from bottom to top, also indicated by color, see
right panel for color code). The right panel plots the $s$ dependence of the
center, $f(y=0,s)$, from which we subtracted
$f(y=\infty,s)=t/(\mu^{-1}+t)=e^{s}/(1+e^{s}).$ It decays as $e^{-3s/2}$, with
a prefactor $A_{1}\approx 1.02$ (best fit, red dashed line). This provides an
estimate of the amplitude $A_{1}$ in (74), which is equal to $c_{0}^{\prime}$
in (43).

In this appendix we characterize in more detail the crossover between the
$t\ll\mu^{-1}$ and $t\gg\mu^{-1}$ asymptotic Ansätze discussed in Section
III.4.1.

###  C.1 Leading term

We first focus on the crossover between leading-order terms. For this, we
perform a change of variables,

$F(x,t)=\frac{1}{t}f(y,s)\,,\,\text{where }y=x/\sqrt{t},s=\ln(t\mu)\,.$ (68)

Here and below, to avoid confusion, we will write $f(y,s)$ to refer to the
time-dependent function just defined, while $f(y)$ stands for the leading term
of the Ansatz (25). Then the KPP equation
$\partial_{t}F=\partial_{x}^{2}F-F^{2}$ is equivalent to

$\partial_{s}f(y,s)=\partial_{y}^{2}f(y,s)+\frac{y}{2}\partial_{y}f(y,s)+f(y,s)-f(y,s)^{2}\,.$
(69)

This is the time-dependent generalization of (26) above. The leading terms of
the two Ansätze — $f(y)$ in (40) (plotted in Figure 3), and $f_{1}(y)=1$ in
(41) — are time-independent solutions of (69). The crossover between them is
described by a time-dependent solution, which tends to $f(y)$ as $s\to-\infty$
and to $f_{1}(y)=1$ as $s\to\infty$ (note that the rescaled time $s=\ln(\mu
t)$ is the logarithm of the physical time, so $s\to-\infty$ and $s\to\infty$
correspond to $t\ll\mu^{-1}$ and $t\gg\mu^{-1}$, respectively). Also, its
large-distance behavior is also fixed as
$f(|y|\to\infty)={t}F_{0}={t}/(\mu^{-1}+t)={e^{s}}/(1+e^{s})$, see (17) above.
Numerical integration of the KPP equation (the original one, not the rescaled
one studied here; see Appendix A for methods) indicates that the crossover
solution with the above boundary conditions is unique. Indeed, we observe that
solutions to the KPP equation with different values of $\mu$ coincide with
each other in the re-scaled coordinate system $(y,s,f)$, as long as $1\lesssim
b^{2}\ll\mu^{-1}$. In Figure 6, we plot the crossover solution obtained
numerically.

Let us characterize the limiting behaviors of the crossover solution:

  * •

At large distances

$f(|y|\to\infty,s)\simeq\frac{e^{s}}{1+e^{s}}\,,$ (70)

which increases from $0$ to $1$ as $s$ increases from $-\infty$ to $+\infty$.

  * •

As $s\to-\infty$, $f(y,s)$ tends to $f(y)$ of (40) and (25), plotted in Figure
3. Linearizing (69) around that fixed point, we find that the difference
$\delta f(y,s)=f(y,s)-f(y)$ satisfies the time-dependent analogue of (27)
above,

$\partial_{s}(\delta f)=-H(\delta
f)\,,\,H=-\partial_{y}^{2}-\frac{y}{2}\partial_{y}+2f(y)-1\,.$ (71)

We checked that $H$ has an eigenfunction $\psi(y)$ with a negative eigenvalue
$-1$, corresponding to an instability. This should be distinguished from the
eigenvalue $\eta>0$ which corresponds to a decaying eigenfunction, see Section
III.3.1 above. A plot of $\psi(y)$ can be found in Figure 7. In particular,
$\psi(|y|\to\infty)$ tends to a nonzero constant; hence, it corresponds to a
non-normalizable wavefunction of the Schrödinger Hamiltonian (28). The
$s\to-\infty$ asymptotics of $\delta f$ is dominated by this unstable mode:

$f(y,s\to-\infty)=f(y)+e^{s}\psi(y)+\dots\,,\,\psi(|y|\to\infty)\simeq 1\,.$
(72)

Note that the $|y|\to\infty$ limit is in agreement with (70).

  * •

As $s\to+\infty$, $f(y,s)$ tends to $f_{1}(y)=1$. Linearization around this
fixed point leads to (71) with the Hamiltonian replaced by

$H=-\partial_{y}^{2}-\frac{y}{2}\partial_{y}+1\,.$ (73)

Its leading eigenfunctions are the constant $\varphi_{0}(y)=1$ with eigenvalue
$1$, and the Gaussian $\varphi_{1}(y)=e^{-y^{2}/4}$ with eigenvalue $3/2$.
Therefore, we have

$f(y,s\to+\infty)=1+A_{0}e^{-s}\varphi_{0}(y)+A_{1}e^{-3s/2}\varphi_{1}(y)+\dots=1-e^{-s}+A_{1}e^{-3s/2}e^{-y^{2}/4}+\dots\,.$
(74)

Here, the amplitude $A_{0}=-1$ has been fixed by comparing with (70). We
estimated numerically that $A_{1}\approx 1.02$, see 6 (right panel). Note that
$A_{1}$ is equal to the coefficient $c_{0}^{\prime}$ in (43).

Figure 7: The unstable eigenfunction $\psi(y)$ of the Hamiltonian $H$ (71)
with eigenvalue $-1$. It tends to $1$ as $|y|\to\infty$.

###  C.2 Subleading term

Figure 8: Crossover between the subleading term of the two asymptotic Ansätze.
The left panel plots the solution to (75) with the asymptotic condition (76),
as a function of $y$, for $s=\ln(t\mu^{-1})\in[-10,10]$ (indicated by color,
see right panel for color code). The right panel plots the $s$ dependence of
the center, $g(y=0,s)$. It decays as $e^{-3s/2}$ for $s\to\infty$, with a
prefactor $A_{g}\approx 2.2$ (best bit, red dashed line).

We now consider a small perturbation around the leading-order crossover
solution studied above, that is, we replace $f(y,s)\to f(y,s)+\varepsilon
g(y,s)$ in (68) and (69), where $\varepsilon$ is a small parameter (to be
fixed below). This leads to a linearized equation for $g$:

$\partial_{s}g(y,s)=-H(s)g(y,s)\,,\,H(s):=-\partial_{y}^{2}-\frac{y}{2}\partial_{y}+2f(y,s)-1\,.$
(75)

It has a time-dependent “Hamiltonian” that interpolates between (71) and (73).
A unique solution is specified by the asymptotic condition

$g(y,s)\simeq e^{-\eta s}\varphi(y)\,,s\to-\infty\,.$ (76)

where $\eta\approx 0.39$ and $\varphi$ are the leading decaying eigenvalue and
eigenfunction of (71) [see below (27) for further discussion, and Figure 3 for
a plot.] The full solution can be again extracted from the numerical
integration of the KPP equation, and is plotted in Figure 8. As $s\to\infty$,
it has the following asmyptotics:

$g(y,s)\simeq A_{g}e^{-3s/2}e^{-y^{2}/4}\,,s\to+\infty\,,\,$ (77)

where we estimated $A_{g}\approx 2.2$.

We now fix the small parameter $\varepsilon$. For this, we match the
$t\ll\mu^{-1}$ (or $s\to-\infty$) behavior of the crossover solution with
perturbation

$F(x,t)=\frac{1}{t}\left[f(x/\sqrt{t},\ln(\mu t))+\varepsilon g(x,\ln(\mu
t))\right]=\frac{1}{t}\left[f(x/\sqrt{t})+\varepsilon(\mu
t)^{-\eta}\varphi(x/\sqrt{t})+\dots\right]$ (78)

with the Ansatz (40), which is

$F(x,t)=\frac{1}{t}\left[f(x/\sqrt{t})+Ct^{-\eta}\varphi(x/\sqrt{t})\right]\,,\,C\simeq
A_{c}b^{2\eta}\,.$

($A_{c}$ is a constant independent of $b,t,\mu$). As a result, we find
$\varepsilon=C\mu^{\eta}\sim b^{2\eta}\mu^{\eta}$. This, combined with (74)
and (77), implies the following $t\gg\mu^{-1}$ (or $s\to+\infty$) asymptotics:

$\displaystyle F(x,t)=$ $\displaystyle\frac{1}{t}\left[f(x/\sqrt{t},\ln(\mu
t))+\varepsilon g(x,\ln(\mu t))\right]$ $\displaystyle=$
$\displaystyle\frac{1}{t}\left[1-(\mu t)^{-1}+\left(A_{1}(\mu
t)^{-3/2}+A_{g}A_{c}t^{-3/2}b^{2\eta}\mu^{\eta-\frac{3}{2}}\right)e^{-\frac{x^{2}}{4t}}+\dots\right]$
(79)

This is exactly (43); in particular, the coefficients in (43)
$c_{0}^{\prime}=A_{1}$ and $c_{1}=A_{g}A_{c}$ are indeed order one constants
that do not depend on $b,t$ and $\mu$.

In conclusion, a more elaborate analysis of the crossover between the two
Ansätze confirms the simple matching argument presented in Section III.4.1.

##  Appendix D Perturbation calculations

###  D.1 Second order

In this appendix, we derive (57). We will work in the regime $1\ll\mu^{-1}\ll
t_{0}\ll b^{2}\ll t$ and will always use them to simplify expressions. First,
by (56) and (50),

$\displaystyle\ell^{(2)}$
$\displaystyle\sim-a^{2}\int_{t_{0}}^{t}\mathrm{d}s\int\mathrm{d}x\mathrm{d}yG(x,t|y,s)F_{1}(y,s)^{2}$
(80) $\displaystyle\sim$
$\displaystyle-a^{2}\int_{t_{0}}^{t}\mathrm{d}s\frac{s^{2}}{t^{2}}\int\mathrm{d}x\mathrm{d}y\,p(x-y,t-s)F_{1}(y,s)^{2}$
(81) $\displaystyle=$
$\displaystyle-a^{2}\int_{t_{0}}^{t}\mathrm{d}s\frac{s^{2}}{t^{2}}\int\mathrm{d}yF_{1}(y,s)^{2}$
(82)

where in the last line, we used the property $\int p(z,u)\mathrm{d}z=1$ of the
standard diffusion propagator. Now, plugging in (54),

$\displaystyle\ell^{(2)}$
$\displaystyle\sim-\frac{1}{\mu^{2}t^{2}}\int_{t_{0}}^{t}\frac{\mathrm{d}s}{s^{2}}\int\mathrm{d}y\int_{-\frac{b}{2}}^{\frac{b}{2}}\mathrm{d}u\int_{-\frac{b}{2}}^{\frac{b}{2}}\mathrm{d}vp(y-u,s-t_{0})p(y-v,s-t_{0})$
(83)
$\displaystyle=-\frac{1}{\mu^{2}t^{2}}\int_{t_{0}}^{t}\frac{\mathrm{d}s}{s^{2}}\left[b\,\mathrm{erf}\left(\frac{b}{2\sqrt{2(s-t_{0})}}\right)-\sqrt{\frac{8(s-t_{0})}{\pi}}\left(1-e^{-\frac{b^{2}}{8(s-t_{0})}}\right)\right]$
(84) $\displaystyle\sim-\frac{1}{\mu^{2}t^{2}t_{0}}\left[\sqrt{2\pi
t_{0}}e^{\frac{b^{2}}{8t_{0}}}\left(1-\text{erf}\left(\frac{b}{2\sqrt{2}\sqrt{t_{0}}}\right)\right)+b-\sqrt{2\pi
t_{0}}\right]$ (85)

The integral of the second line is exact (recall
$\mathrm{erf}(x)=\frac{2}{\sqrt{\pi}}\int^{x}_{0}e^{-z^{2}}\mathrm{d}z$), and
the last line is exact for $t\gg b^{2}$. We can further simplify using
$t_{0}\ll b^{2}$:

$\ell^{(2)}\sim-\frac{1}{\mu^{2}t^{2}}\frac{4}{b}-\frac{1}{\mu^{2}t^{2}}\frac{b}{t_{0}}\,.$
(86)

The second term here is similar to $\ell^{(1)}$ (55), but is $a=1/(\mu t_{0})$
times smaller. Dismissing it gives us the desired result (57).

###  D.2 Higher orders: diagrammatic expansion and power-counting

The perturbation expansion of $\left<\ell(b)\right>$ in Section III.4.2 can be
represented by diagrams, see Fig. 9. This diagrammatic expansion is closely
related to that of the Brownian force model of mean-field avalanches [20, 21].
The diagrams at order $k$ are trees with one root, $k-1$ internal vertices,
and $k$ leaves. The Feynmann rules for calculating the contribution of a
diagram to $\ell^{(k)}$ are the following:

Figure 9: Diagrammatic representations of the perturbation expansion at order
$1$ (left), $2$ (center) and $3$ (right, plus its mirror image).

  1. 1.

To each vertex (i.e., internal vertices, root and leaves) is associated a
space-time point.

  2. 2.

To each edge is associated with one retarded propagator (53)

$G(x,t|y,s)=\frac{(\mu^{-1}+s)^{2}}{(\mu^{-1}+t)^{2}}p(x-y,t-s)\theta(t-s)\,,\,p(z,u)=\frac{e^{\frac{z^{2}}{4\pi
u}}}{\sqrt{4\pi u}}$

between the vertices it connects. The later time $t$ corresponds to the vertex
closer to the root (upper vertex in Figure 9). We will assign the
$(\mu^{-1}+s)^{2}$ to the vertex at time $s$ and
$1/(\mu^{-1}+t)^{2}/\sqrt{t-s}$ to that at time $t$.

  3. 3.

The root is integrated only over space $\in\mathbb{R}$, with time fixed to
$t$.

  4. 4.

The leaves are integrated only over space $x\in[-b/2,b/2]$, with time fixed to
$t_{0}=1/(a\mu)\gg\mu$. Moreover, to each leaf is associated with a factor
$1/(\mu^{-1}+t_{0})$, see (48).

  5. 5.

The internal vertices are integrated over both space $\in\mathbb{R}$ and time
$\in\mathbb{R}$ (then constrained by the causality of the propagator).

  6. 6.

Finally, we multiply the result by $a^{k}$, as per (50).

Now let us count the powers of $t,\mu$ and $b$ of a diagram. First, the root’s
integral can be performed independently of the rest of the diagram, giving a
contribution $1/(\mu^{-1}+t^{2})\sim 1/t^{2}$ (since $t\gg 1/\mu^{-1}$) from
the numerator of the propagator. $t$ will not occur in the rest of the
diagram, except as an upper time limit that will be sent to infinity (since we
are interested in the $t\to\infty$ asymptotics). So the power of $t$ is simply
$1/t^{2}$ for any diagram.

Next, $\mu$ can only occur in the numerator $(\mu^{-1}+t_{0})^{2}$ of the
Green function connecting a leaf (rule 2), and the factor $1/(\mu^{-1}+t_{0})$
associated to it (rule 4). Since there are $k$ leaves, we have a factor
$(\mu^{-1}+t_{0})^{k}\sim t_{0}^{k}\sim a^{-k}\mu^{-k}$ since
$t_{0}\gg\mu^{-1}$. Combined with the $a^{k}$ factor (rule 6), we have a
$\mu^{-k}$ power: this is the $\mu$ dependence. By the way, we see that the
$a$-dependence has been cancelled: the result of the calculation has a well-
defined limit as $a\to 0$.

Finally, we count the powers of $b$. For this, we count the length $[x]$ and
time $[t]$ powers that remain, i.e., those associate with the space-integral
of the leaves, and those associated with the internal vertices (i.e.,
excluding leaves and the root). Each leaf-integral gives a $[x]$ power, which
adds to $[x^{k}]$ (there being $k$ leaves), by rule 4. For each internal
vertex, it has a $[xt]$ power because of its space-time integral, by rule 5.
The three propagator above it (see Figure 9) contributes $[t^{2}]$, and the
two below it contribute $([t^{-1/2}][t^{-2}])^{2}$, by rule 2. Converting
$[t]=[x^{2}]$ (by diffusive scaling of the propagator), we obtain
$[xt][t^{2}]([t^{-1/2}][t^{-2}])^{2}=[xt^{-2}]=[x^{-3}]$ per internal vertex.
Since there are $k-1$ vertices, we have $[x^{-3k+3}]$ for all internal
vertices. Adding both leaf and internal vertex contributions, we get
$[x^{-2k+3}]$, which translates to a $b^{-2k+3}$ dependence on $b$.

In summary, we have shown that every diagram of $\ell^{(k)}$ contributes a
power

$t^{-2}\mu^{-k}b^{-2k+3}\,,$ (87)

as announced in (58) above.

##  Appendix E Numerical methods

###  E.1 Importance sampling

The idea is to sample the realization with a statistical weight biased by the
factor $e^{\mu n(T)}$, where $n(T)$ is the number of particle at the maximal
time step $T$ of the simulation, and $\mu>0$ is a chemical potential, chosen
positive to favor realizations with large number of particles. It should be
neither too small (no effect) nor too large (causes divergence); in practice
$\mu=1/t$ is a good compromise. This bias is then corrected by an inverse
weight $e^{-\mu n(t_{\max})}$ when computing observables.

Thanks to the tree structure of the BBM model, the weighed ensemble can be
sampled directly. To see how, consider the partition function

$Z(\mu,T):=\left<e^{\mu n(T)}\right>\,.$ (88)

By definition, the importance sampling consists in generating configurations
of the statistical ensemble defined by $Z$. By a backward recursion argument,
it is not hard to show

$Z(\mu,T)=\frac{1}{2}(1+Z(\mu,T-1)^{2})\,,\,T>0\,,\,$ (89)

and $Z(\mu,0)=e^{\mu}$. Using this recursion relation we can easily compute
numerically $Z(\mu,t)$ for $t=0,1,\dots,T$. Then, a moment of thought shows
that, to sample $Z(\mu,T)$, we should simulate the branching process with a
sequence of modified weights: at the $k$-th step, each particle branches with
probability $1/(2Z(\mu,T-k+1))$ and annihilates with probability
$Z(\mu,T-k)/(2Z(\mu,T-k+1))$ (we use the sequence of partition functions
backwards in time, since we work with a backward recursion). The random-walk
part of the simulation is not affected by the importance sampling.

###  E.2 Integration and coarse-grain scheme

The numerical integration of the KPP equation involves an integrate-and-
coarse-grain scheme, reminiscent of a real-space renormalization group, in
order to reach large time and space with moderate resource. Indeed, we expect
the solution to becomes smoother and smoother in time. Therefore, we can
gradually coarse-grain the solution, discarding short-distance information,
without losing precision. This is done by iteratively applying the rescaling
$t=\alpha^{2}\tilde{t}$, $x=\alpha\tilde{x},F=\alpha^{-2}\tilde{F}$ (31),
under which the KPP equation is invariant.

More concretely, we approximate the space by a lattice with finite size
$[-A/2,A/2]$ and spacing $\epsilon$. $F$ is represented by an array
$N=A/\epsilon$ data points. We execute the following steps:

  1. 1.

Initial step. Integrate the KPP equation from the initial condition for a
duration of $t_{0}$.

  2. 2.

Coarse-grain. We remove every other data points of the solution just obtained
(so $\alpha=2$), multiply it by $4$, and pad both ends of the array by the
edge value repeated $N/4$ times:

$(F_{1},\dots,F_{N})\leftarrow(\underbrace{F_{1},\dots,F_{1}}_{\text{{\hbox{N/4}}
times}},F_{1},F_{3},F_{5}\dots,F_{N-1},\underbrace{F_{N-1},\dots,F_{N-1}}_{\text{{\hbox{N/4}}
times}})$ (90)

As a result, the lattice remains the same but it now represents the rescaled
$\tilde{x}$ variable. Note that the time is now rescaled to $t_{0}/4$.

  3. 3.

Integrate. Evolve the solution by the KPP equation for a duration of
$t_{0}-t_{0}/4$ so that the (rescaled) time is again $t_{0}$.

We can then carry on by looping steps 2 and 3 indefinitely. After $k$
iterations, the algorithm outputs the rescaled solutin
$\tilde{F}(\tilde{x},\tilde{t})$ at $\tilde{t}=t_{0}$. It is related to the
physical solution in original coordinates $F=\tilde{F}4^{-k}$,
$x=\tilde{x}2^{k}$, $t=t_{0}4^{k}$. Hence, we obtain the solution at a space-
time scale that is an exponential of the computation time; the memory cost is
constant. The geometric quantities are calculated using (18) and (20). To
handle the $\partial_{b}$ derivative involved in the cluster number (20), we
integrate an auxiliary system for $G=\partial_{b}F$ along with $F$ using the
same techniques. Integrating $G$ gives the cluster number directly.

For the scheme to be accurate, the time $t_{0}$ should be large enough: the
solution at $t=t_{0}$ should be sufficient smooth so that the coarse-grain
step does not erase too much information. The lattice is then chosen to afford
a good approximation of the solution for $t\in[0,t_{0}]$. In practice, a
moderate set of parameters $t_{0}=20$, $\epsilon=0.25$ and $A=10^{3}$ is
adequate for our purposes. The integration of the KPP equation is performed by
a first-order Trotterization scheme (with $\delta t=0.1$), in which the
diffusion and nonlinear terms are treated in momentum and real space,
respectively.

## References

  * Zhang  _et al._ [1990] Y. C. Zhang, M. Serva, and M. Polikarpov, Diffusion reproduction processes, Journal of Statistical Physics 58, 849 (1990).
  * Tsimring  _et al._ [1996] L. S. Tsimring, H. Levine, and D. A. Kessler, RNA Virus Evolution via a Fitness-Space Model, Phys. Rev. Lett. 76, 4440 (1996).
  * Houchmandzadeh [2002] B. Houchmandzadeh, Clustering of diffusing organisms, Phys. Rev. E 66, 052902 (2002).
  * Houchmandzadeh [2008] B. Houchmandzadeh, Neutral clustering in a simple experimental ecological community, Phys. Rev. Lett. 101, 078103 (2008).
  * Houchmandzadeh [2009] B. Houchmandzadeh, Theory of neutral clustering for growing populations, Phys. Rev. E 80, 051920 (2009).
  * Bailey  _et al._ [1975] N. T. Bailey _et al._ ,  _The mathematical theory of infectious diseases and its applications_ (Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE., 1975).
  * Kendall [1956] D. G. Kendall, Deterministic and stochastic epidemics in closed populations, in  _Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Volume 4: Contributions to Biology and Problems of Health_ (University of California Press, 1956) pp. 149–165.
  * Dumonteil  _et al._ [2013] E. Dumonteil, S. N. Majumdar, A. Rosso, and A. Zoia, Spatial extent of an outbreak in animal epidemics, Proceedings of the National Academy of Sciences 110, 4239 (2013).
  * Meyer  _et al._ [1996] M. Meyer, S. Havlin, and A. Bunde, Clustering of independently diffusing individuals by birth and death processes, Phys. Rev. E 54, 5567 (1996).
  * Lawson and Jensen [2007] D. J. Lawson and H. J. Jensen, Neutral evolution in a biological population as diffusion in phenotype space: Reproduction with local mutation but without selection, Phys. Rev. Lett. 98, 098102 (2007).
  * Pázsit and Pál [2007] I. Pázsit and L. Pál,  _Neutron fluctuations: A treatise on the physics of branching processes_ (Elsevier, 2007).
  * Zoia  _et al._ [2014] A. Zoia, E. Dumonteil, A. Mazzolo, C. de Mulatier, and A. Rosso, Clustering of branching brownian motions in confined geometries, Phys. Rev. E 90, 042118 (2014).
  * de Mulatier  _et al._ [2015] C. de Mulatier, E. Dumonteil, A. Rosso, and A. Zoia, The critical catastrophe revisited, Journal of Statistical Mechanics: Theory and Experiment 2015, P08021 (2015).
  * Ramola  _et al._ [2014] K. Ramola, S. N. Majumdar, and G. Schehr, Universal order and gap statistics of critical branching brownian motion, Phys. Rev. Lett. 112, 210602 (2014).
  * Ramola  _et al._ [2015] K. Ramola, S. N. Majumdar, and G. Schehr, Branching brownian motion conditioned on particle numbers, Chaos, Solitons & Fractals 74, 79 (2015).
  * Kolmogorov  _et al._ [1937] A. Kolmogorov, I. Petrovsky, and N. Piscounov, Etude de l’équation de la diffusion avec croissance de la quantit é de matière et son application à un problème biologique, Bull. Univ. Etat Moscou A 1, 1 (1937).
  * [17] Ed. Perkins, Super-Brownian motion and critical spatial stochastic systems, Can. Math. Bull., 47(2), 280 (2004), see also http://www.math.ubc.ca/~perkins/superbrownianmotionandcriticalspatialsystems.pdf.
  * Slade [2002] G. Slade, Scaling limits and super-brownian motion, Notices AMS 49, 1056 (2002).
  * Le Doussal [2022] P. Le Doussal, Equivalence of mean-field avalanches and branching diffusions: From the brownian force model to the super-brownian motion, Journal of Physics A: Mathematical and Theoretical (2022).
  * Le Doussal and Wiese [2012] P. Le Doussal and K. J. Wiese, Distribution of velocities in an avalanche, EPL (Europhysics Letters) 97, 46004 (2012).
  * Le Doussal and Wiese [2013] P. Le Doussal and K. J. Wiese, Avalanche dynamics of elastic interfaces, Phys. Rev. E 88, 022106 (2013).
  * pie [2022] P. Le Doussal, More on the Brownian force model: avalanche shapes, tip driven, higher d, ArXiv:2203.10544 (2022).
  * Mueller  _et al._ [2017] C. Mueller, L. Mytnik, and E. Perkins, On the boundary of the support of super-Brownian motion, The Annals of Probability 45, 3481 (2017).
  * Cao  _et al._ [2022] X. Cao, P. Le Doussal, and A. Rosso, Clusters in an epidemic model with long-range dispersal, Phys. Rev. Lett. 129, 108301 (2022).
  * H.Brezis  _et al._ [2016] H.Brezis, L.A.Peletier, and D.Terman, A very singular solution of the heat equation with absorption, Arch. Rational Mech. Anal 95, 185 (2016).