# Kinetics of information scrambling in correlated metals:
disorder-driven transition from shock-wave to FKPP dynamics

Camille Aron<EMAIL_ADDRESS>Laboratoire de Physique de l’École Normale
Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris
Cité, F-75005 Paris, France  Institute of Physics, École Polytechnique
Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Éric Brunet
Laboratoire de Physique de l’École Normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France  Aditi
Mitra  Center for Quantum Phenomena, Department of Physics, New York
University, 726 Broadway, New York, NY, 10003, USA

###### Abstract

Quenched disorder slows down the scrambling of quantum information. Using a
bottom-up approach, we formulate a kinetic theory of scrambling in a
correlated metal near a superconducting transition, following the scrambling
dynamics as the impurity scattering rate is increased. Within this framework,
we rigorously show that the butterfly velocity $v$ is bounded by the light
cone velocity $v_{\rm lc}$ set by the Fermi velocity. We analytically identify
a disorder-driven dynamical transition occurring at small but finite disorder
strength between a spreading of information characterized at late times by a
discontinuous shock wave propagating at the maximum velocity $v_{\rm lc}$, and
a smooth traveling wave belonging to the Fisher or Kolmogorov-Petrovsky-
Piskunov (FKPP) class and propagating at a slower, if not considerably slower,
velocity $v$. In the diffusive regime, we establish the relation
$v^{2}/\lambda_{\rm FKPP}\sim D_{\rm el}$ where $\lambda_{\rm FKPP}$ is the
Lyapunov exponent set by the inelastic scattering rate and $D_{\rm el}$ is the
elastic diffusion constant.

Information scrambling refers to the efficient spreading and loss of
information throughout an extended many-body system. Its characterization is
fundamental to the foundations of quantum chaos and has implications in the
development quantum computing technologies. The picture that has emerged from
the exact or numerical solutions to a variety of classical and quantum
thermalizing models with local interactions is a ballistic spreading of
information at a so-called butterfly velocity. In dual-unitary circuits, the
information scrambling occurs precisely on the light rays propagating at the
maximum velocity allowed by causality [1, 2]. In random quantum circuits made
of Haar-distributed unitaries, scrambling dynamics was related to classical
growth processes with slower ballistic fronts that broaden either diffusively
in $d=1$, or according to fluctuations governed by the Kardar-Parisi-Zhang
universality class in $d=2$ and $d=3$ [3, 4].

Alongside these minimal models, as well as classical [5, 6, 7, 8], semi-
classical, large $N$ or holographic models [9, 10, 11, 12], it is essential to
address those questions in realistic situations where exact or numerical
solutions are out of reach. Aleiner, Faoro and Ioffe articulated those
questions in the larger framework of electronic transport by deriving a
quantum kinetic equation for out-of-time-ordered correlators (OTOCs) within a
so-called many-world Keldysh formalism [13]. They argued that the scrambling
dynamics in metals with either phonon or Coulomb interactions are governed by
Fisher or Kolmogorov-Petrovsky-Piskunov (FKPP) equations, resulting in smooth
non-broadening scrambling fronts propagating at a butterfly velocity set by
the Fermi velocity, $v_{\rm F}$. Interestingly, it was argued that the
presence of disorder would significantly reduce this butterfly velocity [13,
14, 15, 16].

Building on the approach of Ref. [13], we worked out the kinetics of quantum
information scrambling in a paradigmatic model of clean interacting metals in
the vicinity of a superconducting phase transition, where electron-electron
interactions are dominated by superconducting fluctuations [17]. We found
scrambling fronts that travel at $v_{\rm F}$ but do not belong to the FKPP
universality class. Remarkably, their late-time spatial profiles develop a
shock-wave discontinuity at the boundary of the light cone.

In this Letter, we investigate the impact of impurity scattering on the
dynamics of the scrambling front. Momentum relaxation is indeed significant in
metals with an elastic timescale typically in the tens of femtoseconds, i.e.
much shorter than the scattering time due to electronic interactions,
especially at low temperatures. We first derive an effective set of two
coupled partial differential equations (PDEs) governing the scrambling
dynamics at large scales. Then, we analytically elucidate how the traveling-
wave solutions that develop at late times are affected by the presence of
disorder, from the clean case to the diffusive regime. We find that the shock-
wave phenomenology found in the clean case is robust against weak disorder.
However, we unravel a dynamical phase transition that causes the scrambling
kinetics to abruptly conform to the FKPP class when the disorder strength
exceeds a critical value, see Fig. 1. We fully characterize the information
scrambling in the FKPP regime by working out the profiles of the traveling
fronts and their butterfly velocities.

##### Model.

$0$$0.5$$1$$0$${\kappa^{*}}$$10$$20$

$v\propto 1/\sqrt{\kappa}$

SHOCK WAVE

FKPP

$v$

$\kappa$$v_{\rm lc}$ Figure 1: Quantum butterfly velocity $v$ as a function of
the disorder strength $\kappa$. For $\kappa<\kappa^{*}$, the front propagates
at the light cone velocity $v_{\rm lc}=1$ in units of $v_{\rm F}/\sqrt{d}$.
For $\kappa>\kappa^{*}$, the solid line is the FKPP prediction made in Eq.
(22). The red marks are numerical results obtained by solving Eqs. (9) for
$d=1$ up to times $\tau=3000$ ($\gamma=1$ i.e. $\kappa^{*}:=1+\gamma=2$).

We consider a system of interacting electrons in $d\geq 2$ dimensions that are
subject to both elastic and inelastic scattering. The former is due to static
non-magnetic impurities and defects. The latter is due to electron-electron
interactions in the Cooper channel. This choice is guided by the relatively
straightforward treatment of superconducting fluctuations within the random
phase approximation. However, our approach can be extended to other models
with different interactions as long as a quasi-particle description is valid.
For concreteness, we have in mind the Hamiltonian

$\displaystyle H=$
$\displaystyle\sum_{{\bm{k}},\,\sigma}\epsilon_{\bm{k}}c^{\dagger}_{{\bm{k}}\sigma}c_{{\bm{k}}\sigma}+\sum_{{\bm{k}}{\bm{k}}^{\prime},\,\sigma}V_{{\bm{k}}-{\bm{k}}^{\prime}}c^{\dagger}_{{\bm{k}}\sigma}c_{{\bm{k}}^{\prime}\sigma}$
$\displaystyle+U\\!\\!\\!\sum_{{\bm{k}}{\bm{k}}^{\prime}{\bm{k}}^{\prime\prime},\,\sigma}\\!\\!c^{\dagger}_{{\bm{k}}\sigma}c^{\dagger}_{-{\bm{k}}+{\bm{k}}^{\prime}\bar{\sigma}}c_{{\bm{k}}^{\prime\prime}+{\bm{k}}^{\prime}\bar{\sigma}}c_{-{\bm{k}}^{\prime\prime}\sigma}.$
(1)

The operator $c_{{\bm{k}}\sigma}^{\dagger}$ creates a fermion with spin
$\sigma=\uparrow$ or $\downarrow$ ($\bar{\sigma}=\downarrow$, $\uparrow$) and
momentum ${\bm{k}}$ in the Brillouin zone. $\epsilon_{\bm{k}}$ is the
dispersion relation. Electronic energies are measured relative to the chemical
potential and $E_{\rm F}$ is the Fermi energy. For simplicity, we shall assume
a spherical Fermi surface, i.e. $\epsilon_{\bm{k}}=0$ when $k\to k_{\rm F}$.

The attractive interaction $U<0$ facilitates superconductivity. In dimensions
$d>2$, this model exhibits a finite-temperature phase transition towards a
superconducting phase associated with the spontaneous breaking of the $U(1)$
symmetry. In $d=2$, it is replaced by a Berezinskii–Kosterlitz–Thouless (BKT)
transition with quasi-long-range order [18]. Here, we work near criticality in
the normal phase of the BCS-type superconductor and in the regime of validity
of the Ginzburg-Levanyuk criterion, where the superconducting fluctuations are
sizable but weakly interacting [19, 20].

The disordered potential $V({\bm{x}})$ is assumed to be short-ranged and
Gaussian-distributed, with covariance $\langle
V({\bm{x}})V({\bm{x}}^{\prime})\rangle-\langle
V({\bm{x}})\rangle^{2}=g\delta({\bm{x}}-{\bm{x}}^{\prime})$ with $g>0$ 111More
generally, we can assume that the higher-order cumulants of the disordered
potential are irrelevant in an renormalization group sense.. We work far from
localization regimes, where the disorder can be treated in a classical
fashion, i.e. not accounting for coherent effects between scattering
trajectories. In practice, the impurity scattering is treated with the Born
approximation to second order in $g$ 222The terms at order $g$ can be absorbed
by a shift of the dispersion relation..

We formulate the kinetic theory of quantum information scrambling by starting
from the quantum kinetic equation on the many-world distribution functions
$F_{\alpha\beta}(t,{\bm{x}};\omega,{\bm{k}})$, where the indices
$\alpha,\beta\in\\{u,d\\}$ span two replicated worlds (up and down) [13, 17].
The intraworld components $\alpha=\beta$ correspond to the standard electronic
distribution functions. We concentrate on the interworld components
$F_{\alpha\neq\beta}$ which are directly related to four-point OTOCs and to
the growth of operators. In the gradient approximation, valid when the
microscopic scales set by $\hbar/E_{\rm F}$ and $1/k_{\rm F}$ are much shorter
than the spatiotemporal variations of $F_{\alpha\beta}$, the dynamics of the
latter are governed by a non-linear partial-integrodifferential equation
reading

$\displaystyle[\partial_{t}+\bm{v}_{\bm{k}}\cdot\bm{\nabla}_{\\!\\!\bm{x}}]F_{\alpha\beta}=I_{\alpha\beta},$
(2)

with the velocity
$\bm{v}_{\bm{k}}:=\bm{\nabla}_{\\!\\!\bm{k}}\epsilon_{\bm{k}}$ and where the
collision integral $I_{\alpha\beta}$ is a non-linear functional of the
$F_{\alpha\beta}$’s collecting contributions from the disordered potential and
the electronic interactions. The kinetic theory is considerably simplified by
working (i) in terms of the first two components of a partial-wave expansion
in ${\bm{k}}$, (ii) on-shell, i.e. $\omega\to\epsilon_{\bm{k}}$, and (iii)
near the Fermi surface, i.e. $k\to k_{\rm F}$. This amounts to working with
the ansatz

$\displaystyle
F_{\alpha\beta}(t,{\bm{x}};\omega,{\bm{k}})\mapsto\epsilon_{\alpha\beta}\left[\phi(t,{\bm{x}})+{\bm{u}}_{\bm{k}}\cdot\bm{\phi}_{1}(t,{\bm{x}})\right],$
(3)

with $\epsilon_{du}=-\epsilon_{ud}=1$,the unit vector
${\bm{u}}_{\bm{k}}:={\bm{k}}/k$, and where $\phi$ and $\bm{\phi}_{1}$ are the
isotropic and first anisotropic corrections, respectively, to the on-shell
interworld distribution function. The validity of the above ansatz is
discussed for the clean case in Ref. [17], and the partial-wave truncation is
all the more accurate in the presence of impurity scattering as it reduces
momentum anisotropy. Following the steps detailed in Appendix A, the kinetic
equation can be brought to a set of coupled non-linear PDEs reading

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{t}\phi+\frac{v_{\rm
F}}{d}\bm{\nabla}_{\\!\\!\bm{x}}\cdot\bm{\phi}_{1}\\!\\!\\!\\!&=\phi(\phi^{2}-1)/\tau_{\rm
sc}\\\\[6.0pt] \partial_{t}\bm{\phi}_{1}+v_{\rm
F}\bm{\nabla}_{\\!\\!\bm{x}}\phi\\!\\!\\!\\!&=\bm{\phi}_{1}(\gamma\phi^{2}-1)/\tau_{\rm
sc}-\bm{\phi}_{1}/\tau_{\rm el}.\end{array}\right.$ (6)

This represents an effective kinetic theory of scrambling in terms of the two
fields $\phi(t,{\bm{x}})$ and $\bm{\phi}_{1}(t,{\bm{x}})$. The dimensionless
parameter $\gamma$ tunes the distance to the superconducting transition:
$\gamma=1$ corresponds to criticality and $0<\gamma<1$ to off-critical regimes
in the normal phase. $\tau_{\rm el}\propto 1/g$ is the elastic timescale due
to scattering on the disordered potential and $\tau_{\rm sc}$ is the timescale
set by the inelastic scattering on the superconducting fluctuations
(Cooperons). These parameters of the model (6) depend in a non-trivial fashion
on those of the original microscopic model (Model.) and have to be understood
as renormalized quantities resulting from a complex cross-feed. For example,
disorder is known to enhance the inelastic scattering rate $1/\tau_{\rm sc}$
at low temperatures [23]. The PDEs (6) have a “correlated-world” solution,
$\phi=1$ and $\bm{\phi}_{1}=0$, corresponding to both replicated worlds
evolving coherently and is expected to be unstable for chaotic systems.
$\phi=\bm{\phi}_{1}=0$ is the “uncorrelated-world” solution, corresponding to
a total loss of coherence between worlds.

At $\tau_{\rm el}\to\infty$, we recover the clean model studied in Ref. [17].
In the non-interacting limit, $\tau_{\rm sc}\to\infty$, and in the diffusive
regime (at late times) where $\partial_{t}\bm{\phi}_{1}$ can be neglected
relative to $\bm{\phi}_{1}/\tau_{\rm el}$, the above coupled PDEs reduce to a
simple diffusion equation: $\partial_{t}\phi-D_{\rm
el}\bm{\nabla}_{\\!\\!\bm{x}}^{2}\phi=0$ with the elastic diffusion
coefficient $D_{\rm el}:={v_{\rm F}^{2}\tau_{\rm el}}/d$. In this Drude limit,
the correlated-world solution $\phi=1$ is stable against local perturbations,
expressing the absence of quantum information scrambling when only disorder is
present. In the generic case with both elastic and inelastic scattering, we
study how the scrambling dynamics evolve as the dimensionless disorder
coupling constant $\kappa:=\tau_{\rm sc}/\tau_{\rm el}\geq 0$ is increased
from the clean to the diffusive metal 333In particular, we do not start from
the Usadel equation as it relies on a diffusive approximation, see Ref. [43].
The analysis is simplified by rescaling time and space, $\tau:=t/\tau_{\rm
sc}\mbox{ and }{\bm{X}}:={\bm{x}}/\ell_{\rm sc}$ with $\ell_{\rm sc}:=v_{\rm
F}\tau_{\rm sc}/\sqrt{d}$, together with
$\bm{\phi}_{1}\mapsto\sqrt{d}\,\bm{\phi}_{1}$. The above PDEs become

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\phi+\bm{\nabla}_{\\!\\!\bm{X}}\cdot\bm{\phi}_{1}\\!\\!\\!\\!&=\phi(\phi^{2}-1)\\\\[6.0pt]
\partial_{\tau}\bm{\phi}_{1}+\bm{\nabla}_{\\!\\!\bm{X}}\phi\\!\\!\\!\\!&=\bm{\phi}_{1}(\gamma\phi^{2}-1)-\kappa\bm{\phi}_{1}.\end{array}\right.$
(9)

The initial conditions are taken as local and spherically symmetric
perturbations to the correlated-world solution:

$\displaystyle\left\\{\begin{array}[]{rl}\phi(\tau=0,{\bm{X}})&=1-\delta\phi_{0}(X)\\\
\bm{\phi}_{1}(\tau=0,{\bm{X}})&=0,\end{array}\right.$ (12)

with the radial coordinate $X:=\|{\bm{X}}\|$ and a perturbation
$0\leq\delta\phi_{0}(X)\ll 1$ which is non-vanishing on small support of
radius $R_{0}$ and $\delta\phi_{0}(X>R_{0})=0$. This guarantees that
subsequent dynamics reduce to an effective one-dimensional problem for
$\phi(\tau,X)$ and $\phi_{1}(\tau,X)$ along the radial direction. The late-
time dynamics marginally depend on the choice of $\delta\phi_{0}(X)$.

The instability of the correlated-world solution against local perturbations
is expected to generate a transient state where both solutions ($\phi\simeq 1$
and $\phi\simeq 0$) are separated by a domain wall, a front, located on a
sphere of growing radius. The profile and the motion of this front determine
the dynamics of the scrambling of quantum information. At late times, the
front is located far from the origin and it is governed by the $d=1$ version
of the PDEs (9) where $X$ is now the radial coordinate and $\phi_{1}$ is the
radial component of $\bm{\phi}_{1}$ [17]. In Appendix B, we provide rigorous
proof that this growth is bounded by a maximal velocity set by the Fermi
velocity, $v_{\rm lc}=1$ in units of $v_{\rm F}/\sqrt{d}$. This ensures that
the dynamics strictly take place within a causal light cone where $v_{\rm lc}$
acts as the effective speed of light. This motivates us to look for traveling
fronts propagating at a constant velocity $v\leq v_{\rm lc}$ and located at
$m_{\tau}\sim v\tau$ by assuming

$\displaystyle\left\\{\begin{array}[]{rl}\phi(\tau,m_{\tau}+z)&\stackrel{{\scriptstyle\tau\to\infty}}{{\longrightarrow}}f(z)\\\
\phi_{1}(\tau,m_{\tau}+z)&\stackrel{{\scriptstyle\tau\to\infty}}{{\longrightarrow}}f_{1}(z).\end{array}\right.$
(15)

##### Fermi shock-wave dynamics.

$0$$0.5$$1$$-8$$-6$$-4$$-2$$0$$2$$\kappa\leq\kappa^{*}$

$f$

$z=(x-v_{\rm lc}t)/\ell_{\rm
sc}-R_{0}$$\kappa=0$$\kappa=\kappa^{*}/2$$\kappa=\kappa^{*}$ Figure 2: Radial
profiles $f(z)$ of the late-time Fermi shock wave traveling at the maximal
velocity $v_{\rm lc}$, computed exactly from the coupled PDEs (9) for
$\kappa=0$, $\kappa^{*}/2$ and $\kappa^{*}$ ($\gamma=1$ i.e. $\kappa^{*}=2$.)
The offset $R_{0}$ is the radius of the initial condition. The front thickness
is of the order of $\ell_{\rm sc}:=v_{\rm F}\tau_{\rm sc}/\sqrt{d}$ in the
original units.

In the spirit of Ref. [17], we first look for traveling-wave solutions
propagating at the maximum velocity $v=v_{\rm lc}=1$ and that are
discontinuous at the boundary of the light cone. We leave the full computation
of the front profile to Appendix B. Here, we simply focus on extracting the
discontinuity by parametrizing the near-front geometry as

$\displaystyle\left\\{\begin{array}[]{rlrl}f(z<0,|z|\ll 1)&\simeq
L+z/\xi,&f(z>0)&=1\\\ f_{1}(z<0,|z|\ll 1)&\simeq-
M-z/\xi_{1},&f_{1}(z>0)&=0,\end{array}\right.$

where $L$, $M$, $\xi$, and $\xi_{1}$ are positive parameters to be determined,
and $m_{\tau}=\tau+R_{0}$. We find a critical disorder strength

$\displaystyle\kappa^{*}=1+\gamma$ (16)

separating two distinct solutions. For $\kappa<\kappa^{*}$, we find a solution
of Eqs. (9) with a traveling discontinuity from $f(0^{-})=L$ to $f(0^{+})=1$,
with

$\displaystyle
L(\kappa<\kappa^{*})=\frac{\sqrt{1+4\kappa^{*}(1+\kappa)}-1}{2\kappa^{*}}$
(17)

and $M=1-L$. We find finite values of $\xi$, indicating that the thickness of
the front is controlled by $\ell_{\rm sc}$. This generalizes the Fermi shock-
wave dynamics identified in Ref. [17] to the weakly-disordered case. We
illustrate this shock-wave profile in Fig. 2. When $\kappa\to\kappa^{*}$ from
below, the discontinuity closes continuously, $L\to 1$, but the slope remains
discontinuous. For $\kappa>\kappa^{*}$, $L=1$ is the only solution and the
finite slope left of the front brutally vanishes, $f^{\prime}(0)=0$, signaling
the sudden death of the Fermi shock wave.

##### FKKP dynamics.

Inspired by the FKPP equation proposed in Ref. [13], we look for smooth
traveling-wave solutions that belong to the FKPP class, with exponential tails
ahead of the front of the form

$\displaystyle\left\\{\begin{array}[]{rl}f(z\gg 1)&\simeq 1-Az\exp(-\mu z)\\\
f_{1}(z\gg 1)&\simeq-Bz\exp(-\mu z),\end{array}\right.$ (20)

where $A$, $B$, and the spatial decay rate $\mu$ are positive parameters. To
reveal the hidden FKPP nature of the PDEs (9), we reformulate them in terms of
$\delta\phi:=1-\phi$ and $\bm{\phi}_{1}$ which are expected to be small and
slowly varying ahead of the front, as per Eq. (20). Standard algebra detailed
in Appendix D yields, for any $\kappa\geq 0$,

$\left[\partial_{\tau}^{2}+(\kappa\\!-\\!\kappa^{*})\partial_{\tau}-\bm{\nabla}_{\\!\\!\bm{X}}^{2}\right]\\!\delta\phi\\!=\\!2(2\\!+\\!\kappa\\!-\\!\kappa^{*})\,\delta\phi-\mathrm{NL},$
(21)

where we collected the non-linear terms under the symbol $\mathrm{NL}$. The
term in $\partial_{\tau}\delta\phi$ is odd under time reversal and its
prefactor changes sign at $\kappa=\kappa^{*}$. When $\kappa>\kappa^{*}$, it
can be loosely interpreted as a dissipative term, bringing Eq. (21) to a
standard FKPP fashion, whereas it acts as a drive when $\kappa<\kappa^{*}$.
The term $\partial^{2}_{\tau}\delta\phi$ can be interpreted as inertia which,
to the best of our knowledge, has not been discussed in the broader context of
FKPP. Following a standard FKPP analysis detailed in App. D.2, we find that
there are no traveling-wave solutions of the form (20) when
$\kappa<\kappa^{*}$. However, for $\kappa>\kappa^{*}$, we now find a front
propagating at the velocity, in units of $v_{\rm F}/\sqrt{d}$,

$\displaystyle v(\kappa>\kappa^{*})$
$\displaystyle=2\sqrt{2}\frac{\sqrt{2+\kappa-\kappa^{*}}}{4+\kappa-\kappa^{*}}\leq
v_{\rm lc}.$ (22)

In the limit $\kappa\to\kappa^{*}$, we recover $v\to v_{\rm lc}$, consistently
with the Fermi shock wave. $v$ decreases monotonously with increasing disorder
strength, and $v\propto 1/\sqrt{\kappa}\to 0$ when $\kappa\to\infty$. In Fig.
1, we compare the front velocities extracted from the numerical solutions of
the PDEs (9) to the FKPP prediction in Eq. (22). The agrement is excellent. In
App. D.4, using results of Refs. [26, 27, 28, 29, 30, 31, 32, 33, 34, 35], we
provide a meticulous numerical analysis that further demonstrates, beyond any
reasonable doubt, the FKPP nature of the front dynamics as soon as
$\kappa>\kappa^{*}$.

Interestingly, the early-time dynamics governed by the linearized version of
Eq. (21) are characterized by an exponential growth of the spatially-
integrated perturbation $M(\tau):=\int\\!{\rm
d}{\bm{X}}\delta\phi(\tau,{\bm{X}})=M(0)\exp(2\tau)$ where the Lyapunov
exponent $\lambda=2$ does not depend on whether the shock-wave or the FKPP
regime is governing the scrambling dynamics. In the FKPP regime, the growth of
$\delta\phi$ ahead of the front is also exponential, as per Eq. (20), and the
scrambling dynamics can be characterized by a rate $\lambda_{\rm
FKPP}=v\mu=4+8/(\kappa-\kappa^{*})$ where $\mu$ is computed in App. D.2.

##### Diffusive regime.

$0$$0.5$$1$$-10$$0$$10$$\kappa\gg\kappa^{*}$

$\hat{f}$

$y=(x-m_{t})/\ell_{\rm diff}$ Figure 3: Radial profile of the late-time FKPP
front at $\kappa\gg 1$, solution of Eq. (23) traveling at the butterfly
velocity $v$ given in Eq. (24). $m_{t}\sim vt$ is the location of the front.
The dashed line is the numerical result obtained by solving Eqs. (9) with
$\kappa=60$ and $\gamma=1$ up to time $\tau=30$. The front thickness is of the
order of $\ell_{\rm diff}:=v_{\rm F}\sqrt{\tau_{\rm sc}\tau_{\rm
el}/2d}\ll\ell_{\rm sc}$.

We now delve into the overdamped regime $\kappa\gg 1$ where the inertial term
of Eq. (21) may be neglected. Keeping only the leading-order terms the right-
hand side (see the details in App. D.3), the scrambling dynamics are now
governed by, back in terms of $\phi$ and in the original units,

$\displaystyle\partial_{t}\phi-D_{\rm
el}\bm{\nabla}_{\\!\\!\bm{x}}^{2}\phi=\phi(\phi^{2}-1)/\tau_{\rm sc},$ (23)

where $1/\tau_{\rm sc}$ is the inelastic scattering rate and $D_{\rm
el}:=v_{\rm F}^{2}\tau_{\rm el}/d$ is the elastic diffusion coefficient. The
dependence on both $\bm{\phi}_{1}$ and $\gamma$ has dropped. This FKPP
equation is the non-integrable Newell-Whitehead equation which was first
studied in the context of non-linear fluid mechanics [36]. Similar equations
with diffusive terms were recently put forward to describe the dynamics of
scrambling of quantum information [13, 10, 37, 38, 39]. The traveling-wave
solutions of Eq. (23) propagate at a butterfly velocity, in the original
units,

$\displaystyle v=2\sqrt{\frac{2}{d}}\sqrt{\frac{\tau_{\rm el}}{\tau_{\rm
sc}}}v_{\rm F}\ll v_{{\rm F}}.$ (24)

Incidentally, this yields the relation

$\displaystyle v^{2}/\lambda_{\rm FKPP}=2D_{\rm el}$ (25)

which concretely connects scrambling dynamics on the one side to a measurable
transport quantity on the other side [40, 41, 9, 42]. The front profile is
conveniently computed by now measuring space in units of
$l_{\text{diff}}:=v_{\rm F}\sqrt{\tau_{\rm sc}\tau_{\rm el}/2d}\ll\ell_{\rm
sc}$. The rescaled front profile $\hat{f}(y):=f(y\,l_{\text{diff}})$ is the
solution to $2\hat{f}^{\prime\prime}+4\hat{f}^{\prime}=\hat{f}(1-\hat{f}^{2})$
with $\hat{f}(-\infty)=0$, $\hat{f}(\infty)=1$, and we can require
$\hat{f}(0)=1/2$. It is a smooth monotonous function that we represent in Fig.
3. The agreement with the numerical solutions of the PDEs (9) computed in the
diffusive regime ($\kappa=60$) is excellent.

##### Discussion.

In the Fermi shock-wave regime, $\kappa<\kappa^{*}$, it is still to be
clarified whether the scrambling front discontinuity could be smoothened by
corrections to the gradient approximation. In the FKPP regime,
$\kappa>\kappa^{*}$, we found that disorder can considerably reduce the
butterfly velocity, corroborating the results of Refs. [15, 14]: a realistic
ratio at room temperatures $\kappa\approx 10^{4}$ yields $v\sim 10^{-2}v_{\rm
F}\sim 10^{4}$ m/s, on par with the typical phonon-mediated sound velocity in
metals. Contrary to the shock-wave velocity, it cannot strictly be interpreted
a (slower) effective speed of light since the small tail ahead of the front,
while providing a quantum-chaotic exponential regime controlled by the
inelastic scattering rate, surreptitiously undermines the causality of the
light-cone structure. In both the Fermi shock-wave and the FKPP regimes, we
found a sharp butterfly front. This is similar to what has been reported for
models with a large local Hibert space such as the $O(N)$ or the Sachdev-Ye-
Kitaev (SYK) models [10, 9], but different from the diffusively broadening
fronts obtained for random quantum circuits [3, 4]. The precise conditions
under which strong quantum fluctuations and strong disorder could generate
relevant perturbations to the FKPP dynamics are still to be elucidated.

##### Acknowledgments.

We would like to thank Alex Kamenev for his initial inspiration which was
crucial in the development of this research project. We are grateful to Giulio
Biroli and Aavishkar A. Patel for insightful discussions. C.A. acknowledges
the support from the French ANR “MoMA” Project No. ANR-19-CE30-0020 and the
support from the Project No. 6004-1 of the Indo-French Centre for the
Promotion of Advanced Research (IFCPAR). This work received support from the
U.S. National Science Foundation Grant No. NSF-DMR 2018358 (A.M.).

## References

  * Bertini  _et al._ [2019] B. Bertini, P. Kos, and T. c. v. Prosen, Exact correlation functions for dual-unitary lattice models in $1+1$ dimensions, Phys. Rev. Lett. 123, 210601 (2019).
  * Claeys and Lamacraft [2020] P. W. Claeys and A. Lamacraft, Maximum velocity quantum circuits, Physical Review Research 2, 033032 (2020).
  * Nahum  _et al._ [2018] A. Nahum, S. Vijay, and J. Haah, Operator spreading in random unitary circuits, Physical Review X 8, 021014 (2018).
  * Bertini and Piroli [2020] B. Bertini and L. Piroli, Scrambling in random unitary circuits: Exact results, Physical Review B 102, 064305 (2020).
  * Das  _et al._ [2018] A. Das, S. Chakrabarty, A. Dhar, A. Kundu, D. A. Huse, R. Moessner, S. S. Ray, and S. Bhattacharjee, Light-cone spreading of perturbations and the butterfly effect in a classical spin chain, Physical Review Letters 121, 024101 (2018).
  * Khemani  _et al._ [2018] V. Khemani, D. A. Huse, and A. Nahum, Velocity-dependent lyapunov exponents in many-body quantum, semiclassical, and classical chaos, Physical Review B 98, 144304 (2018).
  * Ruidas and Banerjee [2021] S. Ruidas and S. Banerjee, Many-body chaos and anomalous diffusion across thermal phase transitions in two dimensions, SciPost Physics 11, 087 (2021).
  * Roy  _et al._ [2023] D. Roy, D. A. Huse, and M. Kulkarni, Out-of-time-ordered correlator in the one-dimensional kuramoto-sivashinsky and kardar-parisi-zhang equations, Phys. Rev. E 108, 054112 (2023).
  * Gu  _et al._ [2017a] Y. Gu, X.-L. Qi, and D. Stanford, Local criticality, diffusion and chaos in generalized Sachdev-Ye-Kitaev models, Journal of High Energy Physics 2017, 125 (2017a).
  * Chowdhury and Swingle [2017] D. Chowdhury and B. Swingle, Onset of many-body chaos in the $O(N)$ model, Physical Review D 96, 065005 (2017).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, Black holes and the butterfly effect, Journal of High Energy Physics 2014, 67 (2014).
  * Roberts and Swingle [2016] D. A. Roberts and B. Swingle, Lieb-robinson bound and the butterfly effect in quantum field theories, Physical Review Letters 117, 091602 (2016).
  * Aleiner  _et al._ [2016] I. L. Aleiner, L. Faoro, and L. B. Ioffe, Microscopic model of quantum butterfly effect: Out-of-time-order correlators and traveling combustion waves, Annals of Physics 375, 378 (2016).
  * Patel  _et al._ [2017] A. A. Patel, D. Chowdhury, S. Sachdev, and B. Swingle, Quantum butterfly effect in weakly interacting diffusive metals, Physical Review X 7, 031047 (2017).
  * Swingle and Chowdhury [2017] B. Swingle and D. Chowdhury, Slow scrambling in disordered quantum systems, Physical Review B 95, 060201 (2017).
  * Liao and Galitski [2018] Y. Liao and V. Galitski, Nonlinear sigma model approach to many-body quantum chaos: Regularized and unregularized out-of-time-ordered correlators, Phys. Rev. B 98, 205124 (2018).
  * Aron  _et al._ [2023] C. Aron, E. Brunet, and A. Mitra, Traveling discontinuity at the quantum butterfly front, SciPost Phys. 15, 042 (2023).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Larkin and Varlamov [2008] A. I. Larkin and A. A. Varlamov, Fluctuation phenomena in superconductors, in  _Superconductivity: Conventional and Unconventional Superconductors_, edited by K. H. Bennemann and J. B. Ketterson (Springer Berlin Heidelberg, Berlin, Heidelberg, 2008) pp. 369–458.
  * Lemonik and Mitra [2018] Y. Lemonik and A. Mitra, Quench dynamics of superconducting fluctuations and optical conductivity in a disordered system, Physical Review B 98, 10.1103/physrevb.98.214514 (2018).
  * Note [1] More generally, we can assume that the higher-order cumulants of the disordered potential are irrelevant in an renormalization group sense.
  * Note [2] The terms at order $g$ can be absorbed by a shift of the dispersion relation.
  * Altshuler and Aronov [1985] B. Altshuler and A. Aronov, in  _Electron–Electron Interactions in Disordered Systems_, Modern Problems in Condensed Matter Sciences, Vol. 10, edited by A. Efros and M. Pollack (Elsevier, 1985) pp. 1–153.
  * Note [3] In particular, we do not start from the Usadel equation as it relies on a diffusive approximation, see Ref. [43].
  * [25] See the supplementary material.
  * van Saarloos [2003] W. van Saarloos, Front propagation into unstable states, Physics Reports 386, 29 (2003).
  * Bramson [1978] M. D. Bramson, Maximal displacement of branching Brownian motion, Communications on Pure and Applied Mathematics 31, 531 (1978).
  * Roberts [2013] M. I. Roberts, A simple path to asymptotics for the frontier of a branching Brownian motion, The Annals of Probability 41, 3518 (2013).
  * Ebert and van Saarloos [2000] U. Ebert and W. van Saarloos, Front propagation into unstable states: universal algebraic convergence towards uniformly translating pulled fronts, Physica D 146, 1 (2000).
  * Berestycki  _et al._ [2016] J. Berestycki, E. Brunet, S. C. Harris, and M. Roberts, Vanishing corrections for the position in a linear model of FKPP fronts, Communications in Mathematical Physics 349, 857 (2016).
  * Nolen  _et al._ [2019] J. Nolen, J.-M. Roquejoffre, and L. Ryzhik, Refined long-time asymptotics for Fisher-KPP fronts, Communications in Contemporary Mathematics 21, 1850072 (2019).
  * Berestycki  _et al._ [2017] J. Berestycki, E. Brunet, and B. Derrida, Exact solution and precise asymptotics of a Fisher-KPP type front, Journal of Physics A: Mathematical and Theoretical 51, 035204 (2017).
  * Berestycki  _et al._ [2018] J. Berestycki, E. Brunet, and B. Derrida, A new approach to computing the asymptotics of the position of Fisher-KPP fronts, Europhysics Letters 122, 10001 (2018).
  * Graham [2019] C. Graham, Precise asymptotics for Fisher–KPP fronts, Nonlinearity 32, 1967 (2019).
  * Brunet and Derrida [1997] E. Brunet and B. Derrida, Shift in the velocity of a front due to a cutoff, Physical Review E 56, 2597 (1997).
  * Newell and Whitehead [1969] A. C. Newell and J. A. Whitehead, Finite bandwidth, finite amplitude convection, Journal of Fluid Mechanics 38, 279 (1969).
  * Grozdanov  _et al._ [2019] S. Grozdanov, K. Schalm, and V. Scopelliti, Kinetic theory for classical and quantum many-body chaos, Physical Review E 99, 012206 (2019).
  * Xu and Swingle [2019] S. Xu and B. Swingle, Locality, quantum fluctuations, and scrambling, Physical Review X 9, 031048 (2019).
  * Zhou  _et al._ [2023] T. Zhou, A. Guo, S. Xu, X. Chen, and B. Swingle, Hydrodynamic theory of scrambling in chaotic long-range interacting systems, Physical Review B 107, 014201 (2023).
  * Gu  _et al._ [2017b] Y. Gu, A. Lucas, and X.-L. Qi, Energy diffusion and the butterfly effect in inhomogeneous Sachdev-Ye-Kitaev chains, SciPost Physics 2, 018 (2017b).
  * Blake  _et al._ [2017] M. Blake, R. A. Davison, and S. Sachdev, Thermal diffusivity and chaos in metals without quasiparticles, Physical Review D 96, 106008 (2017).
  * Bilitewski  _et al._ [2018] T. Bilitewski, S. Bhattacharjee, and R. Moessner, Temperature dependence of the butterfly effect in a classical many-body system, Physical Review Letters 121, 250602 (2018).
  * Kamenev [2023] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_, 2nd ed. (Cambridge University Press, 2023).
  * Rammer [2007] J. Rammer,  _Quantum Field Theory of Non-equilibrium States_ (Cambridge University Press, 2007).
  * Fisher [1937] R. A. Fisher, The wave of advance of advantageous genes, Annals of Eugenics 7, 355 (1937).
  * Kolmogorov  _et al._ [1937] A. Kolmogorov, I. Petrovskii, and N. Piscunov, A study of the equation of diffusion with increase in the quantity of matter, and its application to a biological problem, Bulletin of Moscow State University Series A: Mathematics and Mechanics 1, 1 (1937).
  * Brunet [2016] É. Brunet,  _Some aspects of the Fisher-KPP equation and the branching Brownian motion_ , Habilitation à Diriger des Recherches, UPMC (2016).
  * Bramson [1983] M. D. Bramson, Convergence of solutions of the Kolmogorov equation to travelling waves, Memoirs of the American Mathematical Society 44, 10.1090/memo/0285 (1983).
  * Hamel  _et al._ [2013] F. Hamel, J. Nolen, J.-M. Roquejoffre, and L. Ryzhik, A short proof of the logarithmic Bramson correction in Fisher-KPP equations, Networks and Heterogeneous Media 8, 275 (2013).

##  Appendix A Kinetic equations of scrambling in a disordered interacting
metal

Figure 4: Two-world Keldysh contour $\mathcal{C}$: the theory is replicated
into an “up” world dynamics, marked by the index $u$, and the “down” world
dynamics marked by $d$. In blue: the insertion points of the OTOC
$\mathcal{A}(t,{\bm{x}};t^{\prime},{\bm{x}}^{\prime}):=\langle\mathrm{T}_{\mathcal{C}}\,\psi_{d}^{-\dagger}(0,0)\psi_{d}^{+}(t,{\bm{x}})\psi_{u}^{-}(0,0)\psi_{u}^{+\dagger}(t^{\prime},{\bm{x}}^{\prime})\rangle$.

In this Appendix, we derive the set of kinetic equations (4) that govern the
dynamics of quantum information scrambling in the disordered interacting metal
given by the Hamiltonian in Eq. (1). More information on the procedure for the
clean case can be found in Ref. [17].

###  A.1 Many-world formalism

We work within the many-world formalism, sometimes referred to as the
augmented Keldysh formalism. We refer the reader to Ref. [13] for a
pedagogical introduction, see also Ref. [43]. Two-point correlators are
defined as

$\displaystyle{\rm
i}G_{\alpha\beta}^{ab}(t,{\bm{x}};t^{\prime},{\bm{x}}^{\prime})$
$\displaystyle:=\langle\mathrm{T}_{\mathcal{C}}\psi_{\alpha}^{a}(t,{\bm{x}})\psi_{\beta}^{b\dagger}(t^{\prime},{\bm{x}}^{\prime})\rangle,$
(26)

with the Keldysh indices $a,b=+,-$, and the world indices $\alpha,\beta=u,d$.
We set $\hbar=1$. $\mathrm{T}_{\mathcal{C}}$ is the time-ordering operator on
the two-world contour $\mathcal{C}$ represented in Fig. 4. The fermionic
operators $\psi$ and $\psi^{\dagger}$ are written in the Heisenberg picture,
and $\langle\ldots\rangle:=\mathrm{Tr}\left[\ldots\rho_{0}\right]$ where
$\rho_{0}$ is the initial density matrix with $\mathrm{Tr}\,\rho_{0}=1$. The
system is initially prepared in equilibrium at the temperature $T$ with the
Gibbs state $\rho_{0}\sim{\rm e}^{-H/T}$ where we set $k_{\rm B}=1$. The
equilibrium dynamics ensure that, in the absence of an external perturbation
explicitly breaking this invariance, the intraworld dynamics is translational
invariant in space and time, and one may work in Fourier space. This is not
the case for interworld quantities ($\alpha\neq\beta$) since we are precisely
probing their response with respect to local perturbations. Therefore, we use
the Wigner representation for the interworld components, e.g.

$\displaystyle
G^{ab}_{\alpha\beta}(\omega,{\bm{k}};t,{\bm{x}}):=\int\\!\\!{\rm
d}{\bm{x}}^{\prime}\\!\\!\int\\!\\!{\rm d}t^{\prime}\,{\rm e}^{{\rm i}(\omega
t^{\prime}-{\bm{k}}\cdot{\bm{x}}^{\prime})}\,G^{ab}_{\alpha\beta}\left(t+\frac{t^{\prime}}{2},{\bm{x}}+\frac{{\bm{x}}^{\prime}}{2};t-\frac{t^{\prime}}{2},{\bm{x}}-\frac{{\bm{x}}^{\prime}}{2}\right).$
(27)

To alleviate the notations, we drop the intraworld indices, e.g.
$G_{uu}=G_{dd}=G$, and we reserve the double indices $\alpha\beta$ to
interworld quantities, i.e. by $\alpha\beta$ we mean $\alpha\neq\beta$
specifically. We work in the Keldysh basis with the retarded ($R$) and Keldysh
($K$) Green’s functions which, after using the properties of the contour
$\mathcal{C}$, read

$\displaystyle\left\\{\begin{array}[]{rl}G^{R}&=G^{++}-G^{+-}\\\
G^{K}&=G^{++}+G^{--}\end{array}\right.\mbox{ and
}\left\\{\begin{array}[]{rl}G^{R}_{\alpha\beta}&=0\\\
G^{K}_{\alpha\beta}&=2G_{\alpha\beta}^{-+}.\end{array}\right.$ (32)

The intraworld Schwinger-Dyson equations yield the standard (single-world)
theory

$\displaystyle\left\\{\begin{array}[]{rl}G^{R}(\omega,{\bm{k}})&=\left[\omega-\epsilon({\bm{k}})-\Sigma^{R}(\omega,{\bm{k}})\right]^{-1}\\\
G^{K}(\omega,{\bm{k}})&=|G^{R}(\omega,{\bm{k}})|^{2}\Sigma^{K}(\omega,{\bm{k}}).\end{array}\right.$
(35)

The interworld Schwinger-Dyson equations read

$\displaystyle\left\\{\begin{array}[]{rl}G_{\alpha\beta}^{R}(\omega,{\bm{k}};t,{\bm{x}})&=0\\\
G^{K}_{\alpha\beta}(\omega,{\bm{k}};t,{\bm{x}})&=G^{R}(\omega,{\bm{k}};t,{\bm{x}})\star\Sigma^{K}_{\alpha\beta}(\omega,{\bm{k}};t,{\bm{x}})\star
G^{R}(\omega,{\bm{k}};t,{\bm{x}})^{*},\end{array}\right.$ (38)

where we introduced the Moyal product $\star:=\exp\left[\frac{{\rm
i}}{2}(\overleftarrow{\partial_{\omega}}\overrightarrow{\partial_{t}}\\!-\\!\overleftarrow{\bm{\nabla}_{\\!\\!\bm{k}}}\cdot\overrightarrow{\bm{\nabla}_{\\!\\!\bm{x}}}\\!-\\!\overleftarrow{\partial_{t}}\overrightarrow{\partial_{\omega}}\\!+\\!\overleftarrow{\bm{\nabla}_{\\!\\!\bm{x}}}\cdot\overrightarrow{\bm{\nabla}_{\\!\\!\bm{k}}})\right]$
where the left (right) arrow designates a derivative operator acting on the
left (right) of the star symbol. The self-energies (the $\Sigma$’s) entering
the Schwinger-Dyson equations (35) and (38) will be discussed below in App.
A.2.

##### Connection to scrambling.

Let us briefly outline the connection between the interworld two-point Green’s
functions and the four-point OTOCs which are commonly used to probe the
dynamical signatures of quantum chaos. An example of such an OTOC is

$\displaystyle\mathcal{A}(t,{\bm{x}};t^{\prime},{\bm{x}}^{\prime}):=$
$\displaystyle\mathrm{Tr}\,\left[\psi^{\dagger}(0)\ {\rm e}^{{\rm
i}Ht}\psi({\bm{x}}){\rm e}^{-{\rm i}Ht}\ \psi(0)\ {\rm e}^{{\rm
i}Ht^{\prime}}\psi^{\dagger}({\bm{x}}^{\prime})\ {\rm e}^{-{\rm
i}Ht^{\prime}}\rho_{0}\right]$ (39) $\displaystyle=$
$\displaystyle\langle\mathrm{T}_{\mathcal{C}}\,\psi_{d}^{-\dagger}(0,0)\psi_{d}^{+}(t,{\bm{x}})\psi_{u}^{-}(0,0)\psi_{u}^{+\dagger}(t^{\prime},{\bm{x}}^{\prime})\rangle.$
(40)

At late times, we expect the decoupling

$\displaystyle\mathcal{A}(t,{\bm{x}};t^{\prime},{\bm{x}}^{\prime})\simeq
n_{0}\,{\rm i}G_{du}^{++}(t,{\bm{x}};t^{\prime},{\bm{x}}^{\prime})\propto
G^{K}_{du}(t,{\bm{x}};t^{\prime},{\bm{x}}^{\prime}),$ (41)

with $n_{0}:=\mathrm{Tr}\,\left[\psi^{\dagger}(0)\psi(0)\rho_{0}\right]$ and
where the interworld Green’s function is computed as the solution of the
interworld Schwinger-Dyson equations in the presence of a local perturbation
to the correlated world solution at ${\bm{x}}=0$ and time $t=0$.

###  A.2 Self-energy contributions

There are two contributions to the electronic self-energy: elastic scattering
on the disordered potential and inelastic scattering on the Cooperons. We
refer the reader to Chap. 8 in Ref. [44], specifically around Eqs (8.180) and
(8.181), for a detailed discussion of the self-energy contributions to the
(intraworld) kinetic equations in a very similar setting. We follow the same
approach for the interworld theory and

$\displaystyle\left\\{\begin{array}[]{rl}\Sigma^{R}&=\Sigma^{R}_{\rm
dis}+\Sigma^{R}_{U}\\\\[3.0pt] \Sigma^{K}&=\Sigma^{K}_{\rm
dis}+\Sigma^{K}_{U}\end{array}\right.\mbox{ and
}\left\\{\begin{array}[]{rl}\Sigma^{R}_{\alpha\beta}&=\Sigma^{R}\\\\[3.0pt]
\Sigma^{K}_{\alpha\beta}&=\Sigma^{K}_{{\rm
dis}\,\alpha\beta}+\Sigma^{K}_{U\,\alpha\beta},\end{array}\right.$ (46)

where the disorder contribution $\Sigma_{\rm dis}$ and the interacting
contribution $\Sigma_{U}$ are detailed below. It is to be noted that both
contributions to the self-energy feed each other. In particular, the
interacting contribution depends implicitly on the disorder strength.

##### Disorder contribution.

Expanding to the second order in the coupling to the disordered potential, we
obtain the self-energy contribution

$\displaystyle\left\\{\begin{array}[]{rl}\Sigma^{R/K}_{{\rm
dis}}(\omega,{\bm{k}})&=g\sum_{{\bm{k}}^{\prime}}G^{R/K}(\omega,{\bm{k}}^{\prime})\\\\[3.0pt]
\Sigma^{K}_{{\rm
dis}\,\alpha\beta}(\omega,{\bm{k}};t,{\bm{x}})&=g\sum_{{\bm{k}}^{\prime}}G^{K}_{\alpha\beta}(\omega,{\bm{k}}^{\prime};t,{\bm{x}}),\end{array}\right.$
(49)

where $g\delta(\bm{x}-{\bm{x}}^{\prime}):=\langle
V({\bm{x}})V({\bm{x}}^{\prime})\rangle-\langle V({\bm{x}})\rangle^{2}$.

##### Interaction contribution.

We treat the finite electronic interaction in the Cooper channel by means of
the random phase approximation (RPA). The intraworld self-energy contributions
($\alpha=\beta$) read

$\displaystyle\left\\{\begin{array}[]{rl}\,\mathrm{Im}\,\Sigma_{U}^{R}(\omega,{\bm{k}})\\!\\!\\!\\!&=\displaystyle-\frac{1}{2}\sum_{{\bm{k}}^{\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\bigg{\\{}\,\mathrm{Im}\,D^{R}(\omega^{\prime},{\bm{k}}^{\prime}){\rm
i}G^{K}(\omega^{\prime}-\omega,{\bm{k}}^{\prime}-{\bm{k}})-{\rm
i}D^{K}(\omega^{\prime},{\bm{k}}^{\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega,{\bm{k}}^{\prime}-{\bm{k}})\bigg{\\}}\\\\[6.0pt]
\Sigma_{\rm U}^{K}(\omega,{\bm{k}})\\!\\!\\!\\!&=2{\rm
i}\,F(\omega)\,\mathrm{Im}\,\Sigma^{R}(\omega,{\bm{k}}),\end{array}\right.$
(52)

where $F(\omega):=\tanh(\omega/2T)$ is dictated by the (fermionic)
fluctuation-dissipation theorem. Within the RPA in the particle-particle
channel, the Cooperon propagators are given by

$\displaystyle\left\\{\begin{array}[]{l}D^{R}(\omega,{\bm{k}})=\left[U^{-1}-\Pi^{R}(\omega,{\bm{k}})\right]^{-1}\\\\[6.0pt]
D^{K}(\omega,{\bm{k}})=2{\rm
i}\,P(\omega)\,\mathrm{Im}\,D^{R}(\omega,{\bm{k}}),\end{array}\right.$ (55)

where $P(\omega):=\coth(\omega/2T)$ is dictated by the (bosonic) fluctuation-
dissipation theorem and with the polarization bubbles read

$\displaystyle\left\\{\begin{array}[]{rl}\,\mathrm{Im}\,\Pi^{R}(\omega,{\bm{k}})\\!\\!\\!\\!&=\displaystyle\sum_{{\bm{k}}^{\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}{\rm
i}G^{K}(\omega^{\prime},{\bm{k}}^{\prime})\,\mathrm{Im}\,G^{R}(\omega-\omega^{\prime},{\bm{k}}-{\bm{k}}^{\prime})\\\\[6.0pt]
\Pi^{K}(\omega,{\bm{k}})\\!\\!\\!\\!&=2{\rm
i}\,P(\omega)\,\mathrm{Im}\,\Pi^{R}(\omega,{\bm{k}}),\end{array}\right.$ (58)

Close to the superconducting transition, the near-critical Cooperon propagator
reads

$\displaystyle D^{R}(\omega,{\bm{k}})\approx\frac{-1/\rho_{\rm F}}{r-{\rm
i}a{\omega}/{T}+\xi^{2}k^{2}+\ldots},$ (59)

where $\rho_{\rm F}$ is the density of states at the Fermi energy, and the
parameter $r\propto(T-T_{\rm c})/T_{\rm c}$ is the detuning from the critical
temperature $T_{\rm c}$. In addition, the remaining parameters are all
positive with $a\sim\mathcal{O}(1)$, $\xi^{2}\sim v_{\rm F}^{2}/T^{2}$, where
$v_{\rm F}$ is the Fermi velocity. At criticality $r\to 0$, the Cooperon
becomes soft with diverging length scale $l\sim 1/r^{\nu}$ and timescale $\sim
l^{z}$ (here $\nu=1/2$, $z=2$), and the propagator is singular at
$\omega=k=0$.

The interworld components read

$\displaystyle\Sigma_{U\,\alpha\beta}^{K}(\omega,{\bm{k}};t,{\bm{x}})$
$\displaystyle=-\frac{{\rm i}}{2}\sum_{{\bm{k}}^{\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}D_{\alpha\beta}^{K}(\omega^{\prime},{\bm{k}}^{\prime};t,{\bm{x}})G_{\beta\alpha}^{K}(\omega^{\prime}-\omega,{\bm{k}}^{\prime}-{\bm{k}};t,{\bm{x}}),$
(60)

with

$\displaystyle D_{\alpha\beta}^{K}(\omega,{\bm{k}};t,{\bm{x}})$
$\displaystyle=D^{R}(\omega,{\bm{k}})\star\Pi_{\alpha\beta}^{K}(\omega,{\bm{k}};t,{\bm{x}})\star
D^{R}(\omega,{\bm{k}})^{*},$ (61)
$\displaystyle\Pi_{\alpha\beta}^{K}(\omega,{\bm{k}};t,{\bm{x}})$
$\displaystyle=\frac{{\rm i}}{2}\sum_{{\bm{k}}^{\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}G_{\alpha\beta}^{K}(\omega^{\prime},{\bm{k}}^{\prime};t,{\bm{x}})G_{\alpha\beta}^{K}(\omega-\omega^{\prime},{\bm{k}}-{\bm{k}}^{\prime};t,{\bm{x}}).$
(62)

###  A.3 Interworld kinetics

We start from the quantum kinetic equation which is a rewrite of the
Schwinger-Dyson equations (38). See, e.g., a standard derivation in Ref [43].
The interworld Keldysh Green’s function ($\alpha\neq\beta$) is parametrized in
terms of the interworld distribution function $F_{\alpha\beta}$:

$\displaystyle
G^{K}_{\alpha\beta}(\omega,{\bm{k}};t,{\bm{x}})=G^{R}(\omega,{\bm{k}})\star
F_{\alpha\beta}(\omega,{\bm{k}};t,{\bm{x}})-F_{\alpha\beta}(\omega,{\bm{k}};t,{\bm{x}})\star
G^{A}(\omega,{\bm{k}}),$ (63)

To alleviate the notations, we now drop the explicit dependence of ${\bm{x}}$
and $t$ except when this obscures the meaning.

##### Gradient approximation.

Expanding the Moyal products to first non-trivial orders, the interworld
quantum kinetic equation reads

$\displaystyle\left[\partial_{t}+\bm{v}_{\bm{k}}\cdot\bm{\nabla}_{\\!\\!\bm{x}}\right]F_{\alpha\beta}(\omega,{\bm{k}})$
$\displaystyle=I_{\alpha\beta}(\omega,{\bm{k}}),$ (64)

with the collision integral

$\displaystyle I_{\alpha\beta}(\omega,{\bm{k}})$ $\displaystyle=2\,{\rm
Im}\Sigma^{R}(\omega,{\bm{k}})\,F_{\alpha\beta}(\omega,{\bm{k}})+{\rm
i}\Sigma^{K}_{\alpha\beta}(\omega,{\bm{k}}).$ (65)

Using the expressions for the self-energies in App. A.2, we have

$\displaystyle I_{\alpha\beta}(\omega,{\bm{k}})$
$\displaystyle=2\sum_{{\bm{k}}^{\prime}}g({\bm{k}}-{\bm{k}}^{\prime})\mathrm{Im}\,G^{R}(\omega,{\bm{k}}^{\prime})\left[F_{\alpha\beta}(\omega,{\bm{k}})-F_{\alpha\beta}(\omega,{\bm{k}}^{\prime})\right]$
$\displaystyle\quad+2\sum_{{\bm{k}}^{\prime}{\bm{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\bm{k}}^{\prime})|^{2}$
$\displaystyle\quad\qquad\times\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\bm{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\bm{k}}^{\prime}-{\bm{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega,{\bm{k}}^{\prime}-{\bm{k}})$
$\displaystyle\quad\qquad\times\bigg{\\{}\left[\tanh\left(\frac{\omega^{\prime\prime}}{2T}\right)+\tanh\left(\frac{\omega^{\prime}-\omega^{\prime\prime}}{2T}\right)\right]\left[\coth\left(\frac{\omega^{\prime}}{2T}\right)+\tanh\left(\frac{\omega-\omega^{\prime}}{2T}\right)\right]F_{\alpha\beta}(\omega,{\bm{k}})$
$\displaystyle\quad\qquad\qquad+F_{\alpha\beta}(\omega^{\prime\prime},{\bm{k}}^{\prime\prime})F_{\alpha\beta}(\omega^{\prime}-\omega^{\prime\prime},{\bm{k}}^{\prime}-{\bm{k}}^{\prime\prime})F_{\beta\alpha}(\omega^{\prime}-\omega,{\bm{k}}^{\prime}-{\bm{k}})\bigg{\\}}.$
(66)

One may check that the collision integral vanishes identically at the
uncorrelated world solution

$\displaystyle F^{\rm uncorr}_{\alpha\beta}(\omega,{\bm{k}})=0,$ (67)

as well as at the correlated world solution

$\displaystyle F^{\rm
corr}_{\alpha\beta}(\omega,{\bm{k}})=\epsilon_{\alpha\beta}+\tanh\left(\omega/2T\right),$
(68)

with $\epsilon_{du}=1$ and $\epsilon_{ud}=-1$.

##### Quasi-particle approximation.

In all practical instances, $F_{\alpha\beta}(\omega,{\bm{k}})$ appears
multiplied by the density of states $\,\mathrm{Im}\,G^{R}(\omega,{\bm{k}})$,
see e.g. Eq. (66). When quasi-particles are well defined, with a dispersion
relation $\epsilon_{\bm{k}}$, the density of states is sharply peaked around
$\omega=\epsilon_{\bm{k}}$ and one may seamlessly exchange
$F_{\alpha\beta}(\omega,{\bm{k}})$ with the on-shell quasi-particle
distribution function
$\tilde{F}_{\alpha\beta}({\bm{k}}):=F_{\alpha\beta}(\omega=\epsilon_{\bm{k}},{\bm{k}})$.
From now on, we use the tilde notation to denote the on-shell prescription
$\omega=\epsilon_{\bm{k}}$.

##### Partial-wave ansatz.

The kinetic equation in (64), with its collision integral in Eq. (66), is
considerably simplified by working with the following ansatz which implements,
altogether, the truncation to the first terms of a partial-wave expansion in
the momentum space, the quasi-particle approximation, and the focusing on the
physics at the Fermi surface. Restoring the explicit dependence on ${\bm{x}}$
and $t$, it reads

$\displaystyle\tilde{F}^{\rm
ansatz}_{\alpha\beta}(t,{\bm{x}};{\bm{u}}_{\bm{k}})$
$\displaystyle=\left[\phi(t,{\bm{x}})+{\bm{u}}_{\bm{k}}\cdot\bm{\phi}_{1}(t,{\bm{x}})\right]\,\tilde{F}_{\alpha\beta}^{\rm
corr}({\bm{k}})\lvert_{k\to k_{\rm F}}$ (69)
$\displaystyle=\epsilon_{\alpha\beta}\left[\phi(t,{\bm{x}})+{\bm{u}}_{\bm{k}}\cdot\bm{\phi}_{1}(t,{\bm{x}})\right],$
(70)

with the unit vector ${\bm{u}}_{\bm{k}}:={\bm{k}}/k$ and where $\phi$ and
$\bm{\phi}_{1}$ are real scalar and vector fields coding, respectively, for
the isotropic component and the first anisotropic corrections to the
interworld distribution function
$\tilde{F}_{\alpha\beta}({\bm{k}})\big{\lvert}_{k\to k_{\rm F}}$. Conversely,
this amounts to performing a projection on the first two terms of a partial-
wave expansion:

$\displaystyle\left\\{\begin{array}[]{rl}\phi(t,{\bm{x}})\\!\\!\\!\\!&=\displaystyle\frac{1}{S_{d-1}}\epsilon_{\alpha\beta}\int{\rm
d}\Omega_{\bm{k}}F_{\alpha\beta}(t,{\bm{x}};\epsilon_{\bm{k}},{\bm{k}})\Big{\lvert}_{k\to
k_{\rm F}}\\\\[6.0pt]
\bm{\phi}_{1}(t,{\bm{x}})\\!\\!\\!\\!&=\displaystyle\frac{d}{S_{d-1}}\epsilon_{\alpha\beta}\int{\rm
d}\Omega_{\bm{k}}{\bm{u}}_{\bm{k}}F_{\alpha\beta}(t,{\bm{x}};\epsilon_{\bm{k}},{\bm{k}})\Big{\lvert}_{k\to
k_{\rm F}},\end{array}\right.$ (73)

where ${\rm d}\Omega_{\bm{k}}$ is the elementary solid angle in the direction
of ${\bm{k}}$, and $S_{d-1}:=\int{\rm d}\Omega_{\bm{k}}$ is the surface area
of the $d-1$-sphere with unit radius. One has $\int{\rm
d}\Omega_{\bm{k}}u_{\bm{k}}^{i}u_{\bm{k}}^{j}=$ $\delta_{ij}\,S_{d-1}/d$. From
Eq. (64), we obtain

$\displaystyle\left[\partial_{t}\phi+{\bm{u}}_{{\bm{k}}}\cdot\partial_{t}\bm{\phi}_{1}+v_{\rm
F}\left({\bm{u}}_{\bm{k}}\cdot\bm{\nabla}_{\\!\\!\bm{x}}\right)\left(\phi+{\bm{u}}_{\bm{k}}\cdot\bm{\phi}_{1}\right)\right]\epsilon_{\alpha\beta}=\tilde{I}_{{\rm
dis}\,\alpha\beta}({\bm{k}})\lvert_{k\to k_{\rm
F}}+\tilde{I}_{U\,\alpha\beta}({\bm{k}})\lvert_{k\to k_{\rm F}},$ (74)

with

$\displaystyle\tilde{I}_{{\rm dis}\,\alpha\beta}({\bm{k}})$
$\displaystyle=g\bm{\phi}_{1}\cdot\sum_{{\bm{k}}^{\prime}}\mathrm{Im}\,G^{R}(\epsilon_{\bm{k}},{\bm{k}}^{\prime})\left[{\bm{u}}_{\bm{k}}\tilde{F}_{\alpha\beta}^{\rm
corr}({\bm{k}})-{\bm{u}}_{{\bm{k}}^{\prime}}\tilde{F}_{\alpha\beta}^{\rm
corr}({\bm{k}}^{\prime})\right].$ (75)

At the Fermi surface, $k\to k_{\rm F}$, the term
$\mathrm{Im}\,G^{R}(\epsilon_{\bm{k}},{\bm{k}}^{\prime})$ selects momenta
${\bm{k}}^{\prime}$ which also live at the Fermi surface. We obtain

$\displaystyle\tilde{I}_{{\rm dis}\,\alpha\beta}({\bm{k}})\lvert_{k\to k_{\rm
F}}$ $\displaystyle=-\epsilon_{\alpha\beta}\pi g\rho_{\rm
F}\,\bm{\phi}_{1}\cdot{\bm{u}}_{\bm{k}}$ (76)
$\displaystyle=-\epsilon_{\alpha\beta}\bm{\phi}_{1}\cdot{\bm{u}}_{\bm{k}}/\tau_{\rm
el},$ (77)

where we made use of $\tilde{F}_{\alpha\beta}({\bm{k}})\lvert_{k\to k_{\rm
F}}=\epsilon_{\alpha\beta}$ with $\epsilon_{du}=-\epsilon_{ud}=1$, and
$\rho_{\rm F}$ is the total density of states at the Fermi energy. In the last
line, we introduced the elastic scattering rate

$\displaystyle\frac{1}{\tau_{\rm el}}:=\pi g\rho_{\rm F}.$ (78)

The use of the Born approximation is legitimate far from the localization
regimes, $1/{\tau_{\rm el}}\ll E_{{\rm F}}$. On the interaction side, using
the results of Ref. [17], we have

$\displaystyle\tilde{I}_{U\,\alpha\beta}({\bm{k}})\lvert_{k\to k_{\rm
F}}=\epsilon_{\alpha\beta}\phi(\phi^{2}-1)/\tau_{\rm
sc}+\epsilon_{\alpha\beta}(\bm{\phi}_{1}\cdot{\bm{u}}_{\bm{k}})(\gamma\phi^{2}-1)/\tau_{\rm
sc}+\ldots$ (79)

where the $\ldots$ stand for higher-order spherical harmonics and we
introduced the (temperature and disorder dependent) inelastic scattering rate

$\displaystyle\frac{1}{\tau_{\rm
sc}}:=-2\,\mathrm{Im}\,\tilde{\Sigma}^{R}(k_{\rm F}),$ (80)

with, explicitly,

$\displaystyle{\rm Im}\,\tilde{\Sigma}_{U}^{R}({\bm{k}})$
$\displaystyle=\sum_{{\bm{k}}^{\prime}{\bm{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\bm{k}}^{\prime})|^{2}\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\bm{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\bm{k}}^{\prime}-{\bm{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\bm{k}},{\bm{k}}^{\prime}-{\bm{k}})$
$\displaystyle\qquad\qquad\qquad\qquad\times\left[\tanh\left(\frac{\omega^{\prime\prime}}{2T}\right)+\tanh\left(\frac{\omega^{\prime}-\omega^{\prime\prime}}{2T}\right)\right]\left[\coth\left(\frac{\omega^{\prime}}{2T}\right)+\tanh\left(\frac{\epsilon_{\bm{k}}-\omega^{\prime}}{2T}\right)\right],$
(81)

and where the parameter $0<\gamma\leq 1$ quantifies the distance to
criticality ($\gamma=1$ at criticality).

We now project Eq. (74) onto the first two terms of the partial-wave
expansion, using

$\displaystyle\left\\{\begin{array}[]{rl}\int{\rm
d}\Omega_{\bm{k}}\tilde{I}_{{\rm dis}\,\alpha\beta}({\bm{k}})\lvert_{k\to
k_{\rm F}}&=0\\\\[6.0pt] \frac{d}{S_{d-1}}\int{\rm
d}\Omega_{\bm{k}}{\bm{u}}_{\bm{k}}\tilde{I}_{{\rm
dis}\,\alpha\beta}({\bm{k}})\lvert_{k\to k_{\rm
F}}&=-\epsilon_{\alpha\beta}\bm{\phi}_{1}/\tau_{\rm
el}\end{array}\right.\mbox{ and }\left\\{\begin{array}[]{rl}\int{\rm
d}\Omega_{\bm{k}}\tilde{I}_{U\,\alpha\beta}({\bm{k}})\lvert_{k\to k_{\rm
F}}&=\epsilon_{\alpha\beta}\phi(\phi^{2}-1)/\tau_{\rm sc}\\\\[6.0pt]
\frac{d}{S_{d-1}}\int{\rm
d}\Omega_{\bm{k}}{\bm{u}}_{\bm{k}}\tilde{I}_{U\,\alpha\beta}({\bm{k}})\lvert_{k\to
k_{\rm F}}&=\epsilon_{\alpha\beta}\bm{\phi}_{1}(\gamma\phi^{2}-1)/\tau_{\rm
sc}.\end{array}\right.$ (86)

Finally, we obtain the set of coupled PDEs (4) which govern the kinetics of
quantum information scrambling,

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{t}\phi+\frac{v_{\rm
F}}{d}\bm{\nabla}_{\\!\\!\bm{x}}\cdot\bm{\phi}_{1}\\!\\!\\!\\!&=\phi(\phi^{2}-1)/\tau_{\rm
sc}\\\\[6.0pt] \partial_{t}\bm{\phi}_{1}+v_{\rm
F}\bm{\nabla}_{\\!\\!\bm{x}}\phi\\!\\!\\!\\!&=\bm{\phi}_{1}(\gamma\phi^{2}-1)/\tau_{\rm
sc}-\bm{\phi}_{1}/\tau_{\rm el}.\end{array}\right.$ (89)

##  Appendix B Light-cone structure of scrambling dynamics

Here below, we carefully demonstrate that the solutions of the coupled PDEs
(5) in $d=1$ are strictly causal in the sense that they are bounded by a light
cone which propagates at velocity $v_{\rm lc}=1$ (i.e. $v_{\rm lc}=v_{\rm
F}/\sqrt{d}$ in the original units). In other words, for any $X$ and any
$\tau>0$, the values of $\phi(\tau,X)$ and $\phi_{1}(\tau,X)$ depend only on
the initial conditions $\phi_{0}(z)$ and ${\phi_{1}}_{0}(z)$ in the range
$z\in[X-\tau,X+\tau]$. Then, with initial conditions such that $\phi_{0}(X)=1$
and ${\phi_{1}}_{0}(X)=0$ for all $|X|>R_{0}$, as studied in the paper, the
solution strictly satisfies

$\phi(\tau,X)=1\quad\text{and}\quad\phi_{1}(\tau,X)=0\quad\text{for
{\hbox{|X|>R_0+\tau}}}.$ (90)

Let us first recall the system of equations (5) in dimension $d=1$:

$\left\\{\begin{array}[]{rl}\partial_{\tau}\phi+\partial_{X}\phi_{1}&=\phi(\phi^{2}-1)\\\
\partial_{\tau}\phi_{1}+\partial_{X}\phi&=\phi_{1}\left(\gamma\phi^{2}-1-\kappa\right).\end{array}\right.$
(91)

where $0\leq\gamma\leq 1$ and $\kappa\geq 0$. More generally, the light cone
property can be proved to hold for any system of the form
$\partial_{\tau}\phi+\partial_{X}\phi_{1}=f(\phi,\phi_{1})$,
$\partial_{\tau}\phi_{1}+\partial_{X}\phi=g(\phi,\phi_{1})$ with “reasonable”
(say, polynomial) functions $f$ and $g$, and for any bounded initial
condition. However, in order to simplify the discussion, we only consider the
specific case (91) and initial conditions such that

$0\leq\phi_{0}(X)\leq 1,\qquad{\phi_{1}}_{0}(X)=0.$ (92)

For such an initial condition, we also show that the solution satisfies the
following properties:

  * •

$0\leq\phi(\tau,X)\leq 1$ for all $X$ and $\tau\geq 0$,

  * •

if $\phi_{0}(X)$ and $\tilde{\phi}_{0}(X)$ are two initial conditions
satisfying Eq. (92) such that $\phi_{0}(X)\leq\tilde{\phi}_{0}(X)$ for all
$X$, then the corresponding solutions satisfy
$\phi(\tau,X)\leq\tilde{\phi}(\tau,X)$ for all $\tau>0$ and all $X$.

###  B.1 A construction of the solutions $\phi$ and $\phi_{1}$

To show these properties, we rewrite Eq. (91) in terms of
$\varphi:=\phi+\phi_{1}$ and $\psi:=\phi-\phi_{1}$:

$\begin{cases}\partial_{\tau}\varphi+\partial_{X}\varphi=\varphi(\gamma\phi^{2}-1-\kappa)+(1-\gamma)\phi^{3}+\kappa\phi\\\
\partial_{\tau}\psi-\partial_{X}\psi=\psi(\gamma\phi^{2}-1-\kappa)+(1-\gamma)\phi^{3}+\kappa\phi.\end{cases}$
(93)

Going respectively into the $X=z+\tau$ and the $X=z-\tau$ frame, this is
equivalent to

$\begin{cases}\dfrac{{\rm d}}{{\rm
d}\tau}\varphi(\tau,z+\tau)=\Big{[}\varphi(\gamma\phi^{2}-1-\kappa)+(1-\gamma)\phi^{3}+\kappa\phi\Big{]}_{X=z+\tau}\\\\[8.61108pt]
\dfrac{{\rm d}}{{\rm
d}\tau}\psi(\tau,z-\tau)=\Big{[}\psi(\gamma\phi^{2}-1-\kappa)+(1-\gamma)\phi^{3}+\kappa\phi\Big{]}_{X=z-\tau}.\end{cases}$
(94)

Then, this is also equivalent to the integral equations

$\begin{cases}\varphi(\tau,z+\tau)=\displaystyle{\rm e}^{\int_{0}^{\tau}{\rm
d}s\,[\gamma\phi^{2}(s,z+s)-1-\kappa]}\varphi_{0}(z)+\int_{0}^{\tau}{\rm
d}u\,{\rm e}^{\int_{u}^{\tau}{\rm
d}s\,[\gamma\phi^{2}(s,z+s)-1-\kappa]}\big{[}(1-\gamma)\phi(u,z+u)^{3}+\kappa\phi(u,z+u)\big{]}\\\\[8.61108pt]
\psi(\tau,z-\tau)=\displaystyle{\rm e}^{\int_{0}^{\tau}{\rm
d}s\,[\gamma\phi^{2}(s,z-s)-1-\kappa]}\psi_{0}(z)+\int_{0}^{\tau}{\rm
d}u\,{\rm e}^{\int_{u}^{\tau}{\rm
d}s\,[\gamma\phi^{2}(s,z-s)-1-\kappa]}\big{[}(1-\gamma)\phi(u,z-u)^{3}+\kappa\phi(u,z-u)\big{]},\end{cases}$
(95)

where $\varphi_{0}(z)=\phi_{0}(z)+{\phi_{1}}_{0}(z)$ and
$\psi_{0}(z)=\phi_{0}(z)-{\phi_{1}}_{0}(z)$. Note that, from Eq. (92),
$\varphi_{0}(z)=\psi_{0}(z)=\phi_{0}(z)$, but we only need to assume that
$\varphi_{0}$ and $\psi_{0}$ are both in $[0,1]$, i.e. that
$\phi_{0}\pm{\phi_{1}}_{0}\in[0,1]$.

Now, define recursively $\varphi_{n}(\tau,X)$ and $\psi_{n}(\tau,X)$ by

$\displaystyle\begin{cases}\varphi_{1}(\tau,X):=0\\\
\psi_{1}(\tau,X):=0,\end{cases}$ (96)

and

$\begin{cases}\varphi_{n+1}(\tau,z+\tau):=\displaystyle{\rm
e}^{\int_{0}^{\tau}{\rm
d}s\,[\gamma\phi_{n}^{2}(s,z+s)-1-\kappa]}\varphi_{0}(z)+\int_{0}^{\tau}{\rm
d}u\,{\rm e}^{\int_{u}^{\tau}{\rm
d}s\,[\gamma\phi_{n}^{2}(s,z+s)-1-\kappa]}\big{[}(1-\gamma)\phi_{n}(u,z+u)^{3}+\kappa\phi_{n}(u,z+u)\big{]}\\\\[8.61108pt]
\psi_{n+1}(\tau,z-\tau):=\displaystyle{\rm e}^{\int_{0}^{\tau}{\rm
d}s\,[\gamma\phi_{n}^{2}(s,z-s)-1-\kappa]}\psi_{0}(z)+\int_{0}^{\tau}{\rm
d}u\,{\rm e}^{\int_{u}^{\tau}{\rm
d}s\,[\gamma\phi_{n}^{2}(s,z-s)-1-\kappa]}\big{[}(1-\gamma)\phi_{n}(u,z-u)^{3}+\kappa\phi_{n}(u,z-u)\big{]},\end{cases}$
(97)

where
$\phi_{n}(\tau,X):=\frac{1}{2}\big{[}\varphi_{n}(\tau,X)+\psi_{n}(\tau,X)\big{]}.$
We shall show in App. B.3 that, for any choice of $T>0$, the series of
functions $\varphi_{n}$ and $\psi_{n}$ converge uniformly over all $X$ and all
$\tau\in[0,T]$ as $n\to\infty$. The limits $\varphi$ and $\psi$ of these
convergent series must then satisfy Eq. (95) for all $X$ and all $\tau$, which
means that they are the solution to the system (93).

###  B.2 Properties of the solution

Once the convergence of $\varphi_{n}$ and $\psi_{n}$ is proved, we have a
construction of the solution $\\{\varphi,\psi\\}$ to the system (93). This
construction allows us to prove all the properties we claimed:

Light cone property

We observe by recursion that $\varphi_{n}$ and $\psi_{n}$ satisfy the light
cone property. Indeed, it is obvious for $n=1$; assume now that the property
holds for $\varphi_{n}$ and $\psi_{n}$, then $\phi_{n}(u,z+u)$ (and,
similarly, $\phi_{n}(s,z+s)$) in the right hand side of Eq. (97) only depends
on the initial condition in the interval $[z,2u]\in[z,2\tau]$. This implies
that $\varphi_{n+1}(\tau,z+\tau)$ only depends on the initial condition in the
interval $[z,2\tau]$, and this is the light cone property for $\varphi_{n+1}$
(and, similarly, for $\psi_{n+1}$).

Then, going to the limit $n\to\infty$, we obtain that $\varphi$ and $\psi$,
and then $\phi$ and $\phi_{1}$, satisfy the light cone property.

Comparison property

Let us consider two initial conditions $\\{\varphi_{0}(z),\psi_{0}(z)\\}$ and
$\\{\tilde{\varphi}_{0}(z),\tilde{\psi}_{0}(z)\\}$, as well as
$\\{\varphi_{n},\psi_{n},\phi_{n}=\tfrac{1}{2}(\varphi_{n}+\psi_{n})\\}$ and
$\\{\tilde{\varphi}_{n},\tilde{\psi}_{n},\tilde{\phi}_{n}=\tfrac{1}{2}(\tilde{\varphi}_{n}+\tilde{\psi}_{n})\\}$,
the functions respectively obtained from these initial conditions with the
recursion (97). If $\varphi_{0}\leq\tilde{\varphi}_{0}$ and
$\psi_{0}\leq\tilde{\psi}_{0}$, it is clear by recursion from Eq. (97) that,
for all $n$, one has $\varphi_{n}\leq\tilde{\varphi}_{n}$,
$\psi_{n}\leq\tilde{\psi}_{n}$, and then $\phi_{n}\leq\tilde{\phi}_{n}$. Note
that argument makes use of $0\leq\gamma\leq 1$ and $\kappa\geq 0$.

Going to the limit $n\to\infty$, we obtain that $\varphi\leq\tilde{\varphi}$,
$\psi\leq\tilde{\psi}$, and then $\phi\leq\tilde{\phi}$.

Bounds on $\phi$

If $\varphi_{0}\geq 0$ and $\psi_{0}\geq 0$, it is immediately clear that, by
recursion, $\varphi_{n}\geq 0$ and $\psi_{n}\geq 0$, which of course also
implies that $\phi_{n}\geq 0$. With slightly more work, it is also clear that
if $\varphi_{0}\leq 1$ and $\psi_{0}\leq 1$, then by recursion
$\varphi_{n}\leq 1$ and $\psi_{n}\leq 1$ (and $\phi_{n}\leq 1$). Note that
argument makes use of $0\leq\gamma\leq 1$ and $\kappa\geq 0$.

Taking the limit $n\to\infty$, we obtain that $0\leq\varphi\leq 1$ and
$0\leq\psi\leq 1$ (and then $0\leq\phi\leq 1$) if $0\leq\varphi_{0}\leq 1$ and
$0\leq\psi_{0}\leq 1$.

Another way to obtain the same result would be to use the comparison property:
observing that $\\{\varphi=0,\psi=0\\}$ and $\\{\varphi=1,\psi=1\\}$ are two
solutions of (93), then any initial condition “in between” must lead to a
solution remaining between 0 and 1.

###  B.3 Convergence of $\varphi_{n}$ and $\psi_{n}$

Let us now explain why the series $\varphi_{n}$ and $\psi_{n}$ converge as
$n\to\infty$. To simplify the argument, let us assume that the initial
condition satisfies $0\leq\varphi_{0}\leq 1$ and $0\leq\psi_{0}\leq 1$; this
implies from App. B.2 that for all $n$, $0\leq\varphi_{n}\leq 1$,
$0\leq\psi_{n}\leq 1$, and therefore $0\leq\phi_{n}\leq 1$. Let us introduce

$\Delta_{n}(\tau)=\sup_{X\in\mathbb{R}}\big{|}\varphi_{n+1}(\tau,X)-\varphi_{n}(\tau,X)\big{|}+\sup_{X\in\mathbb{R}}\big{|}\psi_{n+1}(\tau,X)-\psi_{n}(\tau,X)\big{|}.$
(98)

Below, we shall show that $\Delta_{n}(\tau)\to 0$ sufficiently quickly as
$n\to\infty$ to ensure it is summable; specifically, we shall show that, for
every $T>0$, one has

$\sum_{n\geq 1}\sup_{\tau\in[0,T]}\Delta_{n}(\tau)<\infty.$ (99)

Then, writing

$\varphi_{n}(\tau,X)=\sum_{p=1}^{n-1}[\varphi_{p+1}(\tau,X)-\varphi_{p}(\tau,X)],$
(100)

(remember that $\varphi_{1}(\tau,X)=0$), we see that the series is convergent,
as the running term is bounded by $\Delta_{p}(\tau)$. This proves that
$\varphi_{n}$ (and, similarly, $\psi_{n}$) converges uniformly over all $X$
and all $\tau\in[0,T]$ to a limiting $\varphi$, as required.

To prove Eq. (99), we use Eq. (97) to bound $\varphi_{n+2}-\varphi_{n+1}$ and
$\psi_{n+2}-\psi_{n+1}$ in terms of $\phi_{n+1}-\phi_{n}$. To do so, we use
extensively $0\leq\gamma\leq 1$, $\kappa\geq 0$ and $0\leq\phi_{n}\leq 1$. In
particular, one has $\gamma\phi_{n}^{2}-1-\kappa\leq 0$.

The mean-value theorem gives the bounds

$\big{|}{\rm e}^{X_{n+1}}-{\rm
e}^{X_{n}}\big{|}\leq|X_{n+1}-X_{n}|\quad\text{if {\hbox{X\leq
0}}},\qquad|Y_{n+1}^{\alpha}-Y_{n}^{\alpha}|\leq\alpha|Y_{n+1}-Y_{n}|\quad\text{if
{\hbox{|Y|\leq 1}}}.$ (101)

Furthermore, as can be seen from writing the left-hand side as $\big{|}({\rm
e}^{X_{n+1}}-{\rm e}^{X_{n}})Y_{n+1}+{\rm e}^{X_{n}}(Y_{n+1}-Y_{n})\big{|}$,

$\big{|}{\rm e}^{X_{n+1}}Y_{n+1}-{\rm
e}^{X_{n}}Y_{n}\big{|}\leq|X_{n+1}-X_{n}|+|Y_{n+1}-Y_{n}|\quad\text{if
{\hbox{X\leq 0}} and {\hbox{|Y|\leq 1}}.}$ (102)

Then, we write $\varphi_{n+2}-\varphi_{n+1}$ using Eq. (97) and bound the
different terms arising in this difference using the relations (101) and
(102). In practice, we need (for $\alpha=0,1,3$)

$\displaystyle\left|{\rm e}^{\int_{u}^{\tau}{\rm
d}s[\gamma\phi_{n+1}(s,z+s)^{2}-1-\kappa]}\phi_{n+1}(u,z+u)^{\alpha}-{\rm
e}^{\int_{u}^{\tau}{\rm
d}s[\gamma\phi_{n}(s,z+s)^{2}-1-\kappa]}\phi_{n}(u,z+u)^{\alpha}\right|$ (103)
$\displaystyle\qquad\qquad\leq\int_{u}^{\tau}{\rm
d}s\,\gamma\big{|}\phi_{n+1}(s,z+s)^{2}-\phi_{n}(s,z+s)^{2}\big{|}+\big{|}\phi_{n+1}(u,z+u)^{\alpha}-\phi_{n}(u,z+u)^{\alpha}\big{|}$
$\displaystyle\qquad\qquad\leq 2\gamma\int_{u}^{\tau}{\rm
d}s\,\big{|}\phi_{n+1}(s,z+s)-\phi_{n}(s,z+s)\big{|}+\alpha\big{|}\phi_{n+1}(u,z+u)-\phi_{n}(u,z+u)\big{|}$
$\displaystyle\qquad\qquad\leq\gamma\int_{u}^{\tau}{\rm
d}s\,\Delta_{n}(s)+\tfrac{\alpha}{2}\Delta_{n}(u).$

In the last line, we used
$|\phi_{n+1}-\phi_{n}|=\frac{1}{2}|\varphi_{n+1}+\psi_{n+1}-\varphi_{n}-\psi_{n}|\leq\frac{1}{2}\Delta_{n}$.
Then, from Eq. (97),

$\displaystyle|\varphi_{n+2}(\tau,z+\tau)-\varphi_{n+1}(\tau,z+\tau)|\leq\varphi_{0}(z)\gamma\int_{0}^{\tau}{\rm
d}s\,\Delta_{n}(s)$ (104)
$\displaystyle\qquad\qquad\qquad\qquad+\int_{0}^{\tau}{\rm
d}u\,\bigg{[}(1-\gamma)\Big{(}\gamma\int_{u}^{\tau}{\rm
d}s\,\Delta_{n}(s)+\tfrac{3}{2}\Delta_{n}(u)\Big{)}+\kappa\Big{(}\gamma\int_{u}^{\tau}{\rm
d}s\,\Delta_{n}(s)+\tfrac{1}{2}\Delta_{n}(u)\Big{)}\bigg{]}.$

We already used that $\varphi_{0}(z)\geq 0$. Recalling that
$\varphi_{0}(z)\leq 1$, we obtain

$|\varphi_{n+2}(\tau,z+\tau)-\varphi_{n+1}(\tau,z+\tau)|\leq\tfrac{a}{2}\int_{0}^{\tau}{\rm
d}s\,\Delta_{n}(s)+\tfrac{b}{2}\int_{0}^{\tau}{\rm d}u\int_{u}^{\tau}{\rm
d}s\,\Delta_{n}(s)=\tfrac{1}{2}\int_{0}^{\tau}{\rm d}s\,(a+bs)\Delta_{n}(s),$
(105)

with $a=2\gamma+3(1-\gamma)+\kappa$ and $b=2(1-\gamma+\kappa)\gamma$. The
exact same bound also holds, obviously, for
$\psi_{n+2}(\tau,z-\tau)-\psi_{n+1}(\tau,z-\tau)|$. Altogether, we finally
obtain

$\Delta_{n+1}(\tau)\leq\int_{0}^{\tau}{\rm d}s\,(a+bs)\Delta_{n}(s).$ (106)

Obviously, from its definition, $\Delta_{0}(\tau)\leq 2$. Furthermore, we have
the relation

$\frac{\big{(}a\tau+\tfrac{1}{2}b\tau^{2})^{n+1}}{(n+1)!}=\int_{0}^{\tau}{\rm
d}s\,(a+bs)\frac{\big{(}as+\tfrac{1}{2}bs^{2})^{n}}{n!}.$ (107)

Indeed, as functions of $\tau$, both sides have the same derivative and the
same value (zero) at $\tau=0$. Then, by recursion, one obtains

$\Delta_{n}(\tau)\leq 2\frac{\big{(}a\tau+\tfrac{1}{2}b\tau^{2})^{n}}{n!}.$
(108)

Clearly, Eq. (99) holds, and the proof is complete.

Note that, in practice, we have shown that the system (91) admits a solution
when the initial conditions satisfy Eq. (92). The same methods can be used to
show that the solution is unique.

##  Appendix C Profile of the Fermi shock wave

In this Appendix, we compute the discontinuity of the Fermi shock wave that
develops at the boundary of the light cone and we derive the implicit equation
which determines the complete shock wave profile. We work in the reference
frame of the right-moving front by using

$\displaystyle\left\\{\begin{array}[]{rl}\displaystyle\phi_{+}(\tau,z):=\phi(\tau,z+R_{0}+\tau)&\xrightarrow[\tau\to+\infty]{}f(z)\\\\[8.61108pt]
\displaystyle\phi_{1+}(\tau,z):=\phi_{1}(\tau,z+R_{0}+\tau)&\xrightarrow[\tau\to+\infty]{}f_{1}(z),\end{array}\right.$
(111)

where we recall that $R_{0}$ is the extension of the initial condition. The
light cone property implies that, at all times, $\phi_{+}(\tau,z)=1$ and
$\phi_{1+}(\tau,z)=0$ for $z>0$. Then the limiting shapes also satisfy
$f(z)=1$ and $f_{1}(z)=0$ for $z>0$. From the coupled PDEs (5), one obtains
the evolution equations for $\phi_{+}$ and $\phi_{1+}$:

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\phi_{+}+\partial_{z}(\phi_{1+}-\phi_{+})&=\phi_{+}(\phi_{+}^{2}-1)\\\\[8.61108pt]
\partial_{\tau}\phi_{1+}-\partial_{z}(\phi_{1+}-\phi_{+})&=\phi_{1+}(\gamma\phi_{+}^{2}-1-\kappa).\end{array}\right.$
(114)

As $\tau\to\infty$, requiring that $\partial_{\tau}\phi_{+}\to 0$ and
$\partial_{\tau}\phi_{1+}\to 0$ leads to the relations

$\displaystyle\partial_{z}(f_{1}-f)=f(f^{2}-1)=f_{1}[1+\kappa-\gamma f^{2}].$
(115)

Hence,

$\displaystyle f_{1}=f\frac{f^{2}-1}{1+\kappa-\gamma f^{2}},$ (116)

and

$\displaystyle\partial_{z}\left(f\frac{f^{2}-1}{1+\kappa-\gamma
f^{2}}-f\right)=f(f^{2}-1).$ (117)

Let us now introduce

$L:=f(0^{-}),\qquad M:=-f_{1}(0^{-}).$ (118)

From Eq. (116), we have

$\displaystyle M=L\frac{1-L^{2}}{1+\kappa-\gamma L^{2}}.$ (119)

Now, integrating the first line of Eq. (114) over $z\in[-\epsilon,\epsilon]$
yields

$\partial_{\tau}\int_{-\epsilon}^{\epsilon}{\rm
d}z\,\phi_{+}(\tau,z)-1-\phi_{1+}(\tau,-\epsilon)+\phi_{+}(\tau,-\epsilon)=\int_{-\epsilon}^{\epsilon}{\rm
d}z\,\phi_{+}(\tau,z)[\phi_{+}(\tau,z)^{2}-1],$ (120)

where we used $\phi_{+}(\tau,\epsilon)=1$ and $\phi_{1+}(\tau,\epsilon)=0$.
Taking $\tau\to\infty$, and then $\epsilon\to 0$ leads to the relation
$M+L=1$. Combining with Eq. (119) gives two solutions for $L$ and $M$
(assuming $L>0$):

$\displaystyle\left\\{\begin{array}[]{rl}L&=1\\\\[6.0pt]
M&=0\end{array}\right.\mbox{ and
}\left\\{\begin{array}[]{rl}L&=\frac{\sqrt{1+4\kappa^{*}(1+\kappa)}-1}{2\kappa^{*}}\\\\[3.0pt]
M&=1-\frac{\sqrt{1+4\kappa^{*}(1+\kappa)}-1}{2\kappa^{*}},\end{array}\right.$
(125)

where we recall that $\kappa^{*}=\gamma+1$. The expression for $L$ in the
second solution is Eq. (9).

For $\kappa>\kappa^{*}$, the second solution in Eq. (125) leads to $L>1$ which
cannot be correct, see Appendix B. This implies that one must have $L=1$ and
$M=0$: for $\kappa>\kappa^{*}$, there is no front at late times at position
$X\simeq\tau$. As is shown in Appendix D, the actual front propagates at a
velocity $v_{\kappa}<1$ and has all the properties of an FKPP front.

For $\kappa<\kappa^{*}$, the second solution in Eq. (125) satisfies $L<1$;
this corresponds to a discontinuous profile where $f(z)$ jumps from $L$ to $1$
across the light cone boundary at $z=0$. Numerically, this is the solution
towards which the system converges for $\kappa<\kappa^{*}$. We can now solve
Eq. (117) for $z<0$ when $\kappa<\kappa^{*}$, using $f(0^{-})=L$ as boundary
condition. This is a first-order ordinary differential equation (ODE) that can
be integrated by means of decomposition into simple fractions. Implicitly,

$\displaystyle\frac{2}{1+\kappa-\gamma f^{2}}-\frac{2}{1+\kappa-\gamma
L^{2}}+\frac{2+\kappa}{1+\kappa}\log\frac{f^{2}}{L^{2}}+\frac{1-\kappa+\gamma}{1+\kappa-\gamma}\log\frac{1-f^{2}}{1-L^{2}}-\frac{3+3\kappa-\gamma}{(1+\kappa)(1+\kappa-\gamma)}\log\frac{1+\kappa-\gamma
f^{2}}{1+\kappa-\gamma L^{2}}=2z.$ (126)

For $\kappa=0$ and $\gamma=1$, the solution explicitly reads $L=M=1/2$,
$f(z<0)=1/\sqrt{1+3{\rm e}^{-z}}$ [17]. As $\kappa\to\kappa^{*}$ from below,
the discontinuity closes continuously as $L\to 1$. However, the front does not
converge to a flat profile. Indeed, linearizing Eq. (126) close to $z=0^{-}$
with $f(z<0)\simeq L+z/\xi$ where $L$ is taken from the second solution in Eq.
(125) and $1/\xi:=f^{\prime}(0^{-})$, we obtain a finite
$\xi(\kappa\to\kappa^{*},\kappa<\kappa^{*})=1/4+\kappa^{*}/2$.

Exactly at the critical disorder $\kappa=\kappa^{*}$, the implicit equation
reads

$\displaystyle\frac{2}{2+\gamma-\gamma
f^{2}}+\frac{3+\gamma}{2+\gamma}\left[\log f^{2}-\log\left(2+\gamma-\gamma
f^{2}\right)\right]=2z+1-\frac{3+\gamma}{2+\gamma}\log 2$ (127)

and a linearization close to $z=0^{-}$ with $f(z<0)\simeq 1+z/\xi$ yields
$\xi(\kappa=\kappa^{*})=1/2+\kappa^{*}$. The factor of two between
$\xi(\kappa\to\kappa^{*},\kappa<\kappa^{*})=\lim_{\kappa\to\kappa^{*-}}\lim_{z\to
0^{-}}1/f^{\prime}(z)=1/4+\kappa^{*}/2$ and $\xi(\kappa=\kappa^{*})=\lim_{z\to
0^{-}}\lim_{\kappa\to\kappa^{*-}}1/f^{\prime}(z)=1/2+\kappa^{*}$ stems from
the existence of a vanishing lengthscale $\propto\kappa-\kappa^{*}$ which
separates two regimes: $f(|z|\ll\kappa-\kappa^{*}\ll 1)$ and
$f(\kappa-\kappa^{*}\ll|z|\ll 1)$. However, further analysis or input on the
model is required to determine which of the solutions in Eq. (125) is selected
at $\kappa=\kappa^{*}$.

##  Appendix D FKPP dynamics

In this Appendix, we discuss how the set of coupled PDEs (5) can exhibit
solutions belonging to the FKPP class [45, 46] for $\kappa>\kappa^{*}$. First,
in App. D.1, we reformulate these equations to better reveal the backbone of
the underlying FKPP physics. Then, in App. D.2, we recall the standard
properties of the solutions to an FKPP equation, and we derive the expression
Eq. (12) of the front velocity. This FKPP prediction is tested at
$\kappa>\kappa^{*}$ in Fig. 1 of the main manuscript. Later, in App. D.3, we
show that the set of coupled PDEs (5) converges to the actual FKPP equation
for $\kappa\gg 1$ (the diffusive limit). Finally in App. D.4, we shall show by
means of extensive numerical simulations that the solutions to the coupled
PDEs (5) have all the expected properties of an FKPP front as soon as
$\kappa>\kappa^{*}$, and not only in the diffusive limit.

###  D.1 Hidden FKPP equation

In this Appendix, we carefully derive Eq. (11) which results from the coupled
PDEs (5). The latter read

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\phi+\bm{\nabla}_{\\!\\!\bm{X}}\cdot\bm{\phi}_{1}\\!\\!\\!\\!&=\phi(\phi^{2}-1)\\\\[6.0pt]
\partial_{\tau}\bm{\phi}_{1}+\bm{\nabla}_{\\!\\!\bm{X}}\phi\\!\\!\\!\\!&=\bm{\phi}_{1}(\gamma\phi^{2}-1)-\kappa\bm{\phi}_{1}.\end{array}\right.$
(130)

Taking the time derivative of the first one and the spatial derivative of the
second one, we have

$\displaystyle\left\\{\begin{array}[]{rl}\partial^{2}_{\tau}\phi+\bm{\nabla}_{\\!\\!\bm{X}}\cdot\partial_{\tau}\bm{\phi}_{1}\\!\\!\\!\\!&=(3\phi^{2}-1)\partial_{\tau}\phi\\\\[6.0pt]
\bm{\nabla}_{\\!\\!\bm{X}}\cdot\partial_{\tau}\bm{\phi}_{1}+\bm{\nabla}_{\\!\\!\bm{X}}^{2}\phi\\!\\!\\!\\!&=\bm{\nabla}_{\\!\\!\bm{X}}\cdot\bm{\phi}_{1}[\gamma\phi^{2}-1-\kappa]+2\gamma\phi\bm{\nabla}_{\\!\\!\bm{X}}\phi\cdot\bm{\phi}_{1}.\end{array}\right.$
(133)

Using the second equation to eliminate the mixed derivative term in the first
one, this yields

$\displaystyle\partial^{2}_{\tau}\phi-\bm{\nabla}_{\\!\\!\bm{X}}^{2}\phi+\bm{\nabla}_{\\!\\!\bm{X}}\cdot\bm{\phi}_{1}[\gamma\phi^{2}-1-\kappa]+2\gamma\phi\bm{\nabla}_{\\!\\!\bm{X}}\phi\cdot\bm{\phi}_{1}=(3\phi^{2}-1)\partial_{\tau}\phi.$
(134)

Now, $\bm{\nabla}_{\\!\\!\bm{X}}\cdot\bm{\phi}_{1}$ can be expressed in terms
of $\phi$ using the first equation in (130), yielding

$\displaystyle\partial^{2}_{\tau}\phi-\bm{\nabla}_{\\!\\!\bm{X}}^{2}\phi+[-\partial_{\tau}\phi+\phi(\phi^{2}-1)][\gamma\phi^{2}-1-\kappa]+2\gamma\phi\bm{\nabla}_{\\!\\!\bm{X}}\phi\cdot\bm{\phi}_{1}=(3\phi^{2}-1)\partial_{\tau}\phi.$
(135)

We switch to the parameter $\kappa^{*}:=1+\gamma$ and reorganize to only have
on the lhs derivative terms with a constant prefactor.

$\left\\{\begin{array}[]{rl}\partial_{\tau}^{2}\phi+(\kappa-\kappa^{*})\partial_{\tau}\phi-\bm{\nabla}_{\\!\\!\bm{X}}^{2}\phi&=\phi(\phi^{2}-1)[(2+\kappa-\kappa^{*})+(1-\kappa^{*})(\phi^{2}-1)]\\\
&\qquad+(2+\kappa^{*})(\phi^{2}-1)\partial_{\tau}\phi-2(\kappa^{*}-1)\phi\bm{\phi}_{1}\cdot\bm{\nabla}_{\\!\\!\bm{X}}\phi\\\\[3.0pt]
\partial_{\tau}\bm{\phi}_{1}+\bm{\nabla}_{\\!\\!\bm{X}}\phi&=\bm{\phi}_{1}[(\kappa^{*}-1)\phi^{2}-1-\kappa].\end{array}\right.$
(136)

The derivative terms on a rhs cancel when $\phi\to 1$ and $\bm{\phi}_{1}\to
0$. Finally, we switch to the field $\delta\phi:=1-\phi$ and obtain the
following set of coupled PDEs

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}^{2}\delta\phi+(\kappa-\kappa^{*})\partial_{\tau}\delta\phi-\bm{\nabla}_{\\!\\!\bm{X}}^{2}\delta\phi&=\delta\phi(1-\delta\phi)(2-\delta\phi)[(2+\kappa-\kappa^{*})-(1-\kappa^{*})\delta\phi(2-\delta\phi)]\\\
&\qquad-(2+\kappa^{*})\delta\phi(2-\delta\phi)\partial_{\tau}\delta\phi-2(\kappa^{*}-1)(1-\delta\phi)\bm{\phi}_{1}\cdot\bm{\nabla}_{\\!\\!\bm{X}}\delta\phi\\\\[3.0pt]
\partial_{\tau}\bm{\phi}_{1}-\bm{\nabla}_{\\!\\!\bm{X}}\delta\phi&=-(2+\kappa-\kappa^{*})\bm{\phi}_{1}-(\kappa^{*}-1)\bm{\phi}_{1}\delta\phi(2-\delta\phi).\end{array}\right.$
(140)

We now re-express the rhs of both equations above by collecting the non-linear
terms under the symbols $\mathrm{NL}$ and $\mathrm{NL}_{1}$. This yields

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}^{2}\delta\phi+(\kappa-\kappa^{*})\partial_{\tau}\delta\phi-\bm{\nabla}_{\\!\\!\bm{X}}^{2}\delta\phi&=2(2+\kappa-\kappa^{*})\delta\phi-\mathrm{NL}\\\\[3.0pt]
\partial_{\tau}\bm{\phi}_{1}-\bm{\nabla}_{\\!\\!\bm{X}}\delta\phi&=-(2+\kappa-\kappa^{*})\bm{\phi}_{1}-\mathrm{NL}_{1}.\end{array}\right.$
(143)

The initial conditions are $\delta\phi(\tau=0,{\bm{X}})=\delta\phi_{0}(X)$,
$\bm{\phi}_{1}(\tau=0,{\bm{X}})=0$, and
$\partial_{\tau}\delta\phi(\tau=0,{\bm{X}})\simeq 2\delta\phi_{0}(X)$. Notice
the dependence on $\bm{\phi}_{1}$ in the first equation which has not been
fully eliminated but is now entering at a non-linear level. Naturally, one may
check that rhs vanishes at the correlated- and uncorrelated-world solutions
given by ($\delta\phi=\bm{\phi}_{1}=0$) and ($\delta\phi=1$,
$\bm{\phi}_{1}=0$), respectively. This concludes the derivation of Eq. (11).

###  D.2 Properties of the FKPP solutions – Short introduction to FKKP-ology

The FKPP equation was introduced independently in 1937 by Fisher [45] and
Kolmogorov-Petrouvsky-Piscounov [46] to describe problems in biology related
to evolution. It has been the subject of many works and studies since. An
extensive review was published in 2003 [26, 47]. We only consider the one-
dimensional $d=1$ case, and focus on $\phi(\tau,X)$ only, regarding $\phi_{1}$
as an intermediate calculation step. We begin by giving a list of defining
characteristics that an evolution equation for $\phi(\tau,X)$ should satisfy
in order to belong to the FKPP class. By evolution equation, here we mean the
actual equation for times $\tau>0$, without the input of the specific initial
condition at time $\tau=0$.

  1. 1.

_The equation admits a stable uniform solution and an unstable uniform
solution._
In our case, $\phi=1$ is an unstable solution (also called the unstable
phase), while $\phi=0$ is a stable solution. (In both phases, we also need to
take $\phi_{1}=0$ but, as already mentioned, we focus on $\phi$ only.)

  2. 2.

_The equation admits traveling wave solutions for all velocities larger or
equal than a critical velocity $v^{*}$._
This means that, for all $v\geq v^{*}$, there exists functions $f_{v}(z)$ with
$f_{v}(\infty)=1$ and $f_{v}(-\infty)=0$ such that $f_{v}(X-v\tau)$ is a
solution. Depending on the initial condition, the front $\phi(\tau,X)$
properly centered might converge to any of these traveling waves.

  3. 3.

_The traveling wave solutions decay exponentially towards the unstable phase
1, with a linear prefactor for the critical traveling wave only._
This means that, there is a spatial decay rate $\mu=\mu(v)>0$ and
$\mu^{*}:=\mu(v^{*})$ such that

$1-f_{v}(z)\propto\exp({-\mu z})\text{ as }z\to\infty\text{ for
}v>v^{*},\quad\text{ and }1-f_{v^{*}}(z)\simeq A\,z\exp({-{\mu^{*}}z})\text{
as }z\to\infty.$ (144)

Nota: in the main manuscript, we wrote $\mu$ and $v$ instead $\mu^{*}$ and
$v^{*}$ as there was no possible ambiguity. The traveling waves are defined up
to translation: if $f_{v^{*}}(z)$ is a traveling wave at velocity $v^{*}$,
then so is $f_{v^{*}}(z-\delta)$. In this Appendix, we choose the particular
critical traveling wave such that
$f_{v^{*}}(z)=1-[A\,z+\epsilon(z)]\exp(-\mu^{*}z)$, where $\epsilon(z)\to 0$
as $z\to\infty$. The quantity $A$ can be measured numerically.

  4. 4.

_When $\delta\phi_{0}(X):=1-\phi(0,X)$ decays faster than $\exp({-\mu^{*}X})$
for large $X$, then the front $\phi(\tau,X)$ propagates asymptotically at
velocity $v^{*}$ and its shape converges to the critical traveling wave._
If we call $m_{\tau}$ the position of the front at time $\tau$, this means
that $m_{\tau}\sim v^{*}\tau$ as $\tau\to\infty$ and that

$\phi(\tau,m_{\tau}+\delta+z)\to f_{v^{*}}(z)\quad\text{as
{\hbox{\tau\to\infty}},}$ (145)

where $\delta$ is a constant offset depending on the initial condition and on
the precise definition of $m_{\tau}$.

  5. 5.

_Still assuming that $\delta\phi_{0}(X):=1-\phi(0,X)$ decays faster than
$\exp({-\mu^{*}X})$ for large $X$, the asymptotic position of the front is_

$m_{\tau}=v^{*}\tau-B\log\tau+a-\frac{E}{\sqrt{\tau}}+\mathcal{O}\Big{(}\frac{\log\tau}{\tau}\Big{)}\quad\text{as
{\hbox{\tau\to\infty}}}.$ (146)

The quantities $B$ and $E$ are known (see below) and independent of the
initial condition. The terms $-B\log\tau$ and $-E/\sqrt{\tau}$ are
respectively known as the Bramson term [27, 48] (see also Refs. [49, 28]) and
the Ebert-van Saarloos term [29] (see also Refs. [30, 31]). The next term, of
order $(\log\tau)/\tau$, has also a known coefficient [32, 33, 34]. The
constant offset $a$ depends on the initial condition and on the precise
definitions of $m_{\tau}$.

  6. 6.

_Still assuming that $\delta\phi_{0}(X):=1-\phi(0,X)$ decays faster than
$\exp({-\mu^{*}X})$ for large $X$, one has for the shape of the front:_

$1-\phi(\tau,m_{\tau}+\delta+z)\simeq
A\,z\exp\left(-\mu^{*}z-C\frac{z^{2}}{\tau}\right)\qquad\text{as
{\hbox{\tau\to\infty}} and {\hbox{z\propto\sqrt{\tau}}}},$ (147)

where $A$ and $\delta$ are the same as in (144) and (145). The quantity $C$
can be computed (see below).

In our problem, the initial conditions $\delta\phi_{0}(X)$ are defined on a
finite support of radius $R_{0}$ and thus do decay faster than
$\exp({-\mu^{*}X})$ for large $X$. This implies that points 4., 5., and 6.
above do apply. In practice, in order to identify $v^{*}$, the first step is
to look for solutions satisfying

$1-\phi(\tau,X)=f_{v}(X-v\tau)\propto\exp\left(-\mu(X-v\tau)\right)\quad\text{for
large {\hbox{z:=X-v\tau}}}.$ (148)

For that purpose, we can consider the equation linearized around $\phi=1$, in
the unstable phase. This leads to a relation between $v$ and $\mu$, which we
can write $v=v(\mu)$. A necessary condition to have an FKPP equation (see
point 2. above) is that $v(\mu)$ reaches a minimum when $\mu=\mu^{*}$ for some
finite value $\mu^{*}$:

$v^{*}:=\min_{\mu}v(\mu)=v(\mu^{*}).$ (149)

For the problem at hand, we inject Eq. (148) into Eq. (143), where we recall
that $\delta\phi:=1-\phi$. Neglecting non-linear terms, this leads to

$(\mu v)^{2}+(\kappa-\kappa^{*})\mu v-\mu^{2}=2(2+\kappa-\kappa^{*}).$ (150)

The solution with $v>0$ is

$v(\mu)=\frac{\sqrt{(4+\kappa-\kappa^{*})^{2}+4\mu^{2}}-(\kappa-\kappa^{*})}{2\mu}.$
(151)

For $\kappa>\kappa^{*}$, one can see that $v(\mu)$ reaches a minimum at a
finite $\mu^{*}$. Conversely, for $\kappa\leq\kappa^{*}$, $v(\mu)$ is always
decreasing and, thus, does not reach a minimum. This indicates that the
problem cannot be FKPP if $\kappa\leq\kappa^{*}$, and that it might be FKPP if
$\kappa>\kappa^{*}$.

From now on, let us assume that $\kappa>\kappa^{*}$. A simple computation from
Eq. (151) yields the selected front velocity and spatial decay rate

$v^{*}=\frac{2\sqrt{4+2(\kappa-\kappa^{*})}}{4+\kappa-\kappa^{*}},\qquad\mu^{*}=\frac{4+\kappa-\kappa^{*}}{\kappa-\kappa^{*}}\sqrt{4+2(\kappa-\kappa^{*})}.$
(152)

This is the solution announced in Eq. (12) of the main manuscript. For an FKPP
front, the behavior of $v(\mu)$ around $\mu=\mu^{*}$ controls many properties
of the propagating front. In particular, the quantities $B$ and $E$ appearing
in (146), and the quantity $C$ appearing in (147) are given by

$B=\frac{3}{2\mu^{*}},\qquad
E=3\sqrt{\frac{2\pi}{\mu^{*5}\,v^{\prime\prime}(\mu^{*})}},\qquad
C=\frac{1}{2\mu^{*}v^{\prime\prime}(\mu^{*})}.$ (153)

See for instance Ref. [48] for $B$, and Ref. [29] for $E$. The above
expression for $C$ can be obtained with the methods developed in Ref. [35].

###  D.3 FKPP equation in the diffusive regime, $\kappa\gg 1$.

Here, we derive the asymptotic FKPP theory when the disorder strength is very
large, $\kappa\gg 1$. We start from the re-formulation of the coupled PDEs in
Eqs. (136) which, we recall, read

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}^{2}\phi+(\kappa-\kappa^{*})\partial_{\tau}\phi-\bm{\nabla}_{\\!\\!\bm{X}}^{2}\phi&=\phi(\phi^{2}-1)[(2+\kappa-\kappa^{*})+(1-\kappa^{*})(\phi^{2}-1)]\\\
&\qquad+(2+\kappa^{*})(\phi^{2}-1)\partial_{\tau}\phi-2(\kappa^{*}-1)\phi\bm{\phi}_{1}\cdot\bm{\nabla}_{\\!\\!\bm{X}}\phi\\\\[3.0pt]
\partial_{\tau}\bm{\phi}_{1}+\bm{\nabla}_{\\!\\!\bm{X}}\phi&=\bm{\phi}_{1}[(\kappa^{*}-1)\phi^{2}-1-\kappa].\end{array}\right.$
(156)

When $\kappa$ is large, the exponential spatial decay rate $\mu^{*}$ is of
order $\sqrt{\kappa}$, see Eq. (152). This indicates that the width of the
front is small, of order $1/\sqrt{\kappa}$. This suggests rescaling space by
working with

${\bm{Y}}:=\sqrt{2\kappa}{\bm{X}}.$ (157)

Keeping only the terms of order $\kappa$, we see that the rescaled equation
for $\phi$ converges as $\kappa\to\infty$ to

$\partial_{\tau}\phi-2\bm{\nabla}_{\\!\\!\bm{Y}}^{2}\phi=\phi(\phi^{2}-1),$
(158)

which is the FKPP equation. Notice that the dependence on $\gamma$, the
distance to the superconducting transition, has dropped. A standard analysis
of the above FKPP equation, following the steps in App. D.2, yields a front
traveling at the velocity $4$ on the ${\bm{Y}}$ scale, and thus with a
velocity of order $2\sqrt{2/\kappa}$ on the ${\bm{X}}$ scale, in agreement
with Eq. (152). The radial profile of the critical traveling wave
corresponding to Eq. (158) must satisfy the ODE

$2\hat{f}^{\prime\prime}+4\hat{f}^{\prime}+\hat{f}(\hat{f}^{2}-1)=0.$ (159)

The prediction is that the late-time shape of the front, $f_{v^{*}}(z)$
defined in Eq. (145), converges after rescaling to

$f_{v^{*}}(Y/\sqrt{2\kappa})\to\hat{f}(Y)\qquad\text{as
{\hbox{\kappa\to\infty}}}.$ (160)

This is numerically confirmed in Fig. 3 of the main manuscript. (Nota: the
solutions to Eq. (159) are defined up to translation.)

###  D.4 Numerical evidence of FKPP dynamics at $\kappa>\kappa^{*}$ – FKPP-
ometry

Here, we present numerical simulations in $d=1$ that demonstrate that
$\phi(\tau,X)$ behaves as an FKPP front as soon as $\kappa>\kappa^{*}$, and
not only when $\kappa\gg\kappa^{*}$ where we showed in App. D.3 that the
equations for $\phi$ were converging to the actual FKPP equation.
Specifically, we show that although the equations for $\phi$ look quite
different from the FKPP equation at finite $\kappa>\kappa^{*}$, the FKPP
properties listed in Eqs. (145), (146) and (147) hold. This is in addition to
the numerical results presented in Fig. 1 of the main manuscript where we
tested the FKPP prediction for the front velocity when $\kappa>\kappa^{*}$.

#### Discretization and integration schemes

Checking that the late-time Eq. (146) holds requires a very accurate
prediction for the velocity of the front. However, as shown below, the
discretization scheme of the equations as well as their integration scheme
slightly impact the value of that velocity. This implies that the analytical
prediction for the velocity must take into account the specific numerical
integration scheme that is used. Another difficulty is that FKPP fronts are
extremely sensitive to the exponential tail ahead of the front, at a distance
on the order of $\sqrt{\tau}$ ahead of the front, where $\phi(\tau,X)$ is very
close to 1. This requires a numerical solver that remains very accurate in
that region. For these reasons, we found it easier to write our own solver,
which we now describe.

Taking the sum and difference of the two coupled PDEs (5), we obtain

$\partial_{\tau}(\phi\pm\phi_{1})\pm\partial_{X}(\phi\pm\phi_{1})=g_{\pm}(\phi,\phi_{1}).$
(161)

with the functions
$g_{\pm}(\phi,\phi_{1}):=\phi(\phi^{2}-1)\pm\phi_{1}[\gamma\phi^{2}-1-\kappa]$.
This is equivalent to the following set of coupled first-order ODEs,

$\frac{{\rm d}}{{\rm
d}\tau}\Big{[}(\phi\pm\phi_{1})(\tau,\pm\tau+z)\Big{]}=g_{\pm}\left(\phi(\tau,\pm\tau+z),\phi_{1}(\tau,\pm\tau+z)\right).$
(162)

We discretize them with a simple Euler scheme, using a common discretization
step $a$ for both space and time. This reads

$(\phi\pm\phi_{1})(\tau+a,\pm(\tau+a)+z)=(\phi\pm\phi_{1})(\tau,\pm\tau+z)+ag_{\pm}\left(\phi(\tau,\pm\tau+z),\phi_{1}(\tau,\pm\tau+z)\right),$
(163)

or, writing $X=\pm(\tau+a)+z$,

$(\phi\pm\phi_{1})(\tau+a,X)=(\phi\pm\phi_{1})(\tau,X\mp
a)+ag_{\pm}\left(\phi(\tau,X\mp a),\phi_{1}(\tau,X\mp a)\right),$ (164)

Finally, solving for $\phi$ and $\phi_{1}$,

$\displaystyle\left\\{\begin{array}[]{rl}\phi(\tau+a,X)=&\frac{\phi(\tau,X-a)+\phi(\tau,X+a)}{2}+\frac{\phi_{1}(\tau,X-a)-\phi_{1}(\tau,X+a)}{2}+a\frac{g_{+}\left(\phi(\tau,X-a),\phi_{1}(\tau,X-a)\right)+g_{-}\left(\phi(\tau,X+a),\phi_{1}(\tau,X+a)\right)}{2}\\\\[6.0pt]
\phi_{1}(\tau+a,X)=&\frac{\phi(\tau,X-a)-\phi(\tau,X+a)}{2}+\frac{\phi_{1}(\tau,X-a)+\phi_{1}(\tau,X+a)}{2}+a\frac{g_{+}\left(\phi(\tau,X-a),\phi_{1}(\tau,X-a)\right)-g_{-}\left(\phi(\tau,X+a),\phi_{1}(\tau,X+a)\right)}{2}.\end{array}\right.$
(167)

These equations allow efficient numerical computation of $\phi(na,pa)$ for any
$n\in\mathbb{N}$ and $p\in\mathbb{Z}$ as a function of the discretized initial
condition $\phi(0,pa)$ for $p\in\mathbb{Z}$. Let us make a few practical
remarks:

  * •

In this scheme, the discretization step $a$ is the same for space and time.
This guarantees a strict light-cone structure: if the initial condition is
zero for $X>R_{0}$, then $\phi(\tau,X)=0$ for $X>R_{0}+\tau$.

  * •

At a given time step, the values on even lattice sites only depend on the
values on odd lattice sites at the previous time step, and vice versa. In
other words, the even and odd lattice sites at a given time step are
independent realizations with slightly different initial conditions. Unless
the initial conditions are extremely fine-tuned, this leads to a front shape
that is not smooth. To avoid this issue, we only display even lattice sites at
even time steps.

  * •

FKPP equations are very sensitive to what happens in the region where the
front is close to the unstable phase, here $\phi=1$. As the default
representation of reals in a computer (we use the “long double” type) is much
better at handling numbers close to 0 rather than numbers close to 1, we
simulate $\delta\phi:=1-\phi$ rather than $\phi$ directly.

  * •

To save computation time and to avoid instabilities due to loss of precision
for numbers too close to 0 or 1, we truncate values at each time step, setting
to 0 any value of $\delta\phi$ smaller than $10^{-2000}$, and setting to 1 any
value larger than $1-10^{-16}$.

The discretized equations (167) lead to a front that propagates at a velocity
slightly different from the one predicted by the continuous equations. In
fact, we can redo the analysis of App. D.2, looking for a traveling wave
solution with an exponential decay $1-\phi(\tau,X)\propto{\rm
e}^{-\mu(X-v\tau)}$. This predicts a relation

$v(\mu)=\frac{\log\cosh(\mu a)}{\mu a}+\frac{1}{\mu
a}\log\\!\left[\\!1+\frac{a}{2}\left(\sqrt{(4+\kappa-\kappa^{*})^{2}+4\frac{\tanh(\mu
a)^{2}}{a^{2}}\Big{(}1-a(\kappa-\kappa^{*})-2a^{2}(2+\kappa-\kappa^{*})\Big{)}}\\!-\\!(\kappa-\kappa^{*})\right)\right],$
(168)

which reduces to Eq. (151) when $a\to 0$.

#### Numerical tests of the FKPP predictions

We computed $\phi(\tau,X)$ for $\gamma=1$ (i.e. $\kappa^{*}=2$) and $\kappa=3$
up to times $\tau=10^{5}$ starting from the $C^{1}$ initial condition

$\displaystyle\delta\phi_{0}(X):=1-\phi(0,X)=0.2(1-X^{2})^{2}\,\Theta(1-|X|),$
(169)

where $\Theta(X)$ is the Heaviside step function. We took a discretization
step $a=10^{-3}$ for which the minimum of Eq. (168) predicts

$v^{*}\approx 0.979\,551\,163,\qquad\mu^{*}\approx 12.1774,\qquad
v^{\prime\prime}(\mu^{*})\approx 2.674\,10^{-4}.$ (170)

These predicted values are to be compared with
$v^{*}=\frac{2}{5}\sqrt{6}\approx 0.979\,796$, $\mu^{*}=5\sqrt{6}\approx
12.2474$ and $v^{\prime\prime}(\mu^{*})=5^{-5}\sqrt{2/3}\approx
2.613\,10^{-4}$ obtained from the continuous equations, see Eq. (152).

In practice, we locate the front using the following definition:

$m_{\tau}:=2a\sum_{n\geq 0}[1-\phi(\tau,2an)]\simeq\int_{0}^{\infty}{\rm
d}X\,[1-\phi(\tau,X)].$ (171)

This definition of the position only takes into account even lattice sites, as
per the discussion above. As can be seen in Fig. 5 below, it turns out that
Eq. (171) is nearly equivalent to defining $m_{\tau}$ as the position where
the front is equal to 0.6.

In Fig. 5, we plot the shape of the front at early times (up to $\tau=10$),
and we observe that, to the eye, it converges quite quickly to the limiting
shape, as per point 4. of App. D.2. The velocity of the front seems already to
be smaller than 1. In Fig. 6, we plot the shape of the front at larger times
(up to $\tau=10^{5}$), appropriately rescaled such as to check point 6. of
App. D.2. Finally, in Fig. 7, we compare $m_{\tau}$ obtained from the program
with the asymptotic prediction (146) given in point 5.

The data in these three figures unambiguously demonstrate that the function
$\phi(\tau,X)$ has all the properties expected for an FKPP front.

Figure 5: $\delta\phi(\tau,X):=1-\phi(\tau,X)$ as a function of $X$, for
$\tau=0,1,\ldots,10$. The initial condition at $\tau=0$ is given in Eq. (169)
and extends to $R_{0}=1$. On each curve, we placed a black dot at the
coordinate corresponding to the front position $m_{\tau}$ defined in Eq.
(171). To test the convergence of the overall front shape at $\tau=10$, we
superimposed in a dashed line the (yellow) curve for $\tau=9$ shifted by
$m_{10}-m_{9}\simeq 0.96945$. To the eye, the shape of the front has already
converged. The quantity $m_{10}-m_{9}$ gives an estimate for the velocity
which is smaller than the predicted velocity $v^{*}$ from Eq. (170); this is
due to the sublinear terms in Eq. (146). The curve for a given $\tau$ stops at
$X=\tau+R_{0}$ (and is implicitly equal to 0 after that point). This is a
direct consequence of the light cone property. The inset shows a zoom of the
bottom-right corner of the main plot. ($\gamma=1$, $\kappa=3$.)  Figure 6:
Test of the prediction (147) for the shape of the front. (left)
$\big{[}1-\phi(\tau,m_{\tau}+\delta+z)]{\rm e}^{\mu^{*}z}$ as a function of
$z=X-m_{\tau}-\delta$, i.e. the fronts are centered around their positions and
the predicted exponential decay has been expunged in order to study
$\delta\phi:=1-\phi$ in the region where it is small. (right) The same
quantity with both axis now rescaled by $\sqrt{\tau}$ converges to
$Ay\exp({-Cy^{2}})$ (dashed line) as $\tau\to\infty$. The values of $\delta$
and $A$ were adjusted manually and the value of $C$ was computed from Eq.
(153). ($\gamma=1$, $\kappa=3$.)  Figure 7: Test of the prediction (146) for
the position of the front as a function of time. In blue, we plot the position
$m_{\tau}$ minus the first two terms of the expansion (146), which are known,
as a function of $1/\sqrt{\tau}$, for $\tau$ up to $10^{5}$. The result should
asymptotically be given by $a-E/\sqrt{\tau}$, where $E$ is known, but not $a$.
The data is compared to the function $a-E/\sqrt{\tau}+(b\log\tau+c)/\tau$
(dashed line), using the known value of $E$ and fitting for $a$, $b$ and $c$.
The dotted line is the simpler asymptotic expansion $a-E/\sqrt{\tau}$. To
illustrate the precision of this test, the position of the front is equal to
$m_{\tau}\simeq 9793.6719$ at $\tau=10^{4}$; the plotted value is
$m_{\tau}-(v^{*}\tau-B\log\tau)\simeq-0.7052$ while
$a-E/\sqrt{\tau}\simeq-0.7058$. Notably, all the digits of the value of
$v^{*}$ given in Eq. (170) are needed to produce this figure. ($\gamma=1$,
$\kappa=3$.)