Traveling discontinuity at the quantum butterfly front

Camille Aron1,2⋆, Eric Brunet1† and Aditi Mitra3‡

1 Laboratoire de Physique de l’École Normale Supérieure, ENS, Université PSL,

CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France

2 Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL),

CH-1015 Lausanne, Switzerland

3 Center for Quantum Phenomena, Department of Physics, New York University,

726 Broadway, New York, NY, 10003, USA

⋆<EMAIL_ADDRESS>†<EMAIL_ADDRESS>‡<EMAIL_ADDRESS>
## Abstract

We formulate a kinetic theory of quantum information scrambling in the context
of a paradigmatic model of interacting electrons in the vicinity of a
superconducting phase transition. We carefully derive a set of coupled partial
differential equations that effectively govern the dynamics of information
spreading in generic dimensions. Their solutions show that scrambling
propagates at the maximal speed set by the Fermi velocity. At early times, we
find exponential growth at a rate set by the inelastic scattering. At late
times, we find that scrambling is governed by shock-wave dynamics with
traveling waves exhibiting a discontinuity at the boundary of the light cone.
Notably, we find perfectly causal dynamics where the solutions do not spill
outside of the light cone.

###### Contents

  1. 1 Introduction
  2. 2 Many-world formalism
    1. 2.1 Motivations and general idea
    2. 2.2 Interacting fermions
    3. 2.3 Green’s functions in the Keldysh basis
    4. 2.4 Inter-world kinetic equation
  3. 3 Superconducting fluctuations
    1. 3.1 Model
    2. 3.2 Inter-world collision integral
  4. 4 Inter-world kinetics in the near-critical regime
    1. 4.1 Partial-wave ansatz
    2. 4.2 Simplified kinetic equation: coupled PDEs
    3. 4.3 Validating the ansatz numerically
      1. 4.3.1 Early times
      2. 4.3.2 Late times and partial-wave truncation
  5. 5 Dynamics of information scrambling
    1. 5.1 Early-time exponential growth
    2. 5.2 Late-time solutions and discontinuous front
      1. 5.2.1 Traveling front
      2. 5.2.2 Saturation inside the light cone
  6. 6 Discussion and conclusion
  7. A Inter-world collision integral
    1. A.1 Critical case
    2. A.2 Away from criticality
  8. B Early-time exponential growth
  9. C Discontinuous front away from criticality ($\gamma\neq 1$)

##  1 Introduction

Quantum information scrambling is the mechanism by which localized information
in an extended closed quantum many-body system with local interactions flows
to non-local degrees of freedom, becoming practically irretrievable. In
practice, these information dynamics can be conveniently probed by means of
out-of-time-ordered correlators (OTOCs) such as squared commutators of
operators inserted at different space-time points, e.g.
$\mathcal{C}(t,{\boldsymbol{x}}):=$
$\langle[\mathcal{O}(t,{\boldsymbol{x}}),\mathcal{O}(0,0)]^{2}\rangle$. The
effective loss of information is characterized by (i) a ballistic spread of
information, often dubbed as the “quantum butterfly effect”, (ii) a growth
regime, reminiscent of the exponential separation between nearby trajectories
in classical chaotic systems, (iii) a purely quantum saturation regime beyond
a scrambling time $t_{*}$.

The ballistic spread of quantum information has been firmly established on the
basis of the Lieb-Robinson bound [1]. The causal light-cone structure, with a
wavefront traveling at a model-dependent butterfly velocity $v_{\rm B}$, was
confirmed in a wide variety of models, from non-interacting 1$d$ systems to
holographic models. Inside the light cone, the existence of a clearly
delineated exponential growth regime is only expected for semiclassical or
large-$N$ models: $\ \ \mathcal{C}(t,{\boldsymbol{x}})\sim$
$\exp\left[\lambda_{\rm L}(t-t_{*}-|{\boldsymbol{x}}|/v_{\rm B})\right]$ where
$\lambda_{\rm L}$ is the Lyapunov exponent. For truly quantum systems, the
rapid growth concentrated near the light cone boundary is understood to be set
by model-dependent microscopic scales, and significant efforts were made to
compute the particular shape of the butterfly front in a variety of models.

Wavefronts described by power laws, sometimes oscillatory, were found in free
and integrable models [2, 3]. Sharp wavefronts were found in interacting spin
chains [4], non-integrable systems with diffusive transport [5, 6], as well as
large-$N$ or semiclassical models [7, 8] and holographic models [9, 10].
Interestingly, random unitary circuits without conserved quantities, i.e. in
the absence of diffusive transport, were found to develop broad fronts
controlled by a diffusively growing length scale [11, 12, 13, 14]. Notably,
several works [15, 16, 17, 18] have pointed towards an effective description
of the front in terms of partial differential equations belonging to the class
of the Fisher Kolmogorov-Petrovsky–Piskunov (FKPP) equation [19, 20] which is
known to exhibit traveling wave solutions (see Ref. [21] and references
therein).

In this manuscript, we address the dynamics and the geometry of the wavefront
for a locally interacting system in the vicinity of a continuous classical
phase transition corresponding to the spontaneous breaking of the symmetry
associated with a conserved quantity. Practically, we consider a paradigmatic
model of interacting electrons where the interactions are due to strong
superconducting fluctuations. The long wavelength fluctuations close to
criticality produce a separation of scales which we leverage to derive
analytic results. How the results get modified on moving away from criticality
is transparent in our derivation and our approach can be systematized to
address other near-critical quantum many-body systems."

We compute the OTOCs by means of an augmented Keldysh formalism, the so-called
many-world formalism, which was originally proposed in Ref. [22] and recently
used to derive kinetic equations for the spreading of quantum information in
fermionic interacting systems including electron-phonon scattering, electrons-
impurity scattering, as well as electron-electron scattering [15]. The
augmented Keldysh formalism has been used in other recent works as well [23,
24, 25, 26]. This formalism can be conveniently harnessed to the standard
field-theoretic concepts, tools, and approximation schemes that have been
developed over the years in condensed matter theory. Here, we treat the
interaction between the electrons and the superconducting fluctuations by
means of the random-phase approximation (RPA) in the particle-particle
channel.

We carefully derive an effective description of the spreading of quantum
information in terms of a set of coupled partial differential equations which
do not belong to the FKPP class. Notably, we find wavefronts that are
discontinuous at the light cone boundary and that do not feature exponentially
small tails ahead of the front. This strictly causal structure is unlike what
is found in evolutions of the FKPP class, and more generally of equations with
diffusive terms.

Figure 1: Butterfly effect in the clean interacting metal defined in Eqs. (8)
and (33): after a small and localized perturbation at time $\tau=0$, the
space-time dynamics of the inter-world distribution function $F_{du}$ defined
in Eq. (24) follows a causal light-cone structure with a front traveling at
the maximal butterfly velocity $v_{\rm B}=v_{\rm F}/\sqrt{d}$ where $v_{\rm
F}$ is the Fermi velocity and $d$ the spatial dimensions. Here, we sketch
$\phi(\tau,{\boldsymbol{X}})$, the component of $F_{du}$ averaged over the
Fermi surface, to be introduced in Eq. (49). At early times, the exponential
growth regime is governed by the inelastic scattering rate. In the late-time
limit, scrambling is governed by shock-wave dynamics with the traveling front
developing a distinct discontinuity followed by an exponential decay on the
scale of the mean free path $\ell$ set by the electronic interaction.

### Summary and main results

The paper is organized as follows. In Sect. 2, we quickly review the many-
world formalism which is used to compute OTOCs and access quantum chaotic
features of many-body systems. It generalizes the standard Keldysh formalism
by studying two replicas of the theory, the so-called worlds, which are only
correlated through their initial conditions. In particular, we introduce the
inter-world distribution function
$F_{\alpha\beta}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}})$, where
$\alpha\neq\beta$ are the world indices, which quantifies the amount of
correlation between the two worlds.

In Sect. 3, we introduce a paradigmatic model describing interacting electrons
close to a superconducting transition. We avoid the technical challenges of
approaching the critical point from within the superconducting phase and work
in the normal phase where no long-range order develops. We derive the
corresponding kinetic equation for the inter-world distribution which is shown
to be of the form

$\displaystyle\partial_{t}F_{\alpha\beta}+\boldsymbol{v}_{\boldsymbol{k}}\cdot\boldsymbol{\nabla}_{\boldsymbol{x}}F_{\alpha\beta}$
$\displaystyle=I_{\alpha\beta}[F_{\alpha\beta},F_{\beta\alpha}],$

with the non-linear collision integral $I_{\alpha\beta}$ given in Eq. (48).

In Sect. 4, we propose and numerically validate an ansatz for
$F_{\alpha\beta}$ leading to a simplified set of non-linear partial
differential equations (PDEs) involving two fields $\phi(t,{\boldsymbol{x}})$
and $\boldsymbol{\phi}_{1}(t,{\boldsymbol{x}})$. These are the first terms of
the partial-wave expansion in momentum space of
$F_{\alpha\beta}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}})$ evaluated on-
shell and at the Fermi surface, i.e. at $\omega=\epsilon_{\boldsymbol{k}}$ and
$k\to k_{\rm F}$. In terms of dimensionless quantities ${\boldsymbol{X}}$ for
space and $\tau$ for time, the PDEs read

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

where the parameter $\gamma$ effectively encodes the distance from the
superconducting phase transition: $\gamma=1$ corresponds to criticality and
$0<\gamma<1$ corresponds to the off-critical regime in the normal phase.

In Sect. 5, starting from a generic localized initial perturbation, we
analytically solve for the dynamics of $\phi$ and $\boldsymbol{\phi}_{1}$, in
any dimension $d$, at criticality as well as away from criticality. The
results are sketched in Fig. 1. The relaxation of the inter-world distribution
is found to strictly occur within a light cone growing from the initial
perturbation at a constant butterfly velocity $v_{\rm B}=v_{\rm F}/\sqrt{d}$
where $v_{\rm F}$ is the Fermi velocity. We work out the early-time dynamics
with an exponential growth of scrambling which is controlled by the inelastic
scattering rate. Importantly, in the late-time regime, we find that scrambling
is governed by shock-wave dynamics, with a traveling wave that develops a
discontinuous radial front of the form $F_{\alpha\beta}\sim
f_{+}\left(\frac{|{\boldsymbol{x}}|-v_{\rm B}t}{\ell/\sqrt{d}}\right)$ which
extends over a length scale set by the mean free path $\ell$ related to the
scattering of the electrons by superconducting fluctuations. Inside the light
cone, $f_{+}$ dies off exponentially away from its boundary
($|{\boldsymbol{x}}|-v_{\rm B}t<0$), and $f_{+}=1$ outside the light cone
($|{\boldsymbol{x}}|-v_{\rm B}t>0$). Notably, we find that $f_{+}$ is
discontinuous precisely at the boundary ($|{\boldsymbol{x}}|-v_{\rm B}t=0$).
We also work out explicitly the exponential falloff governing the approach to
the saturation regime within the bulk of the light cone.

We conclude in Sect. 6 by discussing the relations of our results to previous
works and by giving future directions.

##  2 Many-world formalism

Figure 2: 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$. The location of the operators $\psi$ and
$\psi^{\dagger}$ correspond to the OTOC in Eq. (4).

###  2.1 Motivations and general idea

Let us first motivate the use of the so-called many-world formalism and review
its basic functioning. Dynamical signature of quantum chaos can be found in
OTOCs of the type (we set $\hbar=1$)

$\displaystyle\mathrm{Tr}\,\left[\psi(0)\ {\rm e}^{{\rm
i}Ht}\psi({\boldsymbol{x}}){\rm e}^{-{\rm i}Ht}\ \psi^{\dagger}(0)\ {\rm
e}^{{\rm i}Ht^{\prime}}\psi^{\dagger}({\boldsymbol{x}}^{\prime})\ {\rm
e}^{-{\rm i}Ht^{\prime}}\rho_{0}\right],$ (4)

where $\rho_{0}$ is the initial density matrix at time $t=0$ which is
normalized as $\mbox{Tr}\,\rho_{0}=1$, $H$ is the Hamiltonian generating the
dynamics, and $\psi({\boldsymbol{x}})$ is a local operator at position
${\boldsymbol{x}}$. We have in mind fermionic annihilation/creation operators,
but the discussion can be easily adapted to the bosonic case. Here, the four
operators $\psi(0)$, $\psi({\boldsymbol{x}})$, $\psi^{\dagger}(0)$, and
$\psi^{\dagger}({\boldsymbol{x}}^{\prime})$ are computed in a non-ordered time
sequence. Like many of the diagnostics of quantum chaos, the convoluted time
structure of OTOCs makes them computable objects which do not, however,
directly correspond to physical observables. This is in contrast to standard
retarded correlators which correspond to response functions to physical
perturbations. Consequently, the computation of OTOCs requires modifying the
standard non-equilibrium Green’s function approach to cope with the out-of-
time ordering.

Here, we quickly review the many-world formalism which generalizes the
standard
Schwinger-Keldysh formalism defined on the two-fold Baym-Kadanoff contour to a
formalism on a four-fold contour, suitable to compute four-point OTOCs. This
was originally introduced in Ref. [22] and we refer the reader to Ref. [15]
for a detailed presentation which we simply follow. The OTOC in Eq. (4)
involves two forward and two backward time-evolution operators. Therefore,
following the standard Schwinger-Keldysh construction [27], this yields the
four-fold time-ordered contour $\mathcal{C}$ depicted in Fig. 2. The forward
(backward) branches are labeled by an index $a=+$ ($a=-$). The first two
branches are said to belong to the “up world” and are labeled by the index
$\alpha=u$. The two other branches, posterior on the contour, correspond to
the so-called “down world”, and are labeled by $\alpha=d$. Notably, the up and
down worlds are identical replicas of the same theory, involving the same
Hamiltonian.

In that language, the OTOC in Eq. (4) can be rewritten as

$\displaystyle\langle\mathrm{T}_{\mathcal{C}}\,\psi_{d}^{-}(0,0)\psi_{d}^{+}(t,{\boldsymbol{x}})\psi_{u}^{-\dagger}(0,0)\psi_{u}^{+\dagger}(t^{\prime},{\boldsymbol{x}}^{\prime})\rangle,$
(5)

where $\mathrm{T}_{\mathcal{C}}$ is the time-ordering operator on the contour
$\mathcal{C}$, the operators $\psi$ and $\psi^{\dagger}$ are now written in
the Heisenberg picture, and
$\langle\ldots\rangle:=\mathrm{Tr}\left[\ldots\rho_{0}\right]$.

In our subsequent study of the quantum butterfly effect, we shall assume
equilibrium conditions: the system is initially prepared in thermal
equilibrium at temperature $T$, i.e. with the Gibbs state $\rho_{0}\sim{\rm
e}^{-H/T}$ where we set $k_{\rm B}=1$, and the subsequent evolution is
unitarily generated by the same Hamiltonian $H$ as the one used in the initial
preparation.

It is useful to introduce the following quantities:

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

with the Keldysh indices $a,b=+,-$, the world indices $\alpha,\beta=u,d$, and
where the mixed operator
$\hat{S}_{0}:=\psi_{u}^{-\dagger}(0,0)\psi_{d}^{-}(0,0)$ results from clubbing
of the two operators in Eq. (5) that are both evaluated at time $t=0$ but at
disconnected locations on the time contour $\mathcal{C}$. The OTOC in Eq. (4)
is recovered by taking $\alpha=d$ and $\beta=u$, while $a=\pm$ and $b=\pm$ can
be chosen arbitrarily. If one is to interpret ${\rm
i}G_{\alpha\beta}^{ab}(t,{\boldsymbol{x}};t^{\prime},{\boldsymbol{x}}^{\prime})$
as a two-point function rather than a four-point function, $\hat{S}_{0}$ has
to be understood as a local modification of the initial condition which, a
priori, acts differently on each world: one particle is added to the up world
at position ${\boldsymbol{x}}=0$ while one particle is removed from the down
world. Other choices of operator content are possible for $\hat{S}_{0}$ and we
shall see that the late-time inter-world dynamics depend very little on this
choice. Later, we shall consider local perturbations around the identity:
$\hat{S}_{0}=1+\delta\phi_{0}\,\psi_{u}^{-\dagger}(0,0)\psi_{d}^{-}(0,0)$ with
the infinitesimal parameter $|\delta\phi_{0}|\ll 1$.

Importantly, an inspection of Eq. (6) for $\hat{S}_{0}=1$ shows that the
intra-world Green’s functions (i.e. $\alpha=\beta$) correspond to the standard
(i.e. single-world) Schwinger-Keldysh two-point correlators:

$\displaystyle G^{ab}_{uu}=G^{ab}_{dd}=G^{ab}\,.$ (7)

From now on, given that the intra-world quantities are identical in both
$\alpha=u,d$ worlds, we simply drop the repeated world indices, e.g.
$G_{\alpha\alpha}\to G$, except when this obscures the meaning. Furthermore,
when using the indices $\alpha\beta$, we specifically mean $\alpha\neq\beta$
unless specified otherwise.

###  2.2 Interacting fermions

For concreteness, and to set the stage for the ensuing developments, we work
in the context of interacting fermions on a $d$-dimensional lattice.
Naturally, this can be easily adapted to other quantum systems. We consider
the generic Hamiltonian

$\displaystyle H$ $\displaystyle=H_{0}+H_{\rm int}\,,$ (8) $\displaystyle
H_{0}$ $\displaystyle=\sum_{{\boldsymbol{k}}\in{\rm
BZ}}\sum_{\sigma}\epsilon_{\boldsymbol{k}}c^{\dagger}_{{\boldsymbol{k}}\sigma}c_{{\boldsymbol{k}}\sigma}\,,$
(9)

where $H_{0}$ is the non-interacting part and the interaction in $H_{\rm int}$
depends on the specific problem at hand, see Eq. (33). The fermionic operator
$c_{{\boldsymbol{k}}\sigma}^{\dagger}$ creates a electrons with spin
$\sigma=\uparrow$ or $\downarrow$ ($\bar{\sigma}=\downarrow$, $\uparrow$) and
momentum ${\boldsymbol{k}}$ in the Brillouin zone (BZ).
$\epsilon_{\boldsymbol{k}}$ is the dispersion relation. The generalization to
multi-band cases is straightforward. For simplicity, we measure electronic
energies relative to the chemical potential, but a finite chemical potential
$\mu$ can be included via the substitution
$\epsilon_{\boldsymbol{k}}\mapsto\epsilon_{\boldsymbol{k}}-\mu$. For
simplicity, we shall assume that the Fermi surface is spherical, i.e.
$\epsilon_{\boldsymbol{k}}=0$ when $k\to k_{\rm F}$. Whenever this does not
harm the understanding, we shall simply drop the spin indices.

###  2.3 Green’s functions in the Keldysh basis

The 16 Green’s functions $G_{\alpha\beta}^{ab}$ are not independent of each
other and one may considerably reduce the redundancies of the formalism. On
the one hand, the causal structure of the contour $\mathcal{C}$ is such that
the inter-world Green’s functions (i.e. $\alpha\neq\beta$) do not depend on
the $\pm$ basis:

$\displaystyle G^{++}_{\alpha\beta}$
$\displaystyle=G^{--}_{\alpha\beta}=G^{+-}_{\alpha\beta}=G^{-+}_{\alpha\beta}\
\mbox{ for }\alpha\neq\beta\,.$ (10)

On the other hand, the intra-world Green’s functions (i.e. $\alpha=\beta$) are
related via [27]

$\displaystyle
G_{\alpha\beta}^{+-}+G_{\alpha\beta}^{-+}=G_{\alpha\beta}^{++}+G_{\alpha\beta}^{--}\
\mbox{ for }\alpha=\beta\,.$ (11)

It is therefore customary to perform a rotation from the $a=\pm$ basis to the
so-called Keldysh basis and work with retarded, advanced, and Keldysh Green’s
functions:

$\displaystyle\left.\begin{array}[]{c}G_{\alpha\beta}^{R}=G_{\alpha\beta}^{++}-G_{\alpha\beta}^{+-}\,,\qquad
G_{\alpha\beta}^{A}=G_{\alpha\beta}^{++}-G_{\alpha\beta}^{-+}\,,\\\\[6.0pt]
G_{\alpha\beta}^{K}=G_{\alpha\beta}^{++}+G_{\alpha\beta}^{--}\,,\end{array}\right.$
(14)

where $G^{A}$ is simply the Hermitian conjugate of $G^{R}$.

The intra-world Green’s functions are the standard Schwinger-Keldysh
correlators, solutions of the Schwinger-Dyson equations which, in thermal
equilibrium and in Fourier space, read

$\displaystyle\left.\begin{array}[]{rl}G^{R}(\omega,{\boldsymbol{k}})\\!\\!\\!\\!&=\left[\omega-\epsilon_{\boldsymbol{k}}-\Sigma^{R}(\omega,{\boldsymbol{k}})\right]^{-1},\\\\[6.0pt]
G^{K}(\omega,{\boldsymbol{k}})\\!\\!\\!\\!&=2{\rm
i}\,F(\omega)\,\mathrm{Im}\,G^{R}(\omega,{\boldsymbol{k}})\,.\end{array}\right.$
(17)

$\Sigma^{R}$ is the retarded component of the self-energy. It is due to the
interaction in $H_{\rm int}$ and can be computed diagrammatically within the
standard Schwinger-Keldysh formalism. The last equality is the expression of
the fermionic fluctuation-dissipation theorem with $F(\omega):=$
$\tanh\left(\omega/2T\right)$.

Concerning inter-world Green’s functions, the relation in Eq. (10) together
with Eq. (14) immediately implies

$\displaystyle G_{ud}^{R}$ $\displaystyle=G_{du}^{R}=0\,.$ (18)

This expresses the fact that while intra-world physics contributes to inter-
world quantities, the opposite, namely that inter-world physics contributes to
intra-world quantities, is strictly forbidden.

To summarize, given that we restrict ourselves to equilibrium physics, we are
to deal with only three independent Greens’ functions: the standard (intra-
world) retarded Green’s function $G^{R}$ which is uniquely specified by the
Hamiltonian $H$ and the temperature $T$, and the two inter-world Keldysh
Green’s functions $G_{ud}^{K}$ and $G_{du}^{K}$ which also depend on the
choice of the operator $\hat{S}_{0}$. For the choice $\hat{S}_{0}=1$, the two
worlds have the same thermal initial conditions, and one may check that
$G_{\alpha\beta}^{K}$ are space- and time-translational invariant and

$\displaystyle\left.\begin{array}[]{rl}G_{ud}^{K}(\omega,{\boldsymbol{k}})\\!\\!\\!\\!&=2{\rm
i}[-1+F(\omega)]\,\mathrm{Im}\,G^{R}(\omega,{\boldsymbol{k}})\,,\\\\[6.0pt]
G_{du}^{K}(\omega,{\boldsymbol{k}})\\!\\!\\!\\!&=2{\rm
i}[+1+F(\omega)]\,\mathrm{Im}\,G^{R}(\omega,{\boldsymbol{k}})\,.\end{array}\right.$
(21)

However, for a generic choice of $\hat{S}_{0}$, $G^{K}_{\alpha\beta}$ is not
guaranteed to be space- and time-translational invariant and it can be
determined via the Schwinger-Dyson equation reading

$\displaystyle\textstyle
G^{K}_{\alpha\beta}(t,{\boldsymbol{x}};t^{\prime},{\boldsymbol{x}}^{\prime})\\!=\\!\\!\int\\!\\!{\rm
d}{\boldsymbol{x}}_{1}\\!\\!\int\\!\\!{\rm
d}{\boldsymbol{x}}_{2}\\!\\!\int\\!\\!{\rm d}t_{1}\\!\\!\int\\!\\!{\rm
d}t_{2}G^{R}(t-t_{1},{\boldsymbol{x}}-{\boldsymbol{x}}_{1})\Sigma^{K}_{\alpha\beta}(t_{1},{\boldsymbol{x}}_{1};t_{2},{\boldsymbol{x}}_{2},)G^{A}(t_{2}-t^{\prime},{\boldsymbol{x}}_{2}-{\boldsymbol{x}}^{\prime}),$
(22)

where $\Sigma^{K}_{\alpha\beta}$ is the Keldysh component of the inter-world
self-energy which can be computed diagrammatically in the many-world
formalism.

###  2.4 Inter-world kinetic equation

It is useful to work in the Wigner representation

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

and parameterize $G^{K}_{\alpha\beta}$ in terms of the real function
$F_{\alpha\beta}$:

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

where we introduced the Moyal product $\star:=\exp\left[\frac{{\rm
i}}{2}(\overleftarrow{\partial_{\omega}}\overrightarrow{\partial_{t}}\\!-\\!\overleftarrow{\boldsymbol{\nabla}_{\boldsymbol{k}}}\cdot\overrightarrow{\boldsymbol{\nabla}_{\boldsymbol{x}}}\\!-\\!\overleftarrow{\partial_{t}}\overrightarrow{\partial_{\omega}}\\!+\\!\overleftarrow{\boldsymbol{\nabla}_{\boldsymbol{x}}}\overrightarrow{\boldsymbol{\nabla}_{\boldsymbol{k}}})\right]$
where the left (right) arrow designates a derivative operator acting on the
left (right) of the star symbol. In analogy to the standard (intra-world)
electronic distribution function $F$, $F_{\alpha\beta}$ is dubbed the inter-
world distribution function. One has
$F_{ud}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}})\in[-2,0]$ and
$F_{du}\in[0,2]$. Assuming that the variations of $F_{\alpha\beta}$ occur on
scales much larger than the microscopic scales involved in $G^{R}$, we may
work in the so-called quasi-classical or gradient approximation which consists
in truncating the derivative expansion to its leading terms.

The interworld distribution $F_{\alpha\beta}$ gives access to the information
scrambling as it is related to the original OTOC in Eq. (5). Let us briefly
outline this connection. At late times, we expect the decoupling

$\displaystyle\langle\mathrm{T}_{\mathcal{C}}\,\psi_{d}^{-}(0,0)\psi_{d}^{+}(t,{\boldsymbol{x}})\psi_{u}^{-\dagger}(0,0)\psi_{u}^{+\dagger}(t^{\prime},{\boldsymbol{x}}^{\prime})\rangle$
$\displaystyle\qquad\simeq(1-n_{0})\,{\rm
i}G_{du}^{++}(t,{\boldsymbol{x}};t^{\prime},{\boldsymbol{x}}^{\prime})\propto
G^{K}_{du}(t,{\boldsymbol{x}};t^{\prime},{\boldsymbol{x}}^{\prime})\propto
F_{du}(t,{\boldsymbol{x}};t^{\prime},{\boldsymbol{x}}^{\prime}),$ (25)

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

Massaging Eq. (24) by acting on both sides with the inverse of $G^{R}$ and
$G^{A}$, and using the Dyson equations (17) and (22), one derives the kinetic
equation on $F_{\alpha\beta}$ which is analogous to the standard (intra-world)
kinetic equation:

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

where the lhs corresponds to the non-interacting physics set by $H_{0}$, with
the velocity
$\boldsymbol{v}_{\boldsymbol{k}}:=\boldsymbol{\nabla}_{\boldsymbol{k}}\epsilon_{\boldsymbol{k}}$,
and the right-hand side (rhs) is the so-called collision integral which stems
from $H_{\rm int}$ and reads

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

Let us recall that Keldysh components and collision integrals depend on space
and time, namely ${\boldsymbol{x}}$ and $t$, through the inter-world
distribution functions. However, here and from now on, we simplify the
notation by dropping the explicit dependence on these objects.

As discussed above, the intra-world quantity $\Sigma^{R}$ cannot depend on
$F_{\alpha\beta}$, however $\Sigma^{K}_{\alpha\beta}$ is expected to be a non-
linear functional of $F_{\alpha\beta}$. The inter-world kinetic equation in
Eq. (26) is therefore a non-linear partial integrodifferential equation. It
has a trivial steady-state solution

$\displaystyle F_{ud}^{\rm uncorr}(\omega,{\boldsymbol{k}})=F_{du}^{\rm
uncorr}(\omega,{\boldsymbol{k}})=0,$ (28)

which reflects a total loss of coherence between the two replicated worlds,
and which is dubbed the “uncorrelated-world” solution. Additionally, one can
easily check that the case of $\hat{S}_{0}=1$ in Eq. (21), where both worlds
evolve coherently, corresponds to another steady state characterized by

$\displaystyle\left.\begin{array}[]{l}F_{ud}^{\rm
corr}(\omega,{\boldsymbol{k}})=-1+\tanh\left(\omega/2T\right)\,,\\\\[6.0pt]
F_{du}^{\rm
corr}(\omega,{\boldsymbol{k}})=+1+\tanh\left(\omega/2T\right)\,,\end{array}\right.$
(31)

and which is dubbed the “correlated-world” solution. As we shall see
explicitly later, the correlated-world solution at $\hat{S}_{0}=1$ is expected
to be unstable against small perturbations of $\hat{S}_{0}$ and the only
stable steady-state is the uncorrelated-world solution.

The aim of this manuscript is to derive and analyze the dynamics of the inter-
world distribution function of a near-critical system of interacting electrons
when the initial condition is of the form

$\displaystyle F_{\alpha\beta}(\omega,{\boldsymbol{k}};t=0,{\boldsymbol{x}})$
$\displaystyle=[1-\delta\phi_{0}({\boldsymbol{x}})]F_{\alpha\beta}^{\rm
corr}(\omega,{\boldsymbol{k}}),$ (32)

where $0\leq\delta\phi_{0}({\boldsymbol{x}})\ll 1$ is an initial perturbation
to the correlated-world solution localized on a compact support around
${\boldsymbol{x}}=0$. For simplicity, we consider a perturbation that is the
same for both $F_{ud}$ and $F_{du}$.

##  3 Superconducting fluctuations

###  3.1 Model

Figure 3: Superconducting self-energy $\Sigma$ and polarization bubble $\Pi$
of the model in Eq. (33) treated in the RPA scheme in the particle-particle
channel. The expressions in the Keldysh formalism are given in Eqs. (43) and
(39), respectively.

Concretely, we consider the standard Hubbard-like electron-electron
interaction restricted to the particle-particle (Cooper) channel. The
interacting piece of the Hamiltonian in Eq. (8) reads

$\displaystyle H_{\rm int}$
$\displaystyle=U\sum_{{\boldsymbol{k}}{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\sum_{\sigma}c^{\dagger}_{{\boldsymbol{k}}\sigma}c^{\dagger}_{-{\boldsymbol{k}}+{\boldsymbol{k}}^{\prime}\bar{\sigma}}c_{{\boldsymbol{k}}^{\prime\prime}+{\boldsymbol{k}}^{\prime}\bar{\sigma}}c_{-{\boldsymbol{k}}^{\prime\prime}\sigma}\,,$
(33)

where $U<0$ is an attractive interaction facilitating superconductivity. In
dimensions $d\geq 2$, this model exhibits a finite-temperature phase
transition towards a superconducting phase associated with the spontaneous
breaking of the $U(1)$ symmetry. Here, we consider the near-critical regime,
above the critical temperature where the $U(1)$ symmetry is not broken but the
superconducting fluctuations are sizable.

The (intra-world) physics of this model is well understood, and we rely on
standard and well-tested methods which we extend to the many-world formalism.
In practice, we decouple the Hubbard interaction in the Cooper channel and
obtain a theory of fermions coupled to bosonic fluctuations. The Cooperon
Green’s functions within RPA in the particle-particle channel, are given by
[28, 29, 30, 31]

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

where $\Pi^{R}$ is the retarded Cooper bubble and the last equality is the
bosonic fluctuation-dissipation theorem with
$P(\omega):=\coth\left(\omega/2T\right)$. The retarded Cooper bubble and the
electronic self-energies within RPA, i.e. limiting ourselves to the one-loop
diagrams depicted in Fig. 3 are:

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

and

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

The real parts of the retarded components can be recovered via the
Kramers–Kronig relation. A self-consistent treatment of Green’s functions,
self-energies, and bubbles ensures that the RPA scheme is a conserving
approximation. It is known to be exact in the large $N$-limit, where $N$ is
the number of electronic orbitals [28, 29, 30, 31].

Close to criticality, the Cooperon propagator reads [32, 28, 29, 30, 31]

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

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
point. 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$.

###  3.2 Inter-world collision integral

Let us now discuss the expressions of inter-world quantities which are
necessary to compute the inter-world collision integral in Eq. (27). As we
already noted, the inter-world retarded components of the Green’s functions
(fermionic and Cooperon), the bubbles and the self-energies simply vanish as a
consequence of the fact that intra-world physics can be expressed
independently of inter-world quantities:
$G^{R}_{\alpha\beta}=D^{R}_{\alpha\beta}=\Pi^{R}_{\alpha\beta}=\Sigma^{R}_{\alpha\beta}=0$
for $\alpha\neq\beta$. Within the RPA scheme, the inter-world Keldysh
components read

$\displaystyle D_{\alpha\beta}^{K}(\omega,{\boldsymbol{k}})$
$\displaystyle=\left|D^{R}(\omega,{\boldsymbol{k}})\right|^{2}\Pi_{\alpha\beta}^{K}(\omega,{\boldsymbol{k}})\,,$
(45) $\displaystyle\Pi_{\alpha\beta}^{K}(\omega,{\boldsymbol{k}})$
$\displaystyle=\frac{{\rm i}}{2}\sum_{{\boldsymbol{k}}^{\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}G_{\alpha\beta}^{K}(\omega^{\prime},{\boldsymbol{k}}^{\prime})G_{\alpha\beta}^{K}(\omega-\omega^{\prime},{\boldsymbol{k}}-{\boldsymbol{k}}^{\prime})\,,$
(46) $\displaystyle\Sigma_{\alpha\beta}^{K}(\omega,{\boldsymbol{k}})$
$\displaystyle=-\frac{{\rm
i}}{2}\sum_{{\boldsymbol{k}}^{\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}D_{\alpha\beta}^{K}(\omega^{\prime},{\boldsymbol{k}}^{\prime})G_{\beta\alpha}^{K}(\omega^{\prime}-\omega,{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\,.$
(47)

Altogether, this yields the inter-world kinetic equation (26) with the
collision integral

$\displaystyle I_{\alpha\beta}(\omega,{\boldsymbol{k}})$
$\displaystyle=2\sum_{{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}$
$\displaystyle\qquad\times\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega,{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})$
$\displaystyle\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,{\boldsymbol{k}})$
$\displaystyle\qquad\qquad+F_{\alpha\beta}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})F_{\alpha\beta}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})F_{\beta\alpha}(\omega^{\prime}-\omega,{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\bigg{\\}}\,.$
(48)

We recall that we omitted the local space and time dependence of
$I_{\alpha\beta}$ and $F_{\alpha\beta}$ to shorten the notations. As a sanity
check, one may verify that the collision integral, and more precisely the term
inside the curly brackets, vanishes at both the uncorrelated-world solution in
Eq. (28) and the correlated-world solution in Eq. (31).

It is worthwhile noting that, contrary to standard intra-world collision
integrals, there are no underlying conservation laws associated with the
inter-world electronic distribution $F_{\alpha\beta}$(such as number, energy,
momentum conservation) that guarantee sum rules such as
$\int{\rm d}\omega\sum_{{\boldsymbol{k}}}I(\omega,{\boldsymbol{k}})=0$. In
turn, this lack of underlying conserved quantities has important consequences
on the relaxation dynamics of the inter-world distribution function. Indeed,
the presence of conservation laws implies a separation of timescales and is
typically synonymous with diffusive dynamics for perturbations that vary
slowly enough.

##  4 Inter-world kinetics in the near-critical regime

In this Section, we propose and numerically validate an ansatz to the inter-
world distribution $F_{\alpha\beta}$ which allows us to derive a much simpler
version of the kinetic equation (26) with its collision integral in Eq. (48).
This is done in the vicinity of the superconducting transition where a clear
separation of energy scales can be made. The resulting effective description
of the information scrambling dynamics consists of a set of coupled PDEs that
we solve analytically in Sect. 5.

###  4.1 Partial-wave ansatz

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

Furthermore, given that relaxation is dominated by the electronic states
around the Fermi level, at energy scales (e.g. temperature) that are much
smaller than the Fermi energy, one may focus on the distribution function
close to the Fermi surface by subsequently setting $k\to k_{\rm F}$.

We now propose to simplify drastically the partial integrodifferential kinetic
equation in Eq. (26) with the following ansatz:

$\displaystyle\tilde{F}^{\rm
ansatz}_{\alpha\beta}(t,{\boldsymbol{x}};{\boldsymbol{k}})$
$\displaystyle=\left[\phi(t,{\boldsymbol{x}})+{\boldsymbol{u}}_{\boldsymbol{k}}\cdot\boldsymbol{\phi}_{1}(t,{\boldsymbol{x}})\right]\,\tilde{F}_{\alpha\beta}^{\rm
corr}({\boldsymbol{k}}),$ (49)

where the unit vector ${\boldsymbol{u}}_{\boldsymbol{k}}:={\boldsymbol{k}}/k$.
The two fields $\phi$ and $\boldsymbol{\phi}_{1}$ can be understood as the
first terms of a partial-wave expansion [27] of $\tilde{F}_{\alpha\beta}$,
accounting for its isotropic and anisotropic contributions in momentum space:

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

where $S_{d-1}:=\int{\rm d}\Omega_{\boldsymbol{k}}$ is the surface area of the
$d-1$-sphere with unit radius and ${\rm d}\Omega_{\boldsymbol{k}}$ is the
elementary solid angle in the direction of ${\boldsymbol{k}}$. We did not
include higher-order terms in the ansatz, e.g. of the form
$u^{i}_{\boldsymbol{k}}u^{j}_{\boldsymbol{k}}\phi_{2}^{ij}$.

Let us give the rationale behind the ansatz proposed in Eq. (49). Firstly, let
us note that standard (intra-world) approaches consist in perturbing the
distribution function around its equilibrium value, $F=F^{\rm eq}+\delta F$,
and linearizing the collision integral accordingly. This approach relies on
the fact that the equilibrium distribution $F^{\rm
eq}(\omega,{\boldsymbol{k}})=\tanh(\omega/2T)$ is a stable steady state of the
(intra-world) kinetic equation, guaranteed by the H-theorem. The inter-world
case is much different as the initial condition set by $F_{\alpha\beta}^{\rm
corr}$ is unstable and one cannot propose a perturbative ansatz. This explains
why $\tilde{F}_{\alpha\beta}^{\rm corr}$ appears multiplicatively in Eq. (49)
and why the collision integral cannot be, a priori, linearized. In that
regard, it is similar to the ansatz used in Ref. [15].

Secondly, close to the Fermi surface which is assumed to be spherical, the
solutions $\tilde{F}^{\rm corr}_{\alpha\beta}({\boldsymbol{k}})$ and
$\tilde{F}^{\rm uncorr}_{\alpha\beta}({\boldsymbol{k}})$ do not depend on the
direction of the momentum ${\boldsymbol{k}}$, but only on its norm $k\approx
k_{\rm F}$. This means that we aim at describing the dynamics of
$\tilde{F}_{\alpha\beta}$ from a momentum-space isotropic and real-space
homogeneous (unstable) solution

$\displaystyle\tilde{F}^{\rm
corr}({\boldsymbol{k}};t=0,{\boldsymbol{x}})\longleftrightarrow\left\\{\begin{array}[]{rl}\phi(t=0,{\boldsymbol{x}})\\!\\!\\!\\!&=1,\\\
\boldsymbol{\phi}_{1}(t=0,{\boldsymbol{x}})\\!\\!\\!\\!&=0,\end{array}\right.$
(55)

to another momentum-space isotropic and real-space homogeneous (stable)
solution

$\displaystyle\tilde{F}^{\rm
uncorr}({\boldsymbol{k}};t\to\infty,{\boldsymbol{x}})\longleftrightarrow\left\\{\begin{array}[]{rl}\phi(t\to\infty,{\boldsymbol{x}})\\!\\!\\!\\!&=0,\\\
\boldsymbol{\phi}_{1}(t\to\infty,{\boldsymbol{x}})\\!\\!\\!\\!&=0\,.\end{array}\right.$
(58)

However, as will become clear below, these dynamics can only proceed by
allowing anisotropy in momentum space to develop in the transient regime
towards the stable steady state. This explains why we included the anisotropic
term in $\boldsymbol{\phi}_{1}$ which can be seen as the minimal ingredient to
allow for spatial relaxation.

Thirdly, let us note that $\tilde{F}^{\rm
ansatz}_{\alpha\beta}(t,{\boldsymbol{x}};{\boldsymbol{k}})$ depends on $k$
only through $\tilde{F}_{\alpha\beta}^{\rm corr}({\boldsymbol{k}})$. If this
is trivially true at the correlated- and uncorrelated-world solutions, we
shall see later that this is also compatible with the dynamics which does not
generate extra dependence on $k$.

Finally, let us note that the fields $\phi$ and $\boldsymbol{\phi}_{1}$ are
common to both $\tilde{F}_{ud}$ and $\tilde{F}_{du}$. This stems from our
choice of initial perturbation in Eq. (32).

###  4.2 Simplified kinetic equation: coupled PDEs

We proceed by injecting the ansatz (49) in the inter-world collision integral
in Eq. (48) and consistently truncating its partial-wave expansion to the two
lowest orders. Because this brings further simplifications, we work in the
near-critical regime of the symmetric (normal) phase where the Cooperon
becomes soft. In Eq. (48), this means that the term
$|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}$ diverges as
$\omega^{\prime}\approx 0$ and $k^{\prime}\approx 0$. The details of the
computation are given in Appendix A. We obtain

$\displaystyle\tilde{I}_{\alpha\beta}({\boldsymbol{k}})=2(1-\phi^{2})\left[\phi+(\boldsymbol{\phi}_{1}\cdot{\boldsymbol{u}}_{\boldsymbol{k}})\right]\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}})\,{\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}})\,.$ (59)

Injecting the ansatz in the kinetic equation (26), dividing by
$\tilde{F}_{\alpha\beta}^{\rm corr}({\boldsymbol{k}})$, we get

$\displaystyle\partial_{t}\phi+v_{\boldsymbol{k}}\left({\boldsymbol{u}}_{\boldsymbol{k}}\cdot\boldsymbol{\nabla}_{\boldsymbol{x}}\right)\phi+{\boldsymbol{u}}_{{\boldsymbol{k}}}\cdot\partial_{t}\boldsymbol{\phi}_{1}+v_{\boldsymbol{k}}\left({\boldsymbol{u}}_{\boldsymbol{k}}\cdot\boldsymbol{\nabla}_{\boldsymbol{x}}\right)({\boldsymbol{u}}_{\boldsymbol{k}}\cdot\boldsymbol{\phi}_{1})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad=2(1-\phi^{2})\left[\phi+(\boldsymbol{\phi}_{1}\cdot{\boldsymbol{u}}_{\boldsymbol{k}})\right]\,{\rm
Im}\,\tilde{\Sigma}^{R}({\boldsymbol{k}})\,.$ (60)

We now project on the momentum-space isotropic and first partial-wave
contributions by acting with $\frac{1}{S_{d-1}}\int{\rm
d}\Omega_{\boldsymbol{k}}$ and $\frac{d}{S_{d-1}}\int{\rm
d}\Omega_{\boldsymbol{k}}{\boldsymbol{u}}_{\boldsymbol{k}}$ on both sides of
the above equation. We use $\int{\rm
d}\Omega_{\boldsymbol{k}}u_{\boldsymbol{k}}^{i}u_{\boldsymbol{k}}^{j}=$
$\delta_{ij}\,S_{d-1}/d$. At the Fermi surface, i.e. eventually setting $k\to
k_{\rm F}$, we obtain the following set of coupled partial differential
equations (PDEs)

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

where $v_{{\rm F}}$ is the Fermi velocity and we defined the timescale which
sets the fermionic lifetime as (recalling that the self-energy only depends on
the norm of ${\boldsymbol{k}}$)

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

Note that the temperature dependence of the original problem enters through
$\tau_{\rm F}$ and $\gamma$.

$\gamma$ is a dimensionless parameter that generalizes the computation
performed at criticality, for which $\gamma=1$, to near critical regimes for
which $\gamma<1$ (see Appendix A.2). After an appropriate rescaling of space
and time,

$\displaystyle\tau:=t/\tau_{\rm F}\ \mbox{ and }\
{\boldsymbol{X}}:={\boldsymbol{x}}\sqrt{d}/(v_{\rm F}\tau_{\rm F})\,,$ (65)

together with $\boldsymbol{\phi}_{1}\mapsto\sqrt{d}\boldsymbol{\phi}_{1}$, the
coupled PDEs now only involve dimensionless quantities

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

Importantly, the generic $d$-dimensional case can be reduced to an effective
one-dimensional case. Indeed, assuming spherically-symmetric initial
conditions, we may work with the radial coordinate $r$: $\phi(\tau,r)$ and
$\boldsymbol{\phi}_{1}=\phi_{1}(\tau,r){\boldsymbol{u}}_{r}$ at all times. The
coupled PDEs now read

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\phi+\partial_{r}\phi_{1}\\!\\!\\!\\!&=-\frac{d-1}{r}\phi_{1}+\phi(\phi^{2}-1),\\\\[6.0pt]
\partial_{\tau}\phi_{1}+\partial_{r}\phi\\!\\!\\!\\!&=\phi_{1}(\gamma\phi^{2}-1)\,.\end{array}\right.$
(71)

The set of PDEs in (63) and the following expressions in Eqs. (68), (71) are
one of the main results of this manuscript. Overall, this represents a
considerable simplification from the original partial integrodifferential
kinetic equation (26) governing the dynamics of the inter-world distribution
function $F_{\alpha\beta}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}})$ with
the collision integral in Eq. (48).

To provide a first intuitive understanding of the previous set of PDEs, let us
briefly neglect spatial inhomogeneities of $\phi$ and $\boldsymbol{\phi}_{1}$,
and work in $d=1$. The equation on $\phi(t)$ becomes an autonomous first-order
ODE, reading

$\displaystyle\partial_{\tau}\phi$
$\displaystyle=\,\phi(\phi^{2}-1)=-V^{\prime}(\phi)\,.$ (72)

This is a gradient descent in the potential
$V(\phi)=\frac{1}{2}\phi^{2}-\frac{1}{4}\phi^{4}$. Reinstating the original
units, the rate of escape from the correlated-world solution at the unstable
extremum $\phi=1$, to the uncorrelated-world solution at the global minimum
$\phi=0$ is $t_{*}:=-\frac{1}{2}\tau_{\rm F}\log\delta\phi_{0}$ where
$\delta\phi_{0}:=1-\phi(0)\ll 1$. At early times, the growth of the
perturbation is exponential,

$\displaystyle 1-\phi(t\ll t_{*})\simeq\exp\left[2(t-t_{*})/\tau_{\rm
F}\right]\,,$ (73)

while it saturates at late times,

$\displaystyle\phi(t\gg
t_{*})\simeq\frac{1}{\sqrt{2}}\exp\left[-(t-t_{*})/\tau_{\rm F}\right]\to
0\,.$ (74)

Again, the decoupling of $\boldsymbol{\phi}_{1}$ in such a spatially
homogeneous setting is evidence that the momentum anisotropy captured by
$\boldsymbol{\phi}_{1}$ is a minimal ingredient necessary to allow for spatial
propagation of the relaxation from the correlated-world solution to the
uncorrelated-world solution. This will be the topic of Sect. 5.

In the general case, i.e. in the presence of spatial inhomogeneities, it is
instructive to compare the inter-world situation to the (standard) intra-world
kinetic equations. When an intra-world distribution function $F$ is associated
with a conserved quantity (e.g. number of particles or energy), the
corresponding hydrodynamic equation is typically expected to display diffusive
behavior. Indeed, the timescale associated with the conserved quantity is much
slower than the other modes: those can be effectively replaced by their local-
equilibrium value in a fixed background of $F$, typically resulting in a
diffusive term of the type $\boldsymbol{\nabla}^{2}_{\boldsymbol{X}}F$. Here,
in the inter-world case, the distribution function $F_{\alpha\beta}$ is not
associated with a conserved quantity (until proven otherwise) and there is no
clear separation of timescales in the PDEs (63) governing the dynamics of
$\phi$ and $\boldsymbol{\phi}_{1}$. Consequently, one cannot a priori apply
the standard hydrodynamic approach, and one must solve for the dynamics $\phi$
and $\boldsymbol{\phi}_{1}$ on an equal footing.

###  4.3 Validating the ansatz numerically

We perform two independent checks of the ansatz proposed in Eq. (49) by
comparing, on the one hand, the solutions of the inter-world kinetic equation
in (26) computed with the full-fledged collision integral in Eq. (48) with, on
the other hand, the solutions of the coupled PDEs in (63). The first check is
performed at early times in a near-critical regime while the second check is
performed at later times and at a finite distance from criticality.

Numerically solving the kinetic equation is a formidable task which we
simplify as much as possible by working in one dimension, $d=1$, with a
regular lattice dispersion $\epsilon_{k}=-\cos(k)$ for $k\in[-\pi,\pi)$, and a
point-like Fermi surface located at the wave-vector $k_{\rm F}=\pi/2$. We
measure energies in units of the half-bandwidth. Note that superconductivity
in $d=1$ is known to be quite different from dimensions $d\geq 2$, with the
RPA treatment that we have set up in Section 3 not suited for $d=1$. However,
the objective here is to put to test the ansatz in conditions that are
qualitatively similar to $d\geq 2$, and not to correctly capture the peculiar
one-dimensional physics. This is the reason why we can afford to work in
$d=1$. We further simplify the computation by working in a non-self-consistent
scheme, with the quasi-particle retarded Green’s function reading

$\displaystyle G^{R}(\omega,k)=\frac{1}{\omega-\epsilon_{k}+{\rm i}\Gamma}\,,$
(75)

and where the Cooperon Green’s function $D^{R}(\omega,k)$ is computed
following Eq. (36). $\Gamma>0$ sets a bare fermionic inverse lifetime. In
practice, $\Gamma$ helps the numerical convergence of our algorithms and we
set it as the smallest energy scale in the problem. Finally, we further reduce
the difficulty by working with on-shell quantities:
$\tilde{F}_{\alpha\beta}(k):=F_{\alpha\beta}(\omega=\epsilon_{k},k)$.

####  4.3.1 Early times

$0$$0.5$$1$$-0.3$$-0.2$$-0.1$$0$$0.1$$0.2$$0.3$(a)$\tau=0$

$\phi$

$X$

$0$$0.5$$1$$-0.3$$-0.2$$-0.1$$0$$0.1$$0.2$$0.3$(b)

$\phi$

$X$

Figure 4: (a) $\phi(\tau,X)$ extracted from the solution
$\tilde{F}_{du}(\tau,X;k_{\rm F})$ to the full-fledged 1D kinetic equation and
plotted at different times $\tau=0,0.02,0.04,\ldots,0.3$. The initial
condition is given in Eq. (76) with $R_{0}=0.05$ and $\delta\phi_{0}=0.98$.
(b) Solution $\phi(\tau,X)$ to the coupled PDEs in (63) at $\gamma=1$ and with
the corresponding initial condition. The parameters are $U=-1$, $T=0.1$, and
$\Gamma=0.01$ (in units of the half-bandwidth), corresponding to a small
detuning from criticality $r=-1/(\rho_{{\rm F}}D^{R}(\omega=0,k=0))\approx
0.1$ defined in Eq. (44). Time and space have been rescaled according to Eq.
(65). No adjustable parameters were used.

The first test consists in numerically solving the kinetic equation (26) with
the full-fledged collision integral in Eq. (48) for as long as we can ensure
numerical convergence of the solutions. In practice, this is challenging and
we can only access early times, i.e. on the order of fractions of $\tau_{\rm
F}$. Therefore, in order to benchmark the ansatz in both regimes $\phi\sim 1$
and $\phi\sim 0$, we work with an initial condition that simultaneously spans
those two regimes. We choose an initial condition with a perturbation of
$\tilde{F}_{\alpha\beta}^{\rm corr}(k)$ in the shape of a Gaussian droplet of
large amplitude $\delta\phi_{0}\lesssim 1$ and width $R_{0}$, and localized
around $X=0$. Explicitly, rescaling time and space according to Eq. (65), we
take the following symmetric initial condition

$\displaystyle\tilde{F}_{\alpha\beta}(\tau=0,X;k)=\left[1-\delta\phi_{0}(X)\right]\tilde{F}_{\alpha\beta}^{\rm
corr}(k)\longleftrightarrow\left\\{\begin{array}[]{rl}\phi(\tau=0,X)\\!\\!\\!\\!&=1-\delta\phi_{0}(X),\\\
{\phi}_{1}(\tau=0,X)\\!\\!\\!\\!&=0\,,\end{array}\right.$ (76c)
$\displaystyle\mbox{ with
}\delta\phi_{0}(X)=\delta\phi_{0}\exp[-X^{2}/(2R_{0}^{2})]\Theta(3R_{0}-|X|)\,,$
(76d)

and where $\Theta(X)$ is the Heaviside step function. We found such a
Gaussian-shaped droplet, defined on the support $[-3R_{0},3R_{0}]$, to be
easier to time-evolve numerically than a semi-circular droplet with sharp
edges. The physical parameters are chosen such as to be close to criticality
(on the disordered side) and to obey the hierarchy $|U|,\epsilon_{\rm F}\gg
T\gg\Gamma$, where $\epsilon_{\rm F}$ is the Fermi energy.

In Fig. 4, we compare the solutions $\phi(\tau,X)$ of the corresponding
coupled PDEs with the solutions $\tilde{F}_{du}(\tau,X;k)$ of the full-fledged
kinetic equation. The comparison is made by extracting the first partial-wave
contributions according to Eq. (52) which, in 1D, simply reads

$\displaystyle\left\\{\begin{array}[]{rl}\textstyle\phi(\tau,X)\\!\\!\\!\\!&=\textstyle\frac{1}{2}\left[\tilde{F}_{du}(\tau,X;k_{\rm
F})+\tilde{F}_{du}(\tau,X;-k_{\rm F})\right],\\\\[6.0pt]
\textstyle\phi_{1}(\tau,X)\\!\\!\\!\\!&=\textstyle\frac{1}{2}\left[\tilde{F}_{du}(\tau,X;k_{\rm
F})-\tilde{F}_{du}(\tau,X;-k_{\rm F})\right]\,,\end{array}\right.$ (79)

up to the time $\tau=0.3$. The qualitative agreement is excellent. Notice the
splitting of the initial central perturbation into both a left-moving and a
right-moving front. We repeated this analysis in a wide range of near-critical
parameters and initial conditions and consistently found excellent agreement,
even at a finite distance from criticality, in the presence of fast bosonic
fluctuations. This validates the ansatz in Eq. (49) at early times.

####  4.3.2 Late times and partial-wave truncation

$0$$0.5$$1$$0$$5$$10$(a)

$\phi$

$\tau$$\phi(\tau)$from kinetic Eq.from coupled PDEs

$0$$0.05$$0.1$$0$$5$$10$(b)

$\phi_{1}$

$\tau$$\phi_{1}(\tau)$from kinetic Eq.from coupled PDEs

Figure 5:  (a) $\phi(\tau)$ extracted from the solution
$\tilde{F}_{du}(\tau;k_{\rm F})$ to the kinetic equation with the spatially
homogeneous and momentum anisotropic initial condition given by Eq. (80) with
$\phi_{0}=0.85$ and ${\phi_{1}}_{0}=0.1$, and compared with the solution
$\phi(\tau)$ to the coupled PDEs in Eq. (84). (b) Same for the anisotropic
component $\phi_{1}(\tau)$ where $\gamma=0.5$ is the only adjustable
parameter. The physical parameters are $U=-2$, $T=0.6$, and $\Gamma=0.2$,
corresponding to a sizable detuning from criticality $r=-1/(\rho_{{\rm
F}}D^{R}(\omega=0,k=0))\approx 3.9$. Time has been rescaled according to Eq.
(65).

Because of the difficulty to produce converged numerical solutions of the
kinetic equation at larger times, we resort to a simpler, yet non-trivial,
benchmark for the ansatz. Let us consider the case of a spatially homogeneous
initial condition but with non-zero anisotropic components in momentum space.
Explicitly, we take the initial condition

$\displaystyle\tilde{F}_{ab}(\tau=0;k)=\left[\phi_{0}+{\phi_{1}}_{0}{\rm
sign}(k)\right]\tilde{F}^{\rm corr}_{ab}(k),$ (80)

on the kinetic equation side and, correspondingly,

$\displaystyle\phi(\tau=0)=\phi_{0}\ \mbox{ and }\
\phi_{1}(\tau=0)={\phi_{1}}_{0},$ (81)

on the side of the coupled PDEs which now read

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\phi\\!\\!\\!\\!&=\phi(\phi^{2}-1),\\\
\partial_{\tau}{\phi}_{1}\\!\\!\\!\\!&={\phi}_{1}(\gamma\phi^{2}-1)\,.\end{array}\right.$
(84)

Note that these coupled ordinary-differential equations correspond to the
discussion around Eq. (72). Notably, Eq. (84) predicts that the relaxation
dynamics of $\phi$ (but not those of $\phi_{1}$) are independent of the
distance to criticality, parameterized by $\gamma$. Therefore, complementary
to the previous benchmark in Sect. 4.3.1, we test the ansatz at a finite
distance from criticality (on the disordered side) by choosing an off-critical
set of physical parameters $U$, $T$, and $\Gamma$.

In Fig. 5, we compare the solutions $\phi(\tau)$ and $\phi_{1}(\tau)$ to the
corresponding coupled PDEs with the solutions $\tilde{F}_{du}(\tau,k)$ to the
kinetic equation. The comparison is made by extracting the first partial-wave
contributions according to Eq. (79). The qualitative agreement is very good
from early times down to late times when the dynamics have converged to the
uncorrelated world solution. The agreement for the dynamics of
$\phi_{1}(\tau)$ was made by manually adjusting the off-critical value for
$\gamma$ given in the caption. This validates the partial-wave truncation
which is made in the ansatz.

##  5 Dynamics of information scrambling

In this Section, we solve the set of coupled PDEs (68) that effectively govern
the dynamics of information scrambling. We first discuss the early times, when
a regime of exponential growth takes place. Later, we solve the geometry for
the late-time traveling front. Finally, we address the saturation regime in
the bulk of the information light cone.

###  5.1 Early-time exponential growth

At early times, the solutions to the inter-world kinetic equation and to the
simplified coupled PDEs are expected to be strongly dependent on the system
parameters and the initial conditions.

Here, we solve the coupled PDEs in the linear regime around the correlated-
world solution, $\phi\approx 1$ and $\boldsymbol{\phi}_{1}\approx 0$. We
expect this linear regime to be all the more valid as the initial perturbation
will be small, hence taking a longer time to reach the non-linear regime. To
that end, we consider the droplet-shaped initial condition

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

where the perturbation $0<\delta\phi_{0}({\boldsymbol{X}})\ll 1$ is non-
vanishing on a small compact ball of radius $R_{0}$ around
${\boldsymbol{X}}=0$. Note that this is a different regime from the numerics
presented in Sect. 4.3.1 where the initial condition was probing the non-
linear regime. To simplify, we consider the case $d=1$. The linearized coupled
PDEs on $\delta\phi:=1-\phi$ and $\phi_{1}$ read

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\delta\phi-\partial_{X}\phi_{1}\\!\\!\\!\\!&=2\delta\phi,\\\
\partial_{\tau}\phi_{1}-\partial_{X}\delta\phi\\!\\!\\!\\!&=(\gamma-1)\phi_{1}\,.\end{array}\right.$
(90)

This yields the following linearized PDE on $\delta\phi(\tau,X)$

$\displaystyle\partial_{\tau}^{2}\delta\phi-(\gamma+1)\partial_{\tau}\delta\phi-\partial^{2}_{X}\delta\phi=2(1-\gamma)\delta\phi\,,$
(91)

with $\delta\phi(\tau=0,X)=\delta\phi_{0}(X)$ and
$\partial_{\tau}\delta\phi(\tau=0,X)=2\delta\phi_{0}(X)$. Integrating the
above equation over the whole space, or equivalently the first equation in
(90), and introducing the integrated perturbation $\delta M(\tau):=\int{\rm
d}X\delta\phi(\tau,X)$, we find an exponential growth of the perturbation,
echoing the onset of chaos:

$\displaystyle\delta M(\tau\ll\tau_{*})=\exp\left[2(\tau-\tau_{*})\right],$
(92)

with the typical timescale to escape from the unstable solution given by
$\tau_{*}:=-\frac{1}{2}\log\delta M_{0}$ where $\delta M_{0}:=\int{\rm
d}X\delta\phi_{0}(X)$. The corresponding growth rate, $\lambda_{\rm
L}=2/\tau_{\rm F}$ in the original units, can be interpreted as a Lyapunov
exponent. Note that the dependence on $\gamma$, the parameter quantifying the
distance to criticality, has dropped. A more sophisticated calculation
restricted to $1\ll\tau\ll\tau^{*}$ and $|X|\ll\tau^{3/4}$ yields the
following solution to Eq. (91)

$\displaystyle\delta\phi(\tau,X)\approx{\rm e}^{2(\tau-\tau_{*})}\frac{{\rm
e}^{-(3-\gamma)X^{2}/4\tau}}{\sqrt{4\pi\tau/(3-\gamma)}}\,.$ (93)

See Appendix B for the details of the computation. The above expression
involves an exponential growth and a diffusion kernel. It is similar to the
solution obtained in the context of the $\mathcal{O}(N)$ model in Ref. [8]
(see Eq. 1.13 therein). Let us provide a quick way to justify this solution.
We first absorb the exponential growth of $\delta\phi$ that was found in the
solution (92) by working in terms of $g(\tau,X):=$ ${\rm
e}^{-2\tau}\delta\phi(\tau,X)$. Using Eq. (91), it obeys

$\displaystyle\partial_{\tau}^{2}g+(3-\gamma)\partial_{\tau}g-\partial^{2}_{X}g=0\,,$
(94)

with $g(\tau=0,X)=\delta\phi_{0}(X)$ and $\partial_{\tau}g(\tau=0,X)=0$. Next,
given that the solution of the above PDE is expected to vary slowly with time,
especially at large times, we may neglect the $\partial_{\tau}^{2}g$ term.
This yields a simple diffusion equation which, once we switch back to working
with $\delta\phi$, is solved by the expression in (93). Note that the solution
in (93) can also be seen as the solution to the diffusion+growth equation:

$\displaystyle\partial_{\tau}\delta\phi-\frac{1}{3-\gamma}\partial_{X}^{2}\delta\phi=2\delta\phi,$
(95)

that appears, notably, in the context of branching Brownian motion where it
describes the evolution of the expected density of particles [33].

###  5.2 Late-time solutions and discontinuous front

Here, we access the late-time inter-world distribution by analytically solving
the coupled PDEs (68) in generic spatial dimensions $d$, simply assuming a
spherically-symmetric initial condition. In particular, we shall show that a
wavefront propagates at a constant butterfly velocity controlled by the Fermi
velocity. The wavefront acts as a light cone that separates two causally
disjoint regions: ahead of the front, the inter-world distribution is the
correlated-world solution, while behind the front the two worlds rapidly
decohere to the uncorrelated-world solution. Notably, at the front, the
distribution function develops a discontinuity.

####  5.2.1 Traveling front

The solution develops a traveling front located on a sphere of increasing
radius. Our goal is to compute its velocity and shape.

The first step is to notice that the late-time position of the front is, by
definition, far from the origin and we may neglect the $1/r$ term in the rhs
of Eq. (71). By doing so, we simply recover the $d=1$ equations. Indeed, given
the large radius of the sphere where the front is located, the problem is
locally flat in the non-radial directions. Therefore, we only need to work out
the $d=1$ case. To simplify the presentation of the computation, we work out
the critical case where $\gamma=1$. However, the generic case for $\gamma\neq
1$ can also be solved by similar techniques and we refer the reader to
Appendix C for the corresponding detailed computation. Introducing the fields

$\displaystyle\varphi:=\frac{\phi+\phi_{1}}{2}\quad\mbox{ and
}\quad\psi:=\phi-\phi_{1}\,,$ (96)

the coupled PDEs can be cast as

$\displaystyle\left\\{\begin{array}[]{l}\partial_{\tau}\varphi+\partial_{X}\varphi=\varphi\,(\phi^{2}-1),\\\
\partial_{\tau}\psi-\partial_{X}\psi=\psi\,(\phi^{2}-1),\end{array}\right.$
(99)

with the initial conditions,

$\displaystyle\left\\{\begin{array}[]{l}\varphi(\tau=0,X)=\phi_{0}(X)/2,\\\
\psi(\tau=0,X)=\phi_{0}(X)\,.\end{array}\right.$ (102)

The method of characteristics gives the implicit solutions

$\displaystyle\left\\{\begin{array}[]{l}\varphi(\tau,X)=\frac{1}{2}\phi_{0}(X-\tau)\rm
e^{-\int_{0}^{\tau}{\rm d}s\left[1-\phi^{2}(s,X-\tau+s)\right]},\\\
\psi(\tau,X)=\phi_{0}(X+\tau){\rm e}^{-\int_{0}^{\tau}{\rm
d}s\left[1-\phi^{2}(s,X+\tau-s)\right]},\end{array}\right.$ (105)

which yield

$\displaystyle\phi(\tau,X)$ $\displaystyle=\frac{1}{2}\phi_{0}(X-\tau){\rm
e}^{-\int_{0}^{\tau}{\rm
d}s\left[1-\phi^{2}(s,X-\tau+s)\right]}+\frac{1}{2}\phi_{0}(X+\tau){\rm
e}^{-\int_{0}^{\tau}{\rm d}s\left[1-\phi^{2}(s,X+\tau-s)\right]}\,.$ (106)

Let us now assume that a right-moving front, traveling at velocity 1 (in units
of $v_{\rm F}/\sqrt{d}$), develops at late times, i.e.

$\displaystyle\phi_{+}(\tau,X):=\phi(\tau,X+\tau)\stackrel{{\scriptstyle\tau\to\infty}}{{\to}}f_{+}(X)\,,$
(107)

pointwise. Naturally, there is also a symmetrical left-moving front. We start
from

$\displaystyle\phi_{+}(\tau,X)$ $\displaystyle=\frac{1}{2}\phi_{0}(X){\rm
e}^{-\int_{0}^{\tau}{\rm
d}s\left[1-\phi_{+}^{2}(s,X)\right]}+\frac{1}{2}\phi_{0}(X+2\tau){\rm
e}^{-\int_{0}^{\tau}{\rm d}s\left[1-\phi_{+}^{2}(s,X+2\tau-2s)\right]}$ (108)
$\displaystyle=\frac{1}{2}\phi_{0}(X){\rm e}^{-\int_{0}^{\tau}{\rm
d}s\left[1-\phi_{+}^{2}(s,X)\right]}+\frac{1}{2}\phi_{0}(X+2\tau){\rm
e}^{-\int_{0}^{\tau}{\rm d}s\left[1-\phi_{+}^{2}(\tau-s,X+2s)\right]}\,,$
(109)

where in the second line we performed a change of dummy variable $s\to\tau-s$.
In the long-time limit $\tau\to\infty$, we get

$\displaystyle f_{+}(X)$ $\displaystyle=\frac{1}{2}\phi_{0}(X){\rm
e}^{-\int_{0}^{\infty}{\rm
d}s\left[1-\phi_{+}^{2}(s,X)\right]}+\frac{1}{2}{\rm
e}^{-\int_{0}^{\infty}{\rm d}s\left[1-f_{+}^{2}(X+2s)\right]}\,,$ (110)

where we used the property $\phi_{0}(X+2\tau)\to 1$. Consistently with the Eq.
(106), we postulate111This assumption can be rigorously proven to be true in
$d=1$. that the front is such that

$\displaystyle f_{+}(X<R_{0})<1\quad\mbox{ and}\quad f_{+}(X>R_{0})=1\,.$
(111)

Let us now work with $X<R_{0}$. Hence,
$\lim\limits_{s\to\infty}1-\phi_{+}^{2}(s,X<R_{0})=1-f_{+}^{2}(X)>0$, which
implies that the first term in Eq. (110) is zero. After the change of
variables $X+2s\to u$, we are left with

$\displaystyle f_{+}(X<R_{0})$ $\displaystyle=\frac{1}{2}{\rm
e}^{-\frac{1}{2}\int_{X}^{R_{0}}{\rm d}u\left[1-f_{+}^{2}(u)\right]}\,.$ (112)

At $X\stackrel{{\scriptstyle X<R_{0}}}{{\longrightarrow}}R_{0}$, this yields
$f_{+}(R_{0})=1/2$. Given that we assume $f_{+}(X>R_{0})=1$, this signals a
discontinuity in $f_{+}(X)$ at $X=R_{0}$. Moreover, Eq. (112) implies that
$f_{+}(X)$ obeys

$\displaystyle f_{+}^{\prime}(X)$
$\displaystyle=\frac{1}{2}\left[1-f_{+}^{2}(X)\right]f_{+}(X),$ (113)

which can be solved by separation of variables, yielding the discontinuous
right-moving front with the shape

$\displaystyle f_{+}(X<R_{0})=\frac{1}{\sqrt{1+3{\rm
e}^{R_{0}-X}}}\stackrel{{\scriptstyle X\to
R_{0}^{-}}}{{\mapsto}}\frac{1}{2}\quad\mbox{ and}\quad f_{+}(X>R_{0})=1\,.$
(114)

Note that, here, the information on the precise shape of the initial
perturbation is lost. Except for a trivial spatial offset in the shape of the
front, $R_{0}$ drops out of the problem. This is expected for generic initial
conditions defined on a compact support $|X|<R_{0}$. Indeed, one can show that
the first non-vanishing derivative $\partial^{(n)}_{X}\phi_{0}(X=R_{0}^{-})>0$
is responsible for the generation of a first-order derivative
$f_{+}(X=R_{0}^{-})>0$ which grows exponentially with time, therefore creating
a discontinuity at large times. This independence of the steady-state solution
with respect to the initial condition is the hallmark of a universal solution.
Notably, as previously discussed in Sec. 5.1, $R_{0}$ and the precise shape of
the initial perturbation are however controlling the timescale for the
traveling front to form. Barring this point, the universal features of the
steady state may be safely accessed by sending $R_{0}\to 0^{+}$ after
$\tau\to\infty$. We obtain the right-moving traveling front

$\displaystyle\left\\{\begin{array}[]{cc}\lim\limits_{\tau\to\infty}\phi(\tau,X+\tau)\\!\\!\\!\\!&=f_{+}(X),\\\
\lim\limits_{\tau\to\infty}\phi_{1}(\tau,X+\tau)\\!\\!\\!\\!&=f_{1+}(X),\end{array}\right.$
(117)

with

$\displaystyle\left\\{\begin{array}[]{cc}f_{+}(X<0)=\frac{1}{\sqrt{1+3{\rm
e}^{-X}}}\stackrel{{\scriptstyle X\to 0^{-}}}{{\mapsto}}\frac{1}{2}\quad\mbox{
and}\quad f_{+}(X>0)=1,\\\ f_{1+}(X<0)=-\frac{1}{\sqrt{1+3{\rm
e}^{-X}}}\stackrel{{\scriptstyle X\to
0^{-}}}{{\mapsto}}-\frac{1}{2}\quad\mbox{ and}\quad
f_{1+}(X>0)=0\,.\end{array}\right.$ (120)

The last equation above follows from $\phi_{1}=\varphi-\psi/2$. As we
discussed above, the velocity and the precise shape of the late-time traveling
front, $f_{+}(X)$, and notably its discontinuity, generalize to the radial
front $f_{+}(r)$ in generic dimension $d$ (assuming spherically-symmetric
initial conditions). In terms of the original inter-world distribution
function, reinstating the original units of time and space, the traveling
front reads

$\displaystyle\lim\limits_{t\to\infty}\tilde{F}_{du}({\boldsymbol{k}};t,r+v_{\rm
B}t)\big{\lvert}_{k\to k_{\rm
F}}=\begin{cases}f_{+}\left(\frac{r\sqrt{d}}{\ell}\right)\left[1-{\boldsymbol{u}}_{\boldsymbol{k}}\cdot{\boldsymbol{u}}_{\boldsymbol{r}}\right]\mbox{
if }r<0,\\\ 1\mbox{ if }r>0,\end{cases}$ (121)

with the butterfly velocity given by $v_{\rm B}:={v_{\rm F}}/{\sqrt{d}}$ and
the radial front shape $f_{+}$ given by Eq. (120) that spans over a length
scale controlled by the mean free path $\ell:=v_{\rm F}\tau_{\rm F}$. We
recall that $v_{\rm F}$ is the Fermi velocity and the scattering time
$\tau_{\rm F}$ was defined in Eq. (64).

$0$$0.5$$1$$-30$$-20$$-10$$0$$10$$20$$30$(a)

$\phi$

$X$

$0$$10$$20$$30$$0$$10$$20$$30$(b)

$X$

$\tau$$X_{\rm front}=\tau+R_{0}$

Figure 6: (a) Solution $\phi(\tau,X)$ to the $d=1$ coupled PDEs in Eq. (68) at
$\gamma=1$ for the droplet initial perturbation in Eq. (125) with amplitude
$\delta\phi_{0}=0.02$ and radius $R_{0}=0.1$. Different times
$\tau=0,2,\ldots,26,28$ are plotted with different colors. (b) Location of the
front extracted from the data in (a) marked by crosses. The line is the
analytical prediction, with a constant butterfly velocity, $X_{\rm
front}(\tau)=\tau+R_{0}$. No adjustable parameter was used.
$0$$0.5$$1$$-10$$-5$$0$$5$$10$

$\phi$

$X$$\phi(\tau,X\\!-\\!X_{\rm front}(\tau))\\!$$f_{+}(X)\\!$ Figure 7:
Discontinuous late-time front: the solution $\phi(\tau,X-X_{\rm front}(\tau))$
is extracted from the data of Fig. 6 (a) for $\tau=28$ (solid line), and
compared to the analytical prediction $f_{+}(X)$ t in Eq. (120) for the
traveling front, with a discontinuity at $X=0$ from 0.5 to 1 (dashed line). No
adjustable parameter was used.

This solution can also be generalized to cases away from criticality, i.e. for
a generic $\gamma\neq 1$. The qualitative features are found to be very
similar to the critical case $\gamma=1$. In the Appendix C, we show that the
discontinuity of the front is now from $f_{+}(X=0^{-})=L(\gamma)$ to
$f_{+}(X>0)=1$ with

$\displaystyle L(\gamma)=\frac{\sqrt{5+4\gamma}-1}{2(1+\gamma)}\,.$ (122)

This quantity monotonically interpolates from $L=1/2$ for $\gamma=1$ to
$L=1/(\text{golden ratio})$ for $\gamma=0$.

The discontinuous traveling wavefronts in Eqs. (120), (121), and their
generalization to non-critical regimes in Eq. (122) are one of the main
results of this manuscript. To illustrate the analytical solution, we compare
its predictions to the numerical solution of the coupled PDEs at $\gamma=1$
(critical regime). In Fig. 6, we display the numerical solution starting from
the following droplet-shaped initial perturbation of the correlated-world
solution

$\displaystyle\left\\{\begin{array}[]{rl}\phi(\tau=0,X)\\!\\!\\!\\!&=1-\delta\phi_{0}(X)\mbox{
with
}\delta\phi_{0}(X)=\delta\phi_{0}\sqrt{1-(X/R_{0})^{2}}\Theta(R_{0}-|X|),\\\
{\phi}_{1}(\tau=0,X)\\!\\!\\!\\!&=0\,,\end{array}\right.$ (125)

where $\delta\phi_{0}$ sets the amplitude of the droplet, and $R_{0}$ sets its
radius. This illustrates unambiguously the spatial growth of the loss of
quantum coherence as time goes on. The light cone structure of this growth,
with a front traveling at a constant velocity, is demonstrated by simply
extracting the location of the front as a function of time. In Fig. 7, we
illustrate the discontinuous shape of the traveling front by superimposing the
front extracted from the numerical solution of $\phi(\tau,X)$ at the late-time
$\tau=28$ to the exact result given by $f_{+}(X)$ in Eq. (120).

####  5.2.2 Saturation inside the light cone

Inside the light cone and far enough from its boundaries, we have $\phi\ll 1$
and we can neglect the non-linearities in the rhs of the coupled PDEs (99).
Working in the $d=1$ case, the linearized PDEs read

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\phi+\partial_{X}\phi_{1}\\!\\!\\!\\!&=-\phi,\\\
\partial_{\tau}\phi_{1}+\partial_{X}\phi\\!\\!\\!\\!&=-\phi_{1}\,.\end{array}\right.$
(128)

Assuming symmetric initial conditions, i.e.,

$\phi(\tau=0,X)=\phi(\tau=0,-X)\mbox{ and
}\phi_{1}(\tau=0,X)=-\phi_{1}(\tau=0,-X)\,,$

one can simply show that this symmetry is preserved by the entire time
evolution (even in the presence of the non-linear terms that were at play in
the earlier regime). The solutions are of the form

$\displaystyle\left\\{\begin{array}[]{rl}\phi(\tau,X)\\!\\!\\!\\!&={\rm
e}^{-X}f_{\varphi}(X-\tau)+{\rm e}^{X}f_{\varphi}(-X-\tau),\\\
\phi_{1}(\tau,X)\\!\\!\\!\\!&={\rm e}^{-X}f_{\varphi}(X-\tau)-{\rm
e}^{X}f_{\varphi}(-X-\tau),\end{array}\right.$ (131)

where the function $f_{\varphi}(X)$ should in principle be determined by
solving the early-to-intermediate time problem. In practice, we can determine
the asymptotic behavior of the function $f_{\varphi}(X)$ at large $X$ by
requiring a matching to the left side of the late-time front computed
previously:

$\displaystyle\lim\limits_{\tau\to\infty}\phi(\tau,X+\tau)\sim f_{+}(X)\mbox{
for }X<0\,,|X|\gg 1,$ (132)

where $f_{+}(X)$ was computed in Eq. (120). Using the late-time solution in
Eq. (131), we have

$\displaystyle\lim\limits_{\tau\to\infty}{\rm e}^{-X-\tau}f_{\varphi}(X)+{\rm
e}^{X+\tau}f_{\varphi}(-X-2\tau)\sim f_{+}(X)\,.$ (133)

The first term vanishes and we are left with

$\displaystyle\lim\limits_{\tau\to\infty}{\rm
e}^{X+\tau}f_{\varphi}(-X-2\tau)\sim f_{+}(X)\,,$ (134)

which is solved as

$\displaystyle f_{\varphi}(X)$ $\displaystyle\sim{\rm e}^{X/2}\,.$ (135)

This yields, at large times and far inside the boundary of the light cone,
i.e. for $\tau-|X|\gg 1$,

$\displaystyle\left\\{\begin{array}[]{rl}\phi(\tau,X)\\!\\!\\!\\!&\simeq 2{\rm
e}^{-\tau/2}\cosh{(X/2)},\\\ \phi_{1}(\tau,X)\\!\\!\\!\\!&\simeq-2{\rm
e}^{-\tau/2}\sinh{(X/2)}\,.\end{array}\right.$ (138)

These solutions can also be obtained directly from injecting the expression
(120) for the traveling front in Eq. (105), using the spatial symmetry
$\lim\limits_{\tau\to\infty}\phi(\tau,\tau+X)=\lim\limits_{\tau\to\infty}\phi(\tau,-\tau-X)=f_{+}(X)$.
Note that these asymptotic solutions have already lost the information about
the initial condition. Moreover, they imply the relation
$\phi_{1}(\tau,X)=-2\partial_{X}\phi(\tau,X)$ for $\tau-|X|\gg 1$. Once re-
injected in the linearized PDEs, this yields the following linearly-driven
diffusion equation for $\phi(\tau,X)$ only

$\displaystyle\partial_{\tau}\phi-2\partial_{X}^{2}\phi=-\phi\,,$ (139)

where the diffusive constant in front of the $\partial_{X}^{2}\phi$ term is
$D:=2v_{\rm F}^{2}\tau_{F}$ once the original units of time and space are
reinstated. We note that the diffusive lhs of Eq. (139) is similar to the one
of the FKPP-like equation that was derived in a similar context in Ref. [15],
see Eq. (76) therein.

##  6 Discussion and conclusion

Starting from the microscopic Hamiltonian of a $d$-dimensional quantum many-
body system of interacting electrons close to a superconducting phase
transition, we carefully derived the corresponding dynamics of quantum
information scrambling.

Quite expectedly, we found a ballistic spread of information governed by a
non-universal butterfly velocity $v_{\rm B}$. We presented analytical
solutions in the different regimes relevant to quantum chaotic dynamics: the
early exponential growth, the geometry of the late-time front, and the
saturation within the light cone.

Perhaps the most striking result of our work is the fact that the scrambling
of quantum information at late times is governed by shock-wave dynamics.
Scrambling propagates at the maximum velocity allowed by causality and
develops a distinct discontinuity exactly at the boundary of the light cone.
Notably, these dynamics scrupulously respect causality: scrambling does not
leak outside of the light cone. This is different from the findings of
previous works in similar settings, which often exhibited sharp but continuous
fronts preceded by exponential tails. At a formal level, this difference
arises from the fact that our effective dynamics are governed by a set of
coupled PDEs that does not belong to the reaction-diffusion class. While our
full solution does share strong similarities with the diffusive FKPP solutions
in the linearized regimes (at early times or deep in the bulk of the light
cone), we attribute this explicit absence of a diffusive term to the fact that
information dynamics is not directly associated with a conserved quantity,
unlike usual transport which is associated with, say, charge, energy, or
momentum conservation.

Addressing the robustness of the finite spacetime discontinuity in the
traveling front is perhaps one of the most pressing questions raised by our
results. Let us first note that the discontinuity has to be understood on the
scale of the mean free path set by the inelastic scattering, $\ell:=v_{\rm
F}\tau_{\rm F}$, and not on the microscopic scale $1/k_{\rm F}$ where the
kinetic equation breaks down. The discontinuity is unlikely an artifact of the
one-loop RPA scheme (exact in the limit where the number of electronic
orbitals is sent to infinity), nor of our truncated partial-wave expansion in
momentum space (key to the ensuing strict causal structure, and likely to be a
very good approximation in some geometries). However, higher-order terms in
the Moyal product expansion that were neglected in the quasi-classical
approximation, or any source of noise or disorder [34], e.g. caused by elastic
scattering on random impurities, could smoothen the front at the light-cone
boundary and possibly decrease the butterfly velocity. Numerical confirmation
of the shock-wave dynamics starting from a microscopic Hamiltonian is expected
to be difficult as the scrambling front propagates fast and the discontinuity
only develops at late times. This means that one should simulate large systems
of linear size $L\sim 10\,\ell$ until long times $t_{\rm max}\sim
10\,\tau_{{\rm F}}$. Tensor network approaches in 1d may be suited to meet the
challenge [4].

## Acknowledgements

The authors are grateful to Adam Nahum and Tianci Zhou for insightful
discussions. AM thanks the hospitality of École Normale Supérieure (Paris)
where this work was initiated.

##### Funding information

CA acknowledges the support from the French ANR “MoMA” project
ANR-19-CE30-0020 and the support from the project 6004-1 of the Indo-French
Centre for the Promotion of Advanced Research (IFCPAR). This work was
supported by the US National Science Foundation Grant NSF-DMR 2018358 (AM). AM
acknowledges the Aspen Center for Physics where part of this work was
performed and which is supported by the US National Science Foundation Grant
PHY-1607611.

##  Appendix A Inter-world collision integral

In this Appendix, we carefully derive the contributions of the interworld
collision integral to the coupled PDEs (63). We recall the collision integral,
onshell:

$\displaystyle\tilde{I}_{\alpha\beta}({\boldsymbol{k}})$
$\displaystyle=2\,{\rm
Im}\,\tilde{\Sigma}^{R}({\boldsymbol{k}})\,\tilde{F}_{\alpha\beta}({\boldsymbol{k}})+{\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}}),$ (140)

where

$\displaystyle 2\,{\rm Im}\,\tilde{\Sigma}^{R}({\boldsymbol{k}})$
$\displaystyle=2\sum_{{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}$
$\displaystyle\qquad\times\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})$
$\displaystyle\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_{\boldsymbol{k}}-\omega^{\prime}}{2T}\right)\right],$
(141)

and

$\displaystyle{\rm i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})$
$\displaystyle=2\sum_{{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})$
$\displaystyle\qquad\times\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\tilde{F}_{\alpha\beta}({\boldsymbol{k}}^{\prime\prime})\tilde{F}_{\alpha\beta}({\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})\tilde{F}_{\beta\alpha}({\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}).$
(142)

Now let us inject the following ansatz in Eq. (140),

$\displaystyle\tilde{F}^{\rm ansatz}_{\alpha\beta}({\boldsymbol{k}})$
$\displaystyle=\left[\phi+{\boldsymbol{u}}_{\boldsymbol{k}}\cdot\boldsymbol{\phi}_{1}\right]\tilde{F}_{\alpha\beta}^{\rm
corr}({\boldsymbol{k}}).$ (143)

The term in ${\rm Im}\tilde{\Sigma}^{R}({\boldsymbol{k}})$ is straight-
forward:

$\displaystyle 2\,{\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}})\,\left[\phi+{\boldsymbol{u}}_{\boldsymbol{k}}\cdot\boldsymbol{\phi}_{1}\right]\tilde{F}_{\alpha\beta}^{\rm
corr}({\boldsymbol{k}}).$ (144)

Let us treat the term ${\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})$ in Eq. (142). It
produces terms in $\phi^{3}$, $\phi^{2}\phi_{1}$, $\phi\phi_{1}^{2}$, and
$\phi_{1}^{3}$. In practice, consistently with our choice of ansatz which
consists of tracking only the two first multipolar contributions to
$\tilde{F}({\boldsymbol{k}})$, we discard the terms of order $\phi_{1}^{2}$
and $\phi_{1}^{3}$ which yield higher-order multipolar contributions to the
collision integral. The term in $\phi^{3}$ reads

$\displaystyle
2\phi^{3}\sum_{{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})$
$\displaystyle\qquad\qquad\times\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}}^{\prime\prime})\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})\tilde{F}^{\rm
corr}_{\beta\alpha}({\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})$
$\displaystyle=2\phi^{3}\,{\rm Im}\tilde{\Sigma}^{R}({\boldsymbol{k}}),$ (145)

where we used the trigonometric relation

$\displaystyle F^{\rm corr}_{\alpha\beta}(\omega^{\prime\prime})F^{\rm
corr}_{\alpha\beta}(\omega^{\prime}-\omega^{\prime\prime})F^{\rm
corr}_{\beta\alpha}(\omega^{\prime}-\omega)$
$\displaystyle+\underbrace{\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]}_{\geq
0}F^{\rm corr}_{\alpha\beta}(\omega)=0.$ (146)

Let us now evaluate the term of order $\phi^{2}\phi_{1}$. We have

$\displaystyle{\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})=-2\phi^{3}\,{\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}})F_{\alpha\beta}^{\rm
corr}({\boldsymbol{k}})$
$\displaystyle+2\phi^{2}\sum_{{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})$
(147)
$\displaystyle\times\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\biggl{\\{}\boldsymbol{\phi}_{1}\cdot({\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}})\biggr{\\}}F^{\rm
corr}_{\alpha\beta}(\omega^{\prime\prime})F^{\rm
corr}_{\alpha\beta}(\omega^{\prime}-\omega^{\prime\prime})F^{\rm
corr}_{\beta\alpha}(\omega^{\prime}-\omega).$

The expression can be simplified by use of Eq. (A), yielding

$\displaystyle{\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})=-2\phi^{3}\,{\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}})\tilde{F}_{\alpha\beta}^{\rm
corr}({\boldsymbol{k}})$ $\displaystyle-2\phi^{2}\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}})\sum_{{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})$
$\displaystyle\times\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\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_{\boldsymbol{k}}-\omega^{\prime}}{2T}\right)\right]$
$\displaystyle\times\biggl{\\{}\boldsymbol{\phi}_{1}\cdot({\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}})\biggr{\\}}.$
(148)

###  A.1 Critical case

Close to criticality, the Cooperon propagator reads

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

with the positive parameters $r\propto(T-T_{\rm c})/T_{\rm c}$,
$a\sim\mathcal{O}(1)$, $\xi^{2}\sim v_{\rm F}^{2}/T^{2}$. At criticality $r\to
0$, the Cooperon becomes soft with a diverging length scale $l\sim 1/r^{\nu}$
(here $\nu=1/2$), and the propagator is singular at $\omega=k=0$. In this
case, we can approximate the term
${\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}}\approx-{\boldsymbol{u}}_{\boldsymbol{k}}$.
Thus we have

$\displaystyle{\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})=-2\phi^{3}\,{\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}})F_{\alpha\beta}^{\rm
corr}({\boldsymbol{k}})$ (150)
$\displaystyle+2\phi^{2}(\boldsymbol{\phi}_{1}\cdot{\boldsymbol{u}}_{\boldsymbol{k}})\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}})\sum_{{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})$
$\displaystyle\times\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\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_{\boldsymbol{k}}-\omega^{\prime}}{2T}\right)\right].$

Performing the sums on $\omega^{\prime\prime}$ and
${\boldsymbol{k}}^{\prime\prime}$, we get

$\displaystyle{\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})=-2\phi^{3}\,{\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}})F_{\alpha\beta}^{\rm
corr}({\boldsymbol{k}})$
$\displaystyle+4\phi^{2}(\boldsymbol{\phi}_{1}\cdot{\boldsymbol{u}}_{\boldsymbol{k}})\sum_{{\boldsymbol{k}}^{\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}{\rm
Im}\Pi^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})\left[\coth\left(\frac{\omega^{\prime}}{2T}\right)+\tanh\left(\frac{\epsilon_{\boldsymbol{k}}-\omega^{\prime}}{2T}\right)\right]$
$\displaystyle\times\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}}).$ (151)

The above can now be written as

$\displaystyle{\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})=-2\phi^{2}\left[\phi+(\boldsymbol{\phi}_{1}\cdot{\boldsymbol{u}}_{\boldsymbol{k}})\right]\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}}){\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}}).$ (152)

Altogether, we obtain

$\displaystyle\tilde{I}_{\alpha\beta}({\boldsymbol{k}})=2(1-\phi^{2})\underbrace{\left[\phi+(\boldsymbol{\phi}_{1}\cdot{\boldsymbol{u}}_{\boldsymbol{k}})\right]\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}})}_{\tilde{F}^{\rm
ansatz}_{\alpha\beta}({\boldsymbol{k}})}\,{\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}}).$ (153)

###  A.2 Away from criticality

In the Subsection A.1 above, we have treated the critical case which yields
the following coupled PDEs

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

where $v_{{\rm F}}$ is the Fermi velocity and we defined the timescale
$\tau_{\rm F}$ as $1/\tau_{\rm F}:=-2\,\mathrm{Im}\,\tilde{\Sigma}^{R}(k_{\rm
F})$.

Away from criticality, assuming $1/\tau_{\rm F}>0$, we separate the expression
of ${\rm i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})$ in Eq. (148)
into the critical expression computed in Eq. (152) and the rest by simply
writing

$\displaystyle{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{\boldsymbol{k}}=-{\boldsymbol{u}}_{\boldsymbol{k}}+({\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{\boldsymbol{k}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}}).$
(157)

Explicitly, we have

$\displaystyle{\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}({\boldsymbol{k}})=-2\phi^{2}\left[\phi+(\boldsymbol{\phi}_{1}\cdot{\boldsymbol{u}}_{\boldsymbol{k}})\right]\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}}){\rm
Im}\tilde{\Sigma}^{R}({\boldsymbol{k}})$
$\displaystyle-2\phi^{2}\tilde{F}^{\rm
corr}_{\alpha\beta}({\boldsymbol{k}})\sum_{{\boldsymbol{k}}^{\prime}{\boldsymbol{k}}^{\prime\prime}}\int\frac{{\rm
d}\omega^{\prime}}{2\pi}\frac{{\rm
d}\omega^{\prime\prime}}{2\pi}|D^{R}(\omega^{\prime},{\boldsymbol{k}}^{\prime})|^{2}\,\mathrm{Im}\,G^{R}(\omega^{\prime\prime},{\boldsymbol{k}}^{\prime\prime})\,\mathrm{Im}\,G^{R}(\omega^{\prime}-\omega^{\prime\prime},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime})$
$\displaystyle\times\mathrm{Im}\,G^{R}(\omega^{\prime}-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}})\underbrace{\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_{\boldsymbol{k}}-\omega^{\prime}}{2T}\right)\right]}_{\geq
0}$
$\displaystyle\times\biggl{\\{}\boldsymbol{\phi}_{1}\cdot\left({\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}^{\prime\prime}}+{\boldsymbol{u}}_{\boldsymbol{k}}+{\boldsymbol{u}}_{{\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}}\right)\biggr{\\}},$
(158)

where the first term is the critical expression while the second term collects
the rest. Assuming that
$\epsilon_{\boldsymbol{k}}=\epsilon_{-{\boldsymbol{k}}}$, one may check that
this second term is odd under ${\boldsymbol{k}}\to-{\boldsymbol{k}}$ and
therefore cannot contribute to the projection on the momentum-space isotropic
contribution

$\displaystyle\int{\rm d}\Omega_{\boldsymbol{k}}\mbox{ last term of
Eq.~{}(\ref{eq:Sigma_K_away_from_crit})}=0.$ (159)

This guarantees that the rhs of the first equation in the coupled PDE in (156)
is also valid away from criticality. However, in the absence of a similar
symmetry argument, we expect the projection to the first partial wave to be
non-vanishing, i.e.

$\displaystyle\int{\rm
d}\Omega_{\boldsymbol{k}}{\boldsymbol{u}}_{\boldsymbol{k}}\mbox{ last term of
Eq.~{}(\ref{eq:Sigma_K_away_from_crit})}\neq 0,$ (160)

and therefore to give an extra contribution to the term in
$\phi^{2}\boldsymbol{\phi}_{1}$ in the rhs of the second equation in (156).
Now working at the Fermi surface, we can parameterize, without loss of
generality, the amplitude of this contribution relative to the critical case
by use of the dimensionless quantity $\gamma>0$:

$\displaystyle\frac{d}{S_{d-1}}\int{\rm
d}\Omega_{\boldsymbol{k}}{\boldsymbol{u}}_{\boldsymbol{k}}{\rm
i}\tilde{\Sigma}^{K}_{\alpha\beta}(k_{\rm
F}{\boldsymbol{u}}_{\boldsymbol{k}})=\gamma\boldsymbol{\phi}_{1}\phi^{2}\tilde{F}^{\rm
corr}_{\alpha\beta}(k_{\rm F})/\tau_{{\rm F}}.$ (161)

This justifies the rhs in the second line of the coupled PDEs in (63).
$\gamma=1$ corresponds to the critical case and we expect the near-critical
regime to be described by $\gamma<1$.

##  Appendix B Early-time exponential growth

In this Appendix, we solve the early-time regime of the coupled PDEs (68) in
$d=1$,

$\begin{cases}\partial_{\tau}\phi+\partial_{X}\phi_{1}=\phi(\phi^{2}-1)\\\
\partial_{\tau}\phi_{1}+\partial_{X}\phi=\phi_{1}(\gamma\phi^{2}-1),\end{cases}$
(162)

with $0<\gamma\leq 1$ and the initial conditions assumed to be $C^{1}$ at
least,

$\begin{cases}\phi(\tau=0,X)=1-\delta\phi_{0}(X)\text{ with
}\delta\phi_{0}(|X|>R_{0})=0\\\ \phi_{1}(\tau=0,X)=0.\end{cases}$ (163)

Let us first introduce

$\begin{cases}h:=1-\phi-\phi_{1}=\delta\phi-\phi_{1}\\\
k:=1-\phi+\phi_{1}=\delta\phi+\phi_{1},\end{cases}$ (164)

to obtain the linearized equations

$\begin{cases}\displaystyle\partial_{\tau}h+\partial_{X}h=\frac{1+\gamma}{2}h+\frac{3-\gamma}{2}k\\\\[6.0pt]
\displaystyle\partial_{\tau}k-\partial_{X}k=\frac{1+\gamma}{2}k+\frac{3-\gamma}{2}h,\end{cases}$
(165)

with initial condition $h(\tau=0,X)=k(\tau=0,X)=\delta\phi_{0}(X)$. Writing

$\begin{cases}h(\tau,X)={\rm
e}^{\frac{1+\gamma}{2}\tau}\big{[}\delta\phi_{0}(X-\tau)+\frac{3-\gamma}{4}\int_{-\tau}^{\tau}{\rm
d}s\,\delta\phi_{0}(X+s)a(s,\tau)\big{]}\\\\[6.0pt] k(\tau,X)={\rm
e}^{\frac{1+\gamma}{2}\tau}\big{[}\delta\phi_{0}(X+\tau)+\frac{3-\gamma}{4}\int_{-\tau}^{\tau}{\rm
d}s\,\delta\phi_{0}(X+s)b(s,\tau)\big{]},\end{cases}$ (166)

we obtain by direct substitution

$\begin{cases}\partial_{\tau}h+\partial_{X}h-\frac{1+\gamma}{2}h\\!\\!\\!\\!&=\frac{3-\gamma}{2}{\rm
e}^{\frac{1+\gamma}{2}\tau}\big{[}a(\tau,\tau)\delta\phi_{0}(X+\tau)+\frac{1}{2}\int_{-\tau}^{\tau}{\rm
d}s\,\delta\phi_{0}(X+s)(\partial_{\tau}-\partial_{s})a(s,\tau)\big{]}\\\\[6.0pt]
\partial_{\tau}k-\partial_{X}k-\frac{1+\gamma}{2}k\\!\\!\\!\\!&=\frac{3-\gamma}{2}{\rm
e}^{\frac{1+\gamma}{2}\tau}\big{[}b(-\tau,\tau)\delta\phi_{0}(X-\tau)+\frac{1}{2}\int_{-\tau}^{\tau}{\rm
d}s\,\delta\phi_{0}(X+s)(\partial_{\tau}+\partial_{s})b(s,\tau)\big{]},\end{cases}$
(167)

so that the solution to Eq. (165) is obtained if

$\begin{cases}\displaystyle
a(\tau,\tau)=1,\qquad(\partial_{\tau}-\partial_{s})a(s,\tau)=\frac{3-\gamma}{2}b(s,\tau)\\\\[6.0pt]
\displaystyle
b(-\tau,\tau)=1,\qquad(\partial_{\tau}+\partial_{s})b(s,\tau)=\frac{3-\gamma}{2}a(s,\tau).\end{cases}$
(168)

In turn, the solution to Eq. (168) is

$\begin{cases}\displaystyle
a(s,\tau)=I_{0}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}+I_{1}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}\sqrt{\frac{\tau-s}{\tau+s}}\\\\[6.0pt]
\displaystyle
b(s,\tau)=I_{0}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}+I_{1}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}\sqrt{\frac{\tau+s}{\tau-s}},\end{cases}$
(169)

with $I_{0}$ and $I_{1}$ the modified Bessel functions of the first kind.
Indeed, $a(\tau,\tau)=b(-\tau,\tau)=1$ since $I_{0}(0)=1$, and we check the
differential equations in Eq. (168) hold. For the first one:

$\displaystyle(\partial_{\tau}-\partial_{s})a(s,\tau)$
$\displaystyle=\frac{3-\gamma}{2}\left[\sqrt{\frac{\tau+s}{\tau-s}}I_{0}^{\prime}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}+I_{1}^{\prime}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}\right]+\frac{I_{1}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}}{\sqrt{\tau^{2}-s^{2}}}$
$\displaystyle=\frac{3-\gamma}{2}b(s,\tau).$ (170)

where we used

$I_{0}^{\prime}(z)=I_{1}(z)\mbox{ and
}I_{1}^{\prime}(z)=I_{0}(z)-\frac{1}{z}I_{1}(z).$ (171)

Similarly,

$\displaystyle(\partial_{\tau}+\partial_{s})b(s,\tau)$
$\displaystyle=\frac{3-\gamma}{2}\left[\sqrt{\frac{\tau-s}{\tau+s}}I_{0}^{\prime}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}+I_{1}^{\prime}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}\right]+\frac{I_{1}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}}{\sqrt{\tau^{2}-s^{2}}}$
$\displaystyle=\frac{3-\gamma}{2}a(s,\tau).$ (172)

Back in terms of $\delta\phi$, the exact solution of the linearized equations
is

$\displaystyle\delta\phi(\tau,X)$
$\displaystyle=\frac{h(\tau,X)+k(\tau,X)}{2}$ $\displaystyle={\rm
e}^{\frac{1+\gamma}{2}\tau}\bigg{[}\frac{\delta\phi_{0}(X-\tau)+\delta\phi_{0}(X+\tau)}{2}$
(173)
$\displaystyle\qquad\quad+\frac{3-\gamma}{4}\\!\int_{-\tau}^{\tau}\\!\\!\\!\\!{\rm
d}s\,\delta\phi_{0}(X+s)\bigg{(}I_{0}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}+I_{1}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-s^{2}}\Big{)}\frac{\tau}{\sqrt{\tau^{2}-s^{2}}}\bigg{)}\bigg{]}$

We now use that $\delta\phi_{0}$ is 0 except on a ball of radius $R_{0}$
around 0 (with $R_{0}$ of order 1, at most), and we call $\delta
M_{0}:=\int{\rm d}X\delta\phi_{0}(X)$. The first terms in Eq. (173) vanish
except around $X=\pm\tau$. The integrand is non-zero only around $s=-X$.
Hence, for $\tau-|X|\gg 1$, we obtain

$\delta\phi(\tau,X)\simeq\delta M_{0}\frac{3-\gamma}{4}{\rm
e}^{\frac{1+\gamma}{2}\tau}\bigg{(}I_{0}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-X^{2}}\Big{)}+I_{1}\Big{(}\frac{3-\gamma}{2}\sqrt{\tau^{2}-X^{2}}\Big{)}\frac{\tau}{\sqrt{\tau^{2}-X^{2}}}\bigg{)}.$
(174)

Using

$I_{0}(z\gg 1)\simeq I_{1}(z\gg 1)\simeq\frac{{\rm e}^{z}}{\sqrt{2\pi z}}$
(175)

gives, for $\tau-|X|\gg 1$,

$\delta\phi(\tau,X)\simeq\delta
M_{0}\frac{\sqrt{3-\gamma}}{2\sqrt{4\pi}(\tau^{2}-X^{2})^{3/4}}{\rm
e}^{\frac{1+\gamma}{2}\tau+\frac{3-\gamma}{2}\sqrt{\tau^{2}-X^{2}}}\bigg{(}\tau+{\sqrt{\tau^{2}-X^{2}}}\bigg{)}.$
(176)

Recalling that
$\sqrt{\tau^{2}-X^{2}}=\tau-X^{2}/(2\tau)+\mathcal{O}(X^{4}/\tau^{3})$, we see
in particular that for $\tau\gg 1$ and $|X|\ll\tau^{3/4}$, this yields the
following form for the solution of the linear regime

$\delta\phi(\tau,X)\simeq\delta M_{0}\,{\rm e}^{2\tau}\frac{{\rm
e}^{-(3-\gamma)\frac{X^{2}}{4\tau}}}{\sqrt{4\pi\tau/(3-\gamma)}},$ (177)

which is valid as long as $\delta\phi(\tau,X)\ll 1$, i.e. for
$\tau\ll\tau_{*}:=-\frac{1}{2}\log\delta M_{0}$.

##  Appendix C Discontinuous front away from criticality ($\gamma\neq 1$)

In this Appendix, we compute the late-time solution of the following coupled
PDEs

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

for a generic value of the parameter $\gamma$. Let us assume that two
symmetrical fronts, traveling at velocity $\pm 1$, develop at late times. We
work in the reference frame of the right-moving front by using

$\displaystyle\left\\{\begin{array}[]{rl}\phi_{+}(\tau,X)&\\!\\!\\!\\!:=\phi(\tau,X+\tau)\stackrel{{\scriptstyle\tau\to\infty}}{{\rightarrow}}f_{+}(X),\\\
\phi_{1+}(\tau,X)&\\!\\!\\!\\!:=\phi_{1}(\tau,X+\tau)\stackrel{{\scriptstyle\tau\to\infty}}{{\rightarrow}}f_{1+}(X).\end{array}\right.$
(183)

They obey the equations

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\phi_{+}+\partial_{X}(\phi_{1+}-\phi_{+})\\!\\!\\!\\!&=\phi_{+}(\phi_{+}^{2}-1),\\\
\partial_{\tau}\phi_{1+}-\partial_{X}(\phi_{1+}-\phi_{+})\\!\\!\\!\\!&=\phi_{1+}(\gamma\phi_{+}^{2}-1).\end{array}\right.$
(186)

As $\tau\to\infty$, this leads to

$\displaystyle\partial_{X}(f_{1+}-f_{+})=-f_{+}(1-f_{+}^{2})=f_{1+}(1-\gamma
f_{+}^{2}).$ (187)

Hence

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

and

$\displaystyle\partial_{X}\left(f_{+}\frac{1-f_{+}^{2}}{1-\gamma
f_{+}^{2}}+f_{+}\right)=f_{+}(1-f_{+}^{2}).$ (189)

This is a first-order ODE that can be integrated by means of decomposition
into simple fractions. Implicitly:

$\displaystyle-\frac{2}{\gamma
f_{+}^{2}-1}-\frac{(\gamma+1)\log(1-f_{+}^{2})}{\gamma-1}-\frac{(\gamma-3)\log(1-\gamma
f_{+}^{2})}{\gamma-1}+4\log f_{+}=2X+C.$ (190)

The constant $C$ can be determined by solving $f_{+}(R_{0}^{-})$ (i.e. the
discontinuity). We may extract the missing information from “the other side”,
i.e. on the left-moving front. Introduce

$\displaystyle\phi_{-}(\tau,X)=\phi(\tau,X-\tau),\qquad\phi_{1-}(\tau,X)=\phi_{1}(\tau,X-\tau).$
(191)

Using $\partial_{\tau}\phi_{-}=\partial_{\tau}\phi-\partial_{X}\phi$, we have

$\displaystyle\left\\{\begin{array}[]{rl}\partial_{\tau}\phi_{-}+\partial_{X}(\phi_{1-}+\phi_{-})\\!\\!\\!\\!&=-\phi_{-}(1-\phi_{-}^{2}),\\\
\partial_{\tau}\phi_{1-}+\partial_{X}(\phi_{1-}+\phi_{-})\\!\\!\\!\\!&=-\phi_{1-}(1-\gamma\phi_{-}^{2}).\end{array}\right.$
(194)

Then, integrating the first equation:

$\displaystyle\phi_{-}(\tau,X)$ $\displaystyle={\rm e}^{-\int_{0}^{\tau}{\rm
d}s\,(1-\phi_{-}(s,X)^{2})}\left[\phi_{0}(X)-\int_{0}^{\tau}{\rm d}u\,{\rm
e}^{\int_{0}^{u}{\rm
d}s\,(1-\phi_{-}(s,X)^{2})}(\partial_{X}\phi_{1-}+\partial_{X}\phi_{-})(u,X)\right]$
(195) $\displaystyle={\rm e}^{-\int_{0}^{\tau}{\rm
d}s\,(1-\phi_{-}(s,X)^{2})}\phi_{0}(X)-\int_{0}^{\tau}{\rm d}u\,{\rm
e}^{-\int_{u}^{\tau}{\rm
d}s\,(1-\phi_{-}(s,X)^{2})}(\partial_{X}\phi_{1-}+\partial_{X}\phi_{-})(u,X)$
(196) $\displaystyle={\rm e}^{-\int_{0}^{\tau}{\rm
d}s\,(1-\phi_{-}(\tau-s,X)^{2})}\phi_{0}(X)$
$\displaystyle\qquad-\int_{0}^{\tau}{\rm d}u\,{\rm e}^{-\int_{0}^{u}{\rm
d}s\,(1-\phi_{-}(\tau-s,X)^{2})}(\partial_{X}\phi_{1-}+\partial_{X}\phi_{-})(\tau-u,X).$
(197)

Add $2\tau$ to $X$, and use $\phi_{-}(\tau,X+2\tau)=\phi_{+}(\tau,X)$:

$\displaystyle\phi_{+}(\tau,X)$ $\displaystyle={\rm e}^{-\int_{0}^{\tau}{\rm
d}s\,(1-\phi_{+}(\tau-s,X+2s)^{2})}\phi_{0}(X+2\tau)$
$\displaystyle\qquad-\int_{0}^{\tau}{\rm d}u\,{\rm e}^{-\int_{0}^{u}{\rm
d}s\,(1-\phi_{+}(\tau-s,X+2s)^{2})}(\partial_{X}\phi_{1+}+\partial_{X}\phi_{+})(\tau-u,X+2u).$
(198)

When $X>R_{0}$, recalling that we postulate a solution such that
$\phi_{+}(\tau,X>R_{0})=1$, the above equation trivially simplifies to
$1=1-0$. When $X<R_{0}$, we take $\tau>\frac{R_{0}-X}{2}$. Then, since
$\phi_{+}(\tau,X>R_{0})=1$, we have:

$\displaystyle\phi_{+}(\tau,X)$ $\displaystyle={\rm
e}^{-\int_{0}^{\frac{R_{0}-X}{2}}{\rm d}s\,(1-\phi_{+}(\tau-s,X+2s)^{2})}$
$\displaystyle\qquad-\int_{0}^{\frac{R_{0}-X}{2}}{\rm d}u\,{\rm
e}^{-\int_{0}^{u}{\rm
d}s\,(1-\phi_{+}(\tau-s,X+2s)^{2})}(\partial_{X}\phi_{+1}+\partial_{X}\phi_{+})(\tau-u,X+2u).$
(199)

Before taking the limit $\tau\to\infty$, one is to be careful because the
spatial derivatives in the integrand above are unbounded as we expect the
discontinuity to develop: one cannot blindly replace $\phi_{+}$ by $f_{+}$.
Instead, write:

$\displaystyle(\partial_{X}\phi_{1+}+\partial_{X}\phi_{+})(\tau-u,X+2u)$
$\displaystyle=\frac{1}{2}\,\frac{{\rm d}}{{\rm
d}u}\Big{[}(\phi_{1+}+\phi_{+})(\tau-u,X+2u)\Big{]}$
$\displaystyle\qquad+\frac{1}{2}(\partial_{\tau}\phi_{1+}+\partial_{\tau}\phi_{+})(\tau-u,X+2u)$
$\displaystyle=\frac{1}{2}\,\frac{{\rm d}}{{\rm
d}u}\Big{[}(\phi_{1+}+\phi_{+})(\tau-u,X+2u)\Big{]}$
$\displaystyle\qquad-\frac{1}{2}\Big{[}\phi_{+}(1-\phi_{+}^{2})+\phi_{1+}(1-\gamma\phi_{+}^{2})\Big{]}(\tau-u,X+2u),$
(200)

using Eq. (186). Inserting Eq. (200) into Eq. (199) and integrate by parts:

$\displaystyle\phi_{+}(\tau,X)$ $\displaystyle={\rm
e}^{-\int_{0}^{\frac{R_{0}-X}{2}}{\rm
d}s\,(\cdots)}-\frac{1}{2}(\phi_{1+}+\phi_{+})\left(\tau-\frac{R_{0}-X}{2},R_{0}\right){\rm
e}^{-\int_{0}^{\frac{R_{0}-X}{2}}{\rm d}s(\cdots)}$
$\displaystyle\qquad\qquad+\frac{1}{2}(\phi_{1+}+\phi_{+})(\tau,X)+\int_{0}^{\frac{R_{0}-X}{2}}{\rm
d}u\,(\cdots).$ (201)

Crucially, the terms collected in $(\cdots)$ are bounded. Moreover,
$\phi_{1+}$ and $\phi_{+}$ are continuous and
$(\phi_{1+}+\phi_{+})(\tau,R_{0})=1$ for any finite time $\tau$. Therefore, we
can replace $(\phi_{1+}+\phi_{+})(\tau-\frac{R_{0}-X}{2},R_{0})$ by 1 all
times. Later sending $\tau\to\infty$, we obtain

$\displaystyle f_{+}(X)$ $\displaystyle={\rm
e}^{-\int_{0}^{\frac{R_{0}-X}{2}}{\rm d}s\,(\cdots)}-\frac{1}{2}{\rm
e}^{-\int_{0}^{\frac{R_{0}-X}{2}}{\rm
d}s(\cdots)}+\frac{1}{2}(f_{1+}+f_{+})(X)+\int_{0}^{\frac{R_{0}-X}{2}}{\rm
d}u\,(\cdots),$ (202)

and send $X\to R_{0}^{-}$

$\displaystyle f_{+}(R_{0}^{-})$
$\displaystyle=\frac{1}{2}+\frac{1}{2}(f_{1+}+f_{+})(R_{0}^{-})\qquad\mbox{{i.e.}}\qquad
f_{1+}(R_{0}^{-})=f_{+}(R_{0}^{-})-1.$ (203)

Let us call $L:=f_{+}(R_{0}^{-})$. Using Eq. (188), we have

$\displaystyle-L\frac{1-L^{2}}{1-\gamma L^{2}}=L-1.$ (204)

Simplifying by $L-1$

$\displaystyle(1+\gamma)L^{2}+L-1=0.$ (205)

The positive solution is

$\displaystyle L(\gamma)=\frac{\sqrt{5+4\gamma}-1}{2(1+\gamma)}.$ (206)

This quantity monotonically interpolates from $L=1/2$ for $\gamma=1$ to
$L=1/(\text{golden ratio})$ for $\gamma=0$.

## References

  * [1] E. H. Lieb and D. W. Robinson,  _The finite group velocity of quantum spin systems_ ,  Communications in Mathematical Physics 28(3), 251 (1972),  cmp/1103858407.
  * [2] C.-J. Lin and O. I. Motrunich,  _Out-of-time-ordered correlators in a quantum Ising chain_ ,  Physical Review B 97, 144304 (2018),  10.1103/PhysRevB.97.144304.
  * [3] S. Gopalakrishnan, D. A. Huse, V. Khemani and R. Vasseur,  _Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems_ ,  Physical Review B 98, 220303 (2018),  10.1103/PhysRevB.98.220303.
  * [4] S. Xu and B. Swingle,  _Accessing scrambling using matrix product operators_ ,  Nature Physics 16(2), 199 (2020),  10.1038/s41567-019-0712-4.
  * [5] A. Bohrdt, C. B. Mendl, M. Endres and M. Knap,  _Scrambling and thermalization in a diffusive quantum many-body system_ ,  New Journal of Physics 19(6), 063001 (2017),  10.1088/1367-2630/aa719b.
  * [6] A. A. Patel, D. Chowdhury, S. Sachdev and B. Swingle,  _Quantum butterfly effect in weakly interacting diffusive metals_ ,  Physical Review X 7, 031047 (2017),  10.1103/PhysRevX.7.031047.
  * [7] 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(5), 125 (2017),  10.1007/JHEP05(2017)125.
  * [8] D. Chowdhury and B. Swingle,  _Onset of many-body chaos in the $O(N)$ model_,  Physical Review D 96, 065005 (2017),  10.1103/PhysRevD.96.065005.
  * [9] S. H. Shenker and D. Stanford,  _Black holes and the butterfly effect_ ,  Journal of High Energy Physics 2014(3), 67 (2014),  10.1007/JHEP03(2014)067.
  * [10] D. A. Roberts and B. Swingle,  _Lieb-robinson bound and the butterfly effect in quantum field theories_ ,  Physical Review Letter 117, 091602 (2016),  10.1103/PhysRevLetter117.091602.
  * [11] A. Nahum, S. Vijay and J. Haah,  _Operator spreading in random unitary circuits_ ,  Physical Review X 8, 021014 (2018),  10.1103/PhysRevX.8.021014.
  * [12] C. W. von Keyserlingk, T. Rakovszky, F. Pollmann and S. L. Sondhi,  _Operator hydrodynamics, OTOCs, and entanglement growth in systems without conservation laws_ ,  Physical Review X 8, 021013 (2018),  10.1103/PhysRevX.8.021013.
  * [13] A. Chan, A. De Luca and J. T. Chalker,  _Solution of a minimal model for many-body quantum chaos_ ,  Physical Review X 8, 041019 (2018),  10.1103/PhysRevX.8.041019.
  * [14] B. Bertini and L. Piroli,  _Scrambling in random unitary circuits: Exact results_ ,  Physical Review B 102, 064305 (2020),  10.1103/PhysRevB.102.064305.
  * [15] 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),  https://doi.org/10.1016/j.aop.2016.09.006.
  * [16] X. Chen and T. Zhou,  _Quantum chaos dynamics in long-range power law interaction systems_ ,  Physical Review B 100, 064305 (2019),  10.1103/PhysRevB.100.064305.
  * [17] S. Xu and B. Swingle,  _Locality, quantum fluctuations, and scrambling_ ,  Physical Review X 9, 031048 (2019),  10.1103/PhysRevX.9.031048.
  * [18] S. Xu and B. Swingle,  _Scrambling dynamics and out-of-time ordered correlators in quantum many-body systems: a tutorial_ ,  10.48550/ARXIV.2202.07060 (2022).
  * [19] R. A. Fisher,  _The wave of advance of advantageous genes_ ,  Annals of Eugenics 7(4), 355 (1937),  https://doi.org/10.1111/j.1469-1809.1937.tb02153.x.
  * [20] 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(6), 1 (1937).
  * [21] E. Brunet and B. Derrida,  _An exactly solvable travelling wave equation in the Fisher–KPP class_ ,  Journal of Statistical Physics 161, 801 (2015).
  * [22] A. Larkin and Y. N. Ovchinnikov,  _Quasiclassical method in the theory of superconductivity_ ,  Journal of Experimental and Theoretical Physics 28, 960 (1969).
  * [23] M. Ansari and Y. V. Nazarov,  _Keldysh formalism for multiple parallel worlds_ ,  Journal of Experimental and Theoretical Physics 122(3), 389 (2016),  10.1134/s1063776116030134.
  * [24] 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),  10.1103/PhysRevB.98.205124.
  * [25] S. Chaudhuri and R. Loganayagam,  _Probing out-of-time-order correlators_ ,  Journal of High Energy Physics 2019(7), 6 (2019),  10.1007/JHEP07(2019)006.
  * [26] N. Kolganov and D. A. Trunin,  _Classical and quantum butterfly effect in nonlinear vector mechanics_ ,  Phys. Rev. D 106, 025003 (2022),  10.1103/PhysRevD.106.025003.
  * [27] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ ,  Cambridge University Press,  10.1017/CBO9781139003667 (2011).
  * [28] Y. Lemonik and A. Mitra,  _Time-resolved spectral density of interacting fermions following a quench to a superconducting critical point_ ,  Physical Review B 96(10) (2017),  10.1103/physrevb.96.104506.
  * [29] Y. Lemonik and A. Mitra,  _Model predictions for time-resolved transport measurements made near the superfluid critical points of cold atoms and $K_{3}C_{60}$ films_,  Physical Review Letters 121(6) (2018),  10.1103/physrevLetter121.067001.
  * [30] Y. Lemonik and A. Mitra,  _Quench dynamics of superconducting fluctuations and optical conductivity in a disordered system_ ,  Physical Review B 98(21) (2018),  10.1103/physrevb.98.214514.
  * [31] Y. Lemonik and A. Mitra,  _Transport and spectral signatures of transient fluctuating superfluids in the absence of long-range order_ ,  Physical Review B 100(9) (2019),  10.1103/physrevb.100.094503.
  * [32] A. I. Larkin and A. A. Varlamov,  _Fluctuation phenomena in superconductors_ ,  In K.-H.Bennemann and J. Ketterson, eds., _Handbook on Superconductivity: Conventional and Unconventional Superconductors_. Springer (2002).
  * [33] N. Ikeda, M. Nagasawa and S. Watanabe,  _Branching Markov processes II_ ,  Journal of Mathematics of Kyoto University 8(3), 365 (1968),  10.1215/kjm/1250524059.
  * [34] C. Aron, É. Brunet and A. Mitra,  _Kinetics of information scrambling in correlated electrons: disorder-driven transition from shock-wave to FKPP dynamics_ ,  10.48550/arXiv.2305.04958 (2023), 2305.04958.