






\section*{Abstract}
{\bf
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. They exhibit traveling wave solutions that are discontinuous at the boundary of the light cone, and have a perfectly causal structure where the solutions do not spill outside of the light cone.
}


\vspace{10pt}
\noindent\rule{\textwidth}{1pt}
\tableofcontents\thispagestyle{fancy}
\noindent\rule{\textwidth}{1pt}
\vspace{10pt}



\section{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, \textit{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~\cite{LiebRobinson1972}. 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~\cite{Motrunich2018,Vasseur2018}. Sharp wavefronts were found in interacting spin chains~\cite{Swingle2020}, non-integrable systems with diffusive transport~\cite{Knap2017,SwingleSachdev2017}, as well as large-$N$ or semiclassical models~\cite{Stanford2017,Swingle2017} and holographic models~\cite{Shenker2014,Swingle2016}. Interestingly, random unitary circuits without conserved quantities, \textit{i.e.} in the absence of diffusive transport, were found to develop broad fronts controlled by a diffusively growing length scale~\cite{Nahum2018, Sondhi2018,Chalker2018,Piroli2020}. 
Notably, several works~\cite{Aleiner2016,Zhou2019,Swingle2019,Swingle2022} 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~\cite{Fisher1937,KPP1937} which is known to exhibit traveling wave solutions (see Ref.~\cite{Brunet2015} 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 proximity to criticality produces a separation of energy 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.~\cite{Larkin1969} 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~\cite{Aleiner2016}. The augmented Keldysh formalism has also been adopted by others~\cite{Nazarov2016,Galitski2018,Chaudhuri2019,Trunin2022}.
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.

\begin{figure}
    \centering
    \includegraphics[scale=1.3]{sketch.eps}
    \caption{Butterfly effect: 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.~(\ref{eq:F_param}) follows a causal light-cone structure with a front traveling at a constant butterfly velocity $v_{\rm B}$ set by the Fermi velocity.
    Here, we sketch $\phi(\tau,{\boldsymbol{X}})$, the component of $F_{du}$ averaged over the Fermi surface, to be introduced in Eq.~(\ref{eq:Ansatz}).
    In the late-time limit, the traveling front develops a distinct discontinuity followed by an exponential decay on the scale of the mean free path $\ell$.
    }    \label{fig:sktech}
\end{figure}


\subsection*{Summary and main results}
The paper is organized as follows. In Sect.~\ref{sec:formalism}, 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.~\ref{sec:Superconducting_fluctuations}, 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
    \begin{align}  \nonumber
              \partial_t  F_{\alpha\beta}+ \boldsymbol{v}_{\boldsymbol{k}} \cdot \boldsymbol{\nabla}_{\boldsymbol{x}} F_{\alpha\beta}  &= I_{\alpha \beta}[F_{\alpha\beta},F_{\beta\alpha}],
    \end{align}
    with the non-linear collision integral $I_{\alpha \beta}$ given in Eq.~(\ref{eq:Coll_withoutP}).

In Sect.~\ref{sec:ansatz}, 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, \textit{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
    \begin{align}
    \left\{
\begin{array}{rl}  
     \partial_\tau \phi +  \boldsymbol{\nabla}_{\boldsymbol{X}} \cdot  \boldsymbol{\phi}_1 \!\!\!\!& =  \phi (\phi^2-1),
 \\[\medskipamount]
          \partial_\tau  \boldsymbol{\phi}_1 +  \boldsymbol{\nabla}_{\boldsymbol{X}} \phi  \!\!\!\!& = \boldsymbol{\phi}_1  (\gamma\phi^2-1),
\end{array}
\right. \nonumber
    \end{align}
where $\gamma = 1$ corresponds to the critical regime separating the superconducting phase from the normal phase. 
    

In Sect.~\ref{sec:solutions}, 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.~\ref{fig:sktech}.
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 initial exponential growth. Importantly, in the late-time regime, we find 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$), $f_+ = 1$ outside the light cone ($|{\boldsymbol{x}}| - v_{\rm B}t > 0$), and $f_+$ is discontinuous precisely at the boundary $|{\boldsymbol{x}}| - v_{\rm B}t =0$. We also work out explicitly the saturation regime, within the bulk of the light cone.

We conclude in Sect.~\ref{sec:discussion} by discussing the relations of our results to previous works and by giving future directions.



\section{Many-world formalism} \label{sec:formalism}
\begin{figure}
    \centering
    \includegraphics[scale=0.7]{contour.eps}
    \caption{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.~(\ref{eq:OTOC}).}
    \label{fig:manyworld}
\end{figure}

\subsection{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$)
\begin{align} \label{eq:OTOC}
\mathrm{Tr} \left[ \psi(0) \ {\rm e}^{{\rm i} H t} \psi({\boldsymbol{x}}) {\rm e}^{-{\rm i} H t} \ \psi^\dagger(0) \ {\rm e}^{{\rm i} H t'} \psi^\dagger({\boldsymbol{x}}') \ {\rm e}^{-{\rm i} H t'} \rho_0 \right],
\end{align}
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}}')$ 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 \newline 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.~\cite{Larkin1969} and we refer the reader to Ref.~\cite{Aleiner2016} for a detailed presentation which we simply follow.
The OTOC in Eq.~(\ref{eq:OTOC}) involves two forward and two backward time-evolution operators. Therefore, following the standard Schwinger-Keldysh construction~\cite{kamenev2011}, this yields the four-fold time-ordered contour $\mathcal{C}$ depicted in Fig.~\ref{fig:manyworld}. 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.~(\ref{eq:OTOC}) can be rewritten as
\begin{align} \label{eq:OTOC_2}
    \langle \mathrm{T}_\mathcal{C} \,    \psi_d^-(0,0)  \psi_{d}^{+}(t,{\boldsymbol{x}}) \psi_u^{-\dagger}(0,0)      \psi_{u}^{+\dagger}(t',{\boldsymbol{x}}')\rangle,
\end{align}
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$, \textit{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:
\begin{align} \label{eq:GF}
    {\rm i} G_{\alpha\beta}^{ab}(t,{\boldsymbol{x}};t',{\boldsymbol{x}}') &:= \langle \mathrm{T}_\mathcal{C} \psi_{\alpha}^a(t,{\boldsymbol{x}})  \psi_{\beta}^{b\dagger}(t',{\boldsymbol{x}}') \hat S_0 \rangle, 
\end{align}
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.~(\ref{eq:OTOC_2}) that are both evaluated at time $t=0$ but at disconnected locations on the time contour $\mathcal{C}$.
The OTOC in Eq.~(\ref{eq:OTOC}) 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',{\boldsymbol{x}}')$ 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, \textit{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.~(\ref{eq:GF}) for $\hat S_0 = 1$ shows that the intra-world Green's functions (\textit{i.e.} $\alpha = \beta$) correspond to the standard (\textit{i.e.} single-world) Schwinger-Keldysh two-point correlators: 
 \begin{align}
     G^{ab}_{uu} = G^{ab}_{dd} = G^{ab} \,.
 \end{align}

From now on, given that the intra-world quantities are identical in both $\alpha = u,d$ worlds, we simply drop the repeated world indices, \textit{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.

\subsection{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
\begin{align}
H &= H_0 +  H_{\rm int}\,, \label{eq:H_tot}\\
H_0 & = \sum_{{\boldsymbol{k}} \in {\rm BZ}}  \sum_{\sigma} \epsilon_{\boldsymbol{k}}  c^\dagger_{{\boldsymbol{k}}\sigma} c_{{\boldsymbol{k}}\sigma} \,,
\end{align}
where $H_0$ is the non-interacting part and the interaction in $H_{\rm int}$ depends on the specific problem at hand, see Eq.~(\ref{eq:BCS}).
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, \textit{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.

\subsection{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 (\textit{i.e.} $\alpha \neq \beta$) do not depend on the $\pm$ basis:
\begin{align}
  G^{++}_{\alpha\beta} &= G^{--}_{\alpha\beta} = G^{+-}_{\alpha\beta} = G^{-+}_{\alpha\beta} \ \mbox{ for } \alpha\neq \beta \,. \label{eq:G_R_ud}
\end{align}
On the other hand, the intra-world Green's functions (\textit{i.e.} $\alpha = \beta$) are related via~\cite{kamenev2011}
\begin{align}
    G_{\alpha\beta}^{+-} + G_{\alpha\beta}^{-+} = G_{\alpha\beta}^{++} + G_{\alpha\beta}^{--}  \ \mbox{ for } \alpha = \beta \,.
\end{align}
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:
\begin{align}
\left. \label{eq:Keldysh_basis}
\begin{array}{c}  
    G_{\alpha\beta}^R = G_{\alpha\beta}^{++} - G_{\alpha\beta}^{+- }\,, \qquad
     G_{\alpha\beta}^A = G_{\alpha\beta}^{++} - G_{\alpha\beta}^{-+} \,, \\[\medskipamount]
    G_{\alpha\beta}^K = G_{\alpha\beta}^{++}  + G_{\alpha\beta}^{--} \,,
   \end{array}
\right.
\end{align}
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
\begin{align}\left. \label{eq:Dyson1}
\begin{array}{rl}  
G^R(\omega,{\boldsymbol{k}}) &= \left[\omega - \epsilon_{\boldsymbol{k}} -\Sigma^R(\omega,{\boldsymbol{k}}) \right]^{-1} , \\[\medskipamount]
G^K(\omega,{\boldsymbol{k}}) & = 2{\rm i} \, F(\omega) \,\mathrm{Im}\, G^R(\omega,{\boldsymbol{k}}) \,.
\end{array}
\right.
\end{align}
$\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\newline $F(\omega):= \tanh\left(\omega/2T\right)$.

Concerning inter-world Green's functions, the relation in Eq.~(\ref{eq:G_R_ud}) together with Eq.~(\ref{eq:Keldysh_basis}) immediately implies
\begin{align}
  G_{ud}^R &= G_{du}^R =0 \,.
\end{align}
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
\begin{align}
\left. \label{eq:unstable_Dyson}
\begin{array}{rl}  
G_{ud}^K(\omega,{\boldsymbol{k}})\!\!\!\! &=  2{\rm i} [-1 + F(\omega)]\,\mathrm{Im}\, G^R(\omega,{\boldsymbol{k}}) \,, \\[\medskipamount]
G_{du}^K(\omega,{\boldsymbol{k}}) \!\!\!\! &=  2{\rm i} [+1 + F(\omega)]\,\mathrm{Im}\, G^R(\omega,{\boldsymbol{k}}) \,.
\end{array}
\right.
\end{align}
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
\begin{align}
    \textstyle G^K_{\alpha\beta}(t,{\boldsymbol{x}} ; t',{\boldsymbol{x}}') \!= \!\!\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', {\boldsymbol{x}}_2 - {\boldsymbol{x}}'),
    \label{eq:Dyson2}
\end{align}
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.

\subsection{Inter-world kinetic equation}
It is useful to work in the Wigner representation
\begin{align}
    G^K_{\alpha \beta}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}}) := \int\!\! {\rm d} {\boldsymbol{x}}' \!\! \int \!\!{\rm d} t' \,  {\rm e}^{{\rm i} (\omega t' - {\boldsymbol{k}} \cdot {\boldsymbol{x}}')} \,  G^K_{\alpha\beta}\left( t+\frac{t'}{2} ,{\boldsymbol{x}}+\frac{{\boldsymbol{x}}'}{2}; t-\frac{t'}{2},{\boldsymbol{x}}-\frac{{\boldsymbol{x}}'}{2}\right), 
\end{align}
and parameterize $G^K_{\alpha\beta}$ in terms of the real function $F_{\alpha\beta}$:
\begin{align}
    G^K_{\alpha \beta}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}})
    = F_{\alpha \beta}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}}) \star G^R(\omega,{\boldsymbol{k}}) - G^A(\omega,{\boldsymbol{k}}) \star F_{\alpha \beta}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}}),
    \label{eq:F_param}
\end{align}
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 approximation which consists in truncating the derivative expansion to its leading terms.

Massaging Eq.~(\ref{eq:F_param}) by acting on both sides with the inverse of $G^R$ and $G^A$, and using the Dyson equations~(\ref{eq:Dyson1}) and (\ref{eq:Dyson2}), one derives the kinetic equation on $F_{\alpha\beta}$ which is analogous to the standard (intra-world) kinetic equation:
\begin{align} \label{eq:inter_kinetic}
      \left[ \partial_t + \boldsymbol{v}_{\boldsymbol{k}} \cdot \boldsymbol{\nabla}_{\boldsymbol{x}} \right] F_{\alpha\beta}(\omega,{\boldsymbol{k}} ) &= I_{\alpha \beta}(\omega,{\boldsymbol{k}} )\,,
\end{align}      
where the {\sc 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 {\sc rhs} is the so-called collision integral which stems from $H_{\rm int}$ and reads
\begin{align}
     I_{\alpha \beta}(\omega,{\boldsymbol{k}} ) &= 2 \,  {\rm Im} \Sigma^{R}(\omega,{\boldsymbol{k}}) \,F_{\alpha \beta}(\omega,{\boldsymbol{k}})
    + {\rm i} \Sigma^K_{\alpha\beta}(\omega,{\boldsymbol{k}}) \,. \label{eq:Inter_Collision_Integral}
\end{align}
Let us recall that Keldysh components and collision integrals depend on space and time 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.~(\ref{eq:inter_kinetic}) is therefore a partial integrodifferential equation. It has a trivial steady-state solution
\begin{align}
    F_{ud}^ {\rm uncorr}(\omega,{\boldsymbol{k}})  = F_{du}^ {\rm uncorr}(\omega,{\boldsymbol{k}}) = 0,\label{eq:stable_SS}
\end{align}
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.~(\ref{eq:unstable_Dyson}), where both worlds evolve coherently, corresponds to another steady state characterized by
\begin{align}
\left. \label{eq:unstable_SS}
\begin{array}{l} 
F_{ud}^{\rm corr}(\omega,{\boldsymbol{k}}) =  [-1 + \tanh\left(\omega/2T\right)] \,,\\[\medskipamount]
F_{du}^{\rm corr}(\omega,{\boldsymbol{k}}) =  [+1 +  \tanh\left(\omega/2T\right)]  \,,
\end{array}
\right.
\end{align}
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 
\begin{align}
F_{\alpha\beta}(\omega,{\boldsymbol{k}};t=0,{\boldsymbol{x}}) &= [1 - \delta \phi_0({\boldsymbol{x}})]   F_{\alpha\beta}^{\rm corr}(\omega,{\boldsymbol{k}}),  \label{eq:cond_init}
\end{align}
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}$.

\section{Superconducting fluctuations}\label{sec:Superconducting_fluctuations}
\subsection{Model} 
\begin{figure}
    \centering
    \includegraphics[scale=0.45]{Untitled.eps}
    \caption{Superconducting self-energy $\Sigma$ and polarization bubble $\Pi$ of the model in Eq.~(\ref{eq:BCS}) treated in the RPA scheme in the particle-particle channel. The expressions in the Keldysh formalism are given in Eqs.~(\ref{eq:Sigma}) and (\ref{eq:Pi}), respectively.}
    \label{fig:diagrams}
\end{figure}
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.~(\ref{eq:H_tot}) reads
\begin{align} \label{eq:BCS}
    H_{\rm int} &= U \sum_{{\boldsymbol{k}}\bk'{\boldsymbol{k}}''}\sum_{\sigma} c^\dagger_{{\boldsymbol{k}}\sigma}c^\dagger_{-{\boldsymbol{k}}+{\boldsymbol{k}}'\bar\sigma}c_{{\boldsymbol{k}}''+{\boldsymbol{k}}'\bar\sigma}c_{-{\boldsymbol{k}}''\sigma} \,,
\end{align}
where $U<0$ is an attractive interaction facilitating superconductivity. 
In dimensions $d\geq2$, 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 \cite{Lemonik2017,Lemonik2018,Lemonik2018b,Lemonik2019}\begin{align}
\left.  \label{eq:D_R}
\begin{array}{l}  
D^R(\omega,{\boldsymbol{k}}) =   \left[ U^{-1} -  \Pi^R(\omega,{\boldsymbol{k}}) \right]^{-1} ,\\[\medskipamount]
D^K(\omega,{\boldsymbol{k}})  = 2{\rm i} \, P(\omega) \,\mathrm{Im}\, D^R(\omega,{\boldsymbol{k}}) \,,
\end{array}
\right.
\end{align}
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, \textit{i.e.} limiting ourselves to the one-loop diagrams depicted in Fig.~\ref{fig:diagrams} are:
\begin{align}
\left.  \label{eq:Pi}
\begin{array}{rl}  
\,\mathrm{Im}\, \Pi^R(\omega,{\boldsymbol{k}}) \!\!\!\!&=  \displaystyle \sum_{{\boldsymbol{k}}'} \int \frac{{\rm d} \omega'}{2\pi}  {\rm i} G^K(\omega',{\boldsymbol{k}}')  \,\mathrm{Im}\, G^R(\omega-\omega',{\boldsymbol{k}}-{\boldsymbol{k}}') \,, \\[\medskipamount]
\Pi^K(\omega,{\boldsymbol{k}}) \!\!\!\!&= 2{\rm i} \, P(\omega) \,\mathrm{Im}\,\Pi^R(\omega,{\boldsymbol{k}}) \,,
\end{array}
\right.
\end{align}
and
\begin{align}
\left.   \label{eq:Sigma}
\begin{array}{rl} 
\,\mathrm{Im}\,\Sigma^R(\omega,{\boldsymbol{k}}) \!\!\!\!&= \displaystyle - \frac{1}{2} \sum_{{\boldsymbol{k}}'} \int \frac{{\rm d} \omega'}{2\pi} \bigg{\{}
\,\mathrm{Im}\, D^R(\omega',{\boldsymbol{k}}') {\rm i}  G^K(\omega'-\omega,{\boldsymbol{k}}'-{\boldsymbol{k}}) \\
& \qquad\qquad\qquad\qquad -  {\rm i} D^K(\omega',{\boldsymbol{k}}')\,\mathrm{Im}\, G^R(\omega'-\omega,{\boldsymbol{k}}'-{\boldsymbol{k}}) \bigg{\}} \,, \\[\medskipamount]
\Sigma^K(\omega,{\boldsymbol{k}}) \!\!\!\!&=  2{\rm i} \, F(\omega) \,\mathrm{Im}\,\Sigma^R(\omega,{\boldsymbol{k}}) \,.
\end{array}
\right.
\end{align}
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~\cite{Lemonik2017,Lemonik2018,Lemonik2018b,Lemonik2019}.

Close to criticality, the Cooperon propagator reads~\cite{Larkin00,Lemonik2017,Lemonik2018,Lemonik2018b,Lemonik2019}
\begin{align}
    D^R(\omega,{\boldsymbol{k}}) \approx \frac{-1/\rho_{\rm F}}{r - {\rm i} a {\omega}/{T}+ \xi^2 k^2 + \ldots}, 
\end{align}
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\to0$, 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$.

\subsection{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.~(\ref{eq:Inter_Collision_Integral}). 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
\begin{align}
D_{\alpha\beta}^K(\omega,{\boldsymbol{k}}) &= \left| D^R(\omega,{\boldsymbol{k}}) \right|^2 \Pi_{\alpha\beta}^K(\omega,{\boldsymbol{k}}) \,,\\
\Pi_{\alpha\beta}^K(\omega,{\boldsymbol{k}}) &=  \frac{{\rm i}}{2} \sum_{{\boldsymbol{k}}'}  \int \frac{{\rm d} \omega'}{2\pi} G_{\alpha\beta}^K(\omega',{\boldsymbol{k}}') G_{\alpha\beta}^K(\omega-\omega',{\boldsymbol{k}}-{\boldsymbol{k}}')\,, \\
\Sigma_{\alpha\beta}^K(\omega,{\boldsymbol{k}}) &=  -\frac{{\rm i}}{2} \sum_{{\boldsymbol{k}}'}  \int \frac{{\rm d} \omega'}{2\pi} D_{\alpha\beta}^K(\omega',{\boldsymbol{k}}')  G_{\beta\alpha}^K(\omega'-\omega,{\boldsymbol{k}}'-{\boldsymbol{k}}) \,.
\end{align}
Altogether, this yields the inter-world kinetic equation~(\ref{eq:inter_kinetic}) with the collision integral
\begin{align}
     I_{\alpha \beta}(\omega, {\boldsymbol{k}}) 
     & =  2
  \sum_{{\boldsymbol{k}}'{\boldsymbol{k}}''} \int \frac{{\rm d} \omega'}{2\pi} \frac{{\rm d} \omega''}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2 
  \nonumber \\
  & \qquad \times   
  \,\mathrm{Im}\, G^R(\omega'',{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega,{\boldsymbol{k}}'-{\boldsymbol{k}}) \nonumber \\
  & \qquad \times \bigg{\{}  \left[ \tanh\left(\frac{\omega''}{2T}\right) +\tanh\left(\frac{\omega'-\omega''}{2T}\right) \right] \left[ \coth\left(\frac{\omega'}{2T}\right) +\tanh\left(\frac{\omega-\omega'}{2T}\right) \right]  F_{\alpha\beta}(\omega,{\boldsymbol{k}}) \nonumber \\ 
  &\qquad\qquad  + 
   F_{\alpha\beta}(\omega'',{\boldsymbol{k}}'')
   F_{\alpha\beta}(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'')
   F_{\beta\alpha}(\omega'-\omega,{\boldsymbol{k}}'-{\boldsymbol{k}}) \bigg{\}} \,. \label{eq:Coll_withoutP} 
  \end{align}
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.~(\ref{eq:stable_SS}) and the correlated-world solution in Eq.~(\ref{eq:unstable_SS}).

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 \newline $\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.

\section{Inter-world kinetics in the near-critical regime}
\label{sec:ansatz}
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~(\ref{eq:inter_kinetic}) with its collision integral in Eq.~(\ref{eq:Coll_withoutP}).
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.~\ref{sec:solutions}.

\subsection{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 \textit{e.g.} Eq.~(\ref{eq:Coll_withoutP}). 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 (\textit{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.~(\ref{eq:inter_kinetic}) with the following ansatz:
\begin{align} \label{eq:Ansatz}
    \tilde F^{\rm ansatz}_{\alpha\beta}(t,{\boldsymbol{x}};{\boldsymbol{k}}) &= \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}}),
\end{align}
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~\cite{kamenev2011} of $\tilde{F}_{\alpha\beta}$, accounting for its isotropic and anisotropic contributions in momentum space:
\begin{align} \label{eq:phi_phi1}
 \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}}, \\[\medskipamount]
     \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.
\end{align}
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, \textit{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.~(\ref{eq:Ansatz}).
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.~(\ref{eq:Ansatz}) and why the collision integral cannot be, \textit{a priori}, linearized. 
In that regard, it is similar to the ansatz used in Ref.~\cite{Aleiner2016}.

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 
\begin{align}
    \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.
\end{align}
to another momentum-space isotropic and real-space homogeneous (stable) solution
\begin{align}
    \tilde F^{\rm uncorr}({\boldsymbol{k}};t=0,{\boldsymbol{x}}) \longleftrightarrow \left\{ \begin{array}{rl}
      \phi(t=0,{\boldsymbol{x}})   &= 0,  \\
    \boldsymbol{\phi}_1(t=0,{\boldsymbol{x}})  & = 0 \,.
    \end{array}
    \right.
\end{align}
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.~(\ref{eq:cond_init}). 
\subsection{Simplified kinetic equation: coupled PDEs}
We proceed by injecting the ansatz (\ref{eq:Ansatz}) in the inter-world collision integral in Eq.~(\ref{eq:Coll_withoutP}) 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.~(\ref{eq:Coll_withoutP}), this means that the term $|D^R(\omega',{\boldsymbol{k}}')|^2$ diverges as $\omega' \approx 0$ and $k'\approx 0$.
The details of the computation are given in Appendix~\ref{app:HT_proj}. We obtain
\begin{align}
    \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}}) \,.
\end{align}
Injecting the ansatz in the kinetic equation~(\ref{eq:inter_kinetic}), dividing by $\tilde F_{\alpha\beta}^{\rm corr}({\boldsymbol{k}})$, we get
\begin{align}
     &\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) \nonumber \\
     &\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}}) \,.
\end{align}
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, \textit{i.e.} eventually setting $k \to k_{\rm F}$, we obtain the following set of coupled partial differential equations (PDEs)
\begin{align} \label{eq:coupled_PDEs}
\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},  \\[\medskipamount]
          \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.
\end{align}
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}}$)
\begin{align}
\frac{1}{\tau_{\rm F}} := -2 \,\mathrm{Im}\, \tilde \Sigma^R(k_{\rm F}) \,. \label{eq:def_tds}
 \end{align}
$\gamma$ is a dimensionless parameter that generalizes the computation performed at criticality, for which $\gamma=1$, to near critical regimes for which $\gamma\neq1$ (see Appendix~\ref{app:away_from_crit}).
After an appropriate rescaling of space and time, 
\begin{align} \label{eq:rescaling}
 \tau := t / \tau_{\rm F} \ \mbox{ and } \ {\boldsymbol{X}} := {\boldsymbol{x}} \sqrt{d} /(v_{\rm F} \tau_{\rm F})\,,
\end{align} together with $\boldsymbol{\phi}_1 \mapsto \sqrt{d} \boldsymbol{\phi}_1$, the coupled PDEs now only involve dimensionless quantities
\begin{align}
\left\{     \label{eq:coupled_PDEs_adim}
\begin{array}{rl}  
     \partial_\tau \phi +  \boldsymbol{\nabla}_{\boldsymbol{X}} \cdot  \boldsymbol{\phi}_1 \!\!\!\!& =  \phi (\phi^2-1),
 \\[\medskipamount]
          \partial_\tau  \boldsymbol{\phi}_1 +  \boldsymbol{\nabla}_{\boldsymbol{X}} \phi  \!\!\!\!& = \boldsymbol{\phi}_1  (\gamma\phi^2-1) \,.
\end{array}
\right.
\end{align}
Importantly, the generic $(d+1)$-dimensional case can be reduced to an effective $(1+1)$-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
\begin{align}
\left\{ \label{eq:coupled_PDEs^R}
\begin{array}{rl} 
     \partial_\tau \phi +  \partial_r \phi_1  \!\!\!\!& =  - \frac{d-1}{r} \phi_1 + \phi (\phi^2-1),
     \\[\medskipamount]
          \partial_\tau  \phi_1 +  \partial_r \phi  \!\!\!\!& = \phi_1  (\gamma\phi^2-1)  \,.
          \end{array}
\right.
\end{align}

The set of PDEs in (\ref{eq:coupled_PDEs}) and the following expressions in Eqs.~(\ref{eq:coupled_PDEs_adim}),~(\ref{eq:coupled_PDEs^R}) are one of the main results 
of this manuscript. 
Overall, this represents a considerable simplification from the original partial integrodifferential kinetic equation (\ref{eq:inter_kinetic}) governing the dynamics of the inter-world distribution function $F_{\alpha\beta}(\omega,{\boldsymbol{k}};t,{\boldsymbol{x}})$ with the collision integral in Eq.~(\ref{eq:Coll_withoutP}).  

To provide a first intuitive understanding of the previous set of PDEs, let us briefly neglect spatial inhomogeneities of $\phi$ and $\boldsymbol{\phi}_1$. The equation on $\phi(t)$  becomes an `autonomous' first-order ODE, reading
\begin{align}
     \partial_\tau \phi & = \, \phi (\phi^2-1) = - V'(\phi) \,. \label{eq:gradient_descent}
\end{align}
This is a gradient descent in the potential $V(\phi) = \frac12 \phi^2 - \frac14 \phi^4 $. 
Reinstating the original units, the logarithmic 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_* \sim - \frac12 \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, 
\begin{align}
    1-\phi(t\ll t_*) \sim \exp\left[2(t-t_*)/\tau_{\rm F}\right]\,,
\end{align}
while it saturates at late times,
\begin{align}
    \phi(t \gg t_*) \sim \exp\left[ -(t-t_*)/\tau_{\rm F} \right] \to 0\,.
\end{align}
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.~\ref{sec:solutions}.

In the general case, \textit{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 (\textit{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 (\ref{eq:coupled_PDEs}) governing the dynamics of $\phi$ and $\boldsymbol{\phi}_1$. Consequently, one cannot \textit{a priori} apply the standard hydrodynamic approach, and one must solve for the dynamics $\phi$ and $\boldsymbol{\phi}_1$ on an equal footing.

\subsection{Validating the ansatz numerically} 
We perform two independent checks of the ansatz proposed in Eq.~(\ref{eq:Ansatz}) by comparing, on the one hand, the solutions of the inter-world kinetic equation in~(\ref{eq:inter_kinetic}) computed with the full-fledged collision integral in Eq.~(\ref{eq:Coll_withoutP}) with, on the other hand, the solutions of the coupled PDEs in~(\ref{eq:coupled_PDEs}).

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~\ref{sec:Superconducting_fluctuations} 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
\begin{align} \label{eq:G_R_quasiP}
    G^R(\omega,k) = \frac{1}{\omega-\epsilon_k + {\rm i} \Gamma} \,,
\end{align}
and where the Cooperon Green's function $D^R(\omega,k)$ is computed following Eq.~(\ref{eq:D_R}). $\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)$. 

\subsubsection{Early times} \label{sec:ansatz_early}
\begin{figure}
\centering
\hspace{-.5cm}
\begin{subfigure}{.5\linewidth}
  \centering
  \input{kineticEQ_kx.tex}
\end{subfigure}
\hspace{.2cm}
\begin{subfigure}{.5\linewidth}
  \centering
  \input{kineticPDE_1d.tex}
\end{subfigure}
\caption{(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.~(\ref{eq:cond_init_Gauss}) with $\lambda = 0.05$ and $\delta\phi_0=0.98$. (b) Solution $\phi(\tau,X)$ to the coupled PDEs in (\ref{eq:coupled_PDEs}) 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$.
    Time and space have been rescaled according to Eq.~(\ref{eq:rescaling}).
    No adjustable parameters were used.
   }
\label{fig:full-fledged_num}
\end{figure}

The first test consists in numerically solving the kinetic equation (\ref{eq:inter_kinetic}) with the full-fledged collision integral in Eq.~(\ref{eq:Coll_withoutP}) for as long as we can ensure numerical convergence of the solutions. In practice, this is challenging and we can only access early times, \textit{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 $\lambda$, and localized around $X=0$. Explicitly, we take the following symmetric initial condition 
\begin{subequations}\label{eq:cond_init_Gauss}
\begin{align}  
& \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. \\
& \mbox{ with } 
\delta\phi_0 (X) = \delta\phi_0 \exp[ - X^2/(2\lambda^2) ] \Theta(3\lambda-|X|)\,,
\label{eq:droplet}
\end{align}
\end{subequations}
and where $\Theta(X)$ is the Heaviside step function. We found such a Gaussian-shaped droplet, defined on the support $[-3\lambda,3\lambda]$, 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.~\ref{fig:full-fledged_num}, 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.~(\ref{eq:phi_phi1}) which, in 1D, simply reads
\begin{align}
\left\{ \label{eq:extraction}
\begin{array}{rl} 
    \textstyle \phi(\tau,X) \!\!\!\!&= \textstyle \frac12 \left[ \tilde{F}_{du}(\tau,X;k_{\rm F})  + \tilde{F}_{du}(\tau,X;-k_{\rm F})\right], 
    \\[\medskipamount]
        \textstyle \phi_1(\tau,X) \!\!\!\! &= \textstyle \frac12 \left[ \tilde{F}_{du}(\tau,X;k_{\rm F})  - \tilde{F}_{du}(\tau,X;-k_{\rm F})\right]   \,,
        \end{array}
\right.
\end{align}
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, typical of wave-like solutions (as opposed to diffusive solutions).
We repeated this analysis in a wide range of 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.~(\ref{eq:Ansatz}) at early times.

\subsubsection{Late times and partial-wave truncation}
\begin{figure}
\centering
\hspace{-.5cm}
\begin{subfigure}{.5\linewidth}
  \centering
      \input{ansatz.tex}
\end{subfigure}%
\hspace{.2cm}
\begin{subfigure}{.5\linewidth}
  \centering
  \input{ansatz_phi1.tex}
\end{subfigure}
\caption{
        (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.~(\ref{eq:cond_init_F2}) with $\phi_0=0.9$ and ${\phi_{1}}_0 = 0.1$, and compared with the solution $\phi(\tau)$ to the coupled PDEs in Eq.~(\ref{eq:coupled_PDEs_0d}).
        (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.~(\ref{eq:rescaling}).
    }
    \label{fig:ansatz}
\end{figure}

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
\begin{align} \label{eq:cond_init_F2}
\tilde F_{ab}(\tau=0;k) = \left[ \phi_0 + {\phi_{1}}_0 {\rm sign}(k) \right] \tilde{F}^{\rm corr}_{ab}(k),
\end{align}
on the kinetic equation side and, correspondingly,
\begin{align} 
\phi(\tau = 0) = \phi_0 \ \mbox{ and } \ \phi_1(\tau = 0) = {\phi_1}_0,
\end{align}
on the side of the coupled PDEs which now read
\begin{align}
\left\{ \label{eq:coupled_PDEs_0d}
\begin{array}{rl}  
     \partial_\tau \phi \!\!\!\!& =  \phi (\phi^2-1),
     \\
          \partial_\tau  {\phi}_1   \!\!\!\!& = {\phi}_1  (\gamma\phi^2-1) \,.
          \end{array}
\right.
\end{align}
Note that these coupled ordinary-differential equations correspond to the discussion around Eq.~(\ref{eq:gradient_descent}). Notably, the Eq.~(\ref{eq:coupled_PDEs_0d}) predicts that the relaxation dynamics of $\phi$ (but not those of $\phi_1$) are independent  of the distance to criticality, parameterized by $\gamma$. 
Therefore, and complementary to the previous benchmark in Sect.~\ref{sec:ansatz_early}, 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.~\ref{fig:ansatz}, 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.~(\ref{eq:extraction}).
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.

\section{Dynamics of information scrambling}
\label{sec:solutions}

In this Section, we solve the set of coupled PDEs~(\ref{eq:coupled_PDEs_adim}) 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.

\subsection{Early-time exponential growth}  \label{sec:solutions_early}
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 initial condition
\begin{align}
\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.
\end{align}
where the perturbation $0 < \delta \phi_0({\boldsymbol{X}}) \ll 1 $ is non-vanishing on a small compact support of linear size $\lambda$ around ${\boldsymbol{X}}=0$. Note that this is a different regime from the numerics presented in Sect.~\ref{sec:ansatz_early} 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
\begin{align}
\left\{ \label{eq:coupled_PDEs_linearized}
\begin{array}{rl}  
     \partial_\tau \delta \phi - \partial_X \phi_1 \!\!\!\!& =  2 \delta  \phi,
     \\
    \partial_\tau  \phi_1 -  \partial_X \delta \phi  \!\!\!\!& = 0 \,.
              \end{array}
\right.
\end{align}
Note that the dependence on $\gamma$, the parameter quantifying the distance to criticality, has dropped. This yields the following PDE on $\delta \phi(\tau,X)$
  \begin{align} \label{eq:PDE_growth}
      \partial_\tau^2 \delta\phi  - 2 \partial_\tau \delta\phi -  \partial^2_X \delta \phi & =  0\,.
  \end{align}
Integrating the above equation once over the whole space and introducing the integrated perturbation $M(\tau) := \int{\rm d} X \delta\phi(\tau,X)$, we find an exponential growth of the perturbation, echoing the onset of chaos:
\begin{align}
    M(\tau \ll \tau_*) = \exp\left[ 2(\tau - \tau_*) \right],
\end{align}
with the typical timescale to escape from the unstable solution given by $\tau_* \sim -\frac12 \log \delta M_0$ where $\delta M_0 := \int {\rm d} X \delta \phi_0(X)$.
A more sophisticated calculation restricted to $\lambda \ll X \ll \tau \ll \tau_*$  yields
\begin{align}\label{eq:diffu}
\delta \phi(\tau,X)  = {\rm e}^{2\tau} \int \!{\rm d} Y \,
\delta\phi_0(Y)
\frac{{\rm e}^{-(X-Y)^2/2\tau}}{\sqrt{2\pi \tau}}
 \approx {\rm e}^{2(\tau-\tau_*)} \frac{{\rm e}^{-X^2/2\tau}}{\sqrt{2\pi \tau}}\,.
\end{align}
One may check, by direct substitution, that the above expression is indeed a solution to the PDE~(\ref{eq:PDE_growth}), up to corrections on the order of ${\rm e}^{2(\tau-\tau_*)}/\tau^{5/2}$ which are vanishingly small in the linear regime.
The solution in~(\ref{eq:diffu}) involves a diffusion kernel and an exponential growth. It is similar to the one obtained in the context of the $\mathcal{O}(N)$ model in Ref.~\cite{Swingle2017} (see Eq.~1.13 therein). It can be seen as the solution to the linearly-driven diffusion equation 
\begin{align}
    \partial_\tau \delta \phi - \frac12 \partial_X^2 \delta \phi = 2 \delta \phi,
\end{align}
that appears, notably, in the context of branching Brownian motion~\cite{BBM1968}.

\subsection{Late-time solutions and discontinuous front}
Here, we access the late-time inter-world distribution by analytically solving the coupled PDEs (\ref{eq:coupled_PDEs_adim}) 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.

\subsubsection{Traveling front}
\label{sec:propagating_front}

\begin{figure}
\centering
\hspace{-.5cm}
\begin{subfigure}{.5\linewidth}
  \centering
\input{Fronta.tex}
\end{subfigure}%
\hspace{.2cm}
\begin{subfigure}{.5\linewidth}
  \centering
  \input{Frontb.tex}
\end{subfigure}
\caption{(a) Solution $\phi(\tau,X)$ to the $d=1$ coupled PDEs in Eq.~(\ref{eq:coupled_PDEs_adim}) at $\gamma=1$ for the droplet initial perturbation in Eq.~(\ref{eq:droplet}) with amplitude $\delta\phi_0 = 0.02$ and linear size $\lambda = 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 + \lambda$.  No adjustable parameter was used.}
    \label{fig:front}
\end{figure}

\begin{figure}
    \centering
    \input{Frontc.tex}
    \caption{
    Discontinuous late-time front:
     the solution $\phi(\tau,X-X_{\rm front}(\tau))$ is extracted from the data of Fig.~\ref{fig:front}~(a) for $\tau = 28$ (solid line), and compared to the analytical prediction $f_+(X)$ 
 in Eq.~(\ref{eq:front_shape}) for the steady-state front, with a discontinuity at $X=0$ from 0.5 to 1 (dashed line). No adjustable parameter was used.}
    \label{fig:frontc}
\end{figure}

We assume that the solution develops a traveling front located on a sphere of increasing radius. Our goal is to compute its velocity and steady-state 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.~(\ref{eq:coupled_PDEs^R}). 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 \neq1$ can also be solved by similar techniques and we refer the reader to Appendix~\ref{app:alpha_neq_1} for the corresponding detailed computation.
Introducing the fields
\begin{align}
    \varphi := \frac{\phi+\phi_1}{2} \quad \mbox{ and }  \quad \psi := \phi - \phi_1 \,,
\end{align}
the coupled PDEs can be cast as
\begin{align} \label{eq:PDE_coupled}
\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.
\end{align}
with the initial conditions,
\begin{align}
\left\{  
\begin{array}{l} 
    \varphi(\tau=0,X) = \phi_0(X)/2, \\
    \psi(\tau = 0,X) = \phi_0(X)\,.
\end{array}
\right.
\end{align}
The method of characteristics gives the implicit solutions
\begin{align}
\left\{  
\begin{array}{l}  
    \varphi(\tau,X) = \frac12\phi_0(X-\tau)\rm e^{-\int_0^\tau {\rm d} s \left[ 1-\phi^2(X-\tau+s,s) \right]}, \\
    \psi(\tau,X) = \phi_0(X+\tau) {\rm e}^{-\int_0^\tau {\rm d} s \left[ 1-\phi^2(X+\tau-s,s) \right]}, 
    \end{array}
\right.
\end{align}
which yield
\begin{align}
    \phi(\tau,X) &=  \frac12\phi_0(X-\tau) {\rm e}^{-\int_0^\tau {\rm d} s \left[ 1-\phi^2(s,X-\tau+s) \right]} + \frac12  \phi_0(X+\tau) {\rm e}^{-\int_0^\tau {\rm d} s \left[ 1-\phi^2(s,X+\tau-s) \right]}   \,.
\end{align}
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, \textit{i.e.}
\begin{align}
    \phi_+(\tau,X) := \phi(\tau,X+\tau) \stackrel{\tau \to\infty}{\to}  f_+(X)\,,
\end{align}
pointwise. Naturally, there is also a symmetrical left-moving front.
We start from
\begin{align}
    \phi_+(\tau,X) & = \frac12 \phi_0(X) {\rm e}^{-\int_0^\tau {\rm d} s \left[ 1-\phi_+^2(s,X) \right]} + \frac12  \phi_0(X+2\tau) {\rm e}^{-\int_0^\tau {\rm d} s \left[ 1-\phi_+^2(s,X+2\tau-2s) \right]}\\
    & = \frac12\phi_0(X) {\rm e}^{-\int_0^\tau {\rm d} s \left[ 1-\phi_+^2(s,X) \right]} + \frac12  \phi_0(X+2\tau) {\rm e}^{-\int_0^\tau {\rm d} s \left[ 1-\phi_+^2(\tau - s,X+2 s) \right]}  \,,
\end{align}
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
\begin{align}
    f_+(X) & = \frac12\phi_0(X) {\rm e}^{-\int_0^\infty {\rm d} s \left[ 1-\phi_+^2(s,X) \right]}  + \frac12   {\rm e}^{-\int_0^\infty {\rm d} s \left[ 1-f_+^2(X+2s) \right]} \,,   \label{eq:fplus_intermed}
\end{align}
where we used the property $\phi_0(X+2\tau) \to 1$.
We postulate\footnote{This assumption can be rigorously proven to be true.} that the front is such that 
\begin{align}
   \phi_+(\tau, X < \lambda ) < 1 \quad \mbox{ and} \quad \phi_+(\tau, X > \lambda ) = 1,  
\end{align}
at all times, and therefore
\begin{align}
  f_+(X < \lambda )  < 1  \quad \mbox{ and}  \quad  f_+(X > \lambda ) = 1    \,.
\end{align}
Let us now work with $X <\lambda$.
Hence, $ 1- \phi_+^2(X < \lambda) > 0$, which nullifies the first term in Eq.~(\ref{eq:fplus_intermed}). We are left with
\begin{align}
    f_+( X < \lambda) & =  \frac12   {\rm e}^{-\int_0^\infty {\rm d} s \left[ 1-f_+^2(X+2s) \right]}  = \frac12   {\rm e}^{-\frac12 \int_X^\lambda {\rm d} u \left[ 1-f_+^2(u) \right]} \,. \label{eq:fplus}
\end{align}
At $X \stackrel{X<\lambda}{\longrightarrow} \lambda$, this yields $f_+(\lambda) = 1/2$.
Given that we assume $f_+(X > \lambda) = 1$, this signals a discontinuity in $f_+(X)$ at $X = \lambda$.
Moreover, Eq.~(\ref{eq:fplus}) implies that $f_+(X)$ obeys
\begin{align}
    f_+'(X) & =  \frac12 \left[1-f_+^2(X) \right] f_+(X),
\end{align}
which can be solved by separation of variables, yielding the discontinuous right-moving front with the shape
\begin{align}
   f_+(X < \lambda) =    \frac{1}{\sqrt{1+ 3 {\rm e}^{\lambda-X}}} 
    \stackrel{X \to \lambda^-}{\mapsto} \frac12 
    \quad \mbox{ and}  \quad 
    f_+(X > \lambda) = 1\,.
\end{align}
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, $\lambda$ drops out of the problem.
This is expected for generic initial conditions defined on a compact support $|X| < \lambda$. Indeed, one can show that the first non-vanishing derivative $\partial^{(n)}_X \phi_0(X=\lambda^-) > 0$ is responsible for the generation of a first-order derivative $f_+(X = \lambda^-) > 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.
Noteworthy, as previously discussed in Sec.~\ref{sec:solutions_early}, $\lambda$ and the precise shape of the initial perturbation are however controlling the timescale for the steady-state front to form. Barring this point, the universal features of the steady state may be safely accessed by sending $\lambda \to 0^+$ after $\tau\to\infty$. We obtain the right-moving steady-state front
\begin{align}
\left\{ \label{eq:front_solution}
\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.
\end{align}
with
\begin{align} \label{eq:front_shape}
    \left\{
    \begin{array}{cc}
    f_+(X < 0) = \frac{1}{\sqrt{1+ 3{\rm e}^{-X}}} \stackrel{X \to 0^-}{\mapsto} \frac12 \quad \mbox{ and}  \quad  f_+(X > 0) = 1, \\
    f_{1+}(X < 0) = -\frac{1}{\sqrt{1+ 3{\rm e}^{-X}}} \stackrel{X \to 0^-}{\mapsto} - \frac12 \quad \mbox{ and}  \quad  f_{1+}(X > 0) = 0\,.
    \end{array}
    \right.
\end{align}
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
\begin{align} \label{eq:late_time_solution}
  \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} 
\end{align}
with the  butterfly velocity given by $v_{\rm B}:=  {v_{\rm F}}/{\sqrt{d}}$ and the radial front shape $f_+$ given by Eq.~(\ref{eq:front_shape}) 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.~(\ref{eq:def_tds}).

This solution can also be generalized to cases away from criticality, \textit{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~\ref{app:alpha_neq_1}, we show that the discontinuity of the front is now from $f_+(X=0^-) = L(\gamma)$ to $f_+(X > 0) = 1$ with
\begin{align} \label{eq:L_alpha}
L(\gamma)=\frac{\sqrt{5+4\gamma}-1}{2(1+\gamma)} \,.
\end{align}
This quantity interpolates monotonically from $L=1$ (no discontinuity) for $\gamma=-1$ to $L=1/2$ for $\gamma=1$, going through $L=1/(\text{golden ratio})$ for $\gamma=0$. 

The discontinuous traveling wavefronts in Eqs.~(\ref{eq:front_shape}), (\ref{eq:late_time_solution}), and their generalization to non-critical regimes in Eq.~(\ref{eq:L_alpha}) 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.~\ref{fig:front}, we display the numerical solution starting from the following droplet-shaped initial perturbation of the correlated-world solution
\begin{align} 
& 
\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/\lambda)^2} \Theta(\lambda-|X|), 
     \\
          {\phi}_1(\tau=0,X)   \!\!\!\!& = 0  \,,
          \end{array}
\right. 
\label{eq:droplet}
\end{align}
where $\delta\phi_0$ sets the amplitude of the droplet, and $\lambda$ 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.~\ref{fig:frontc}, we illustrate the discontinuous shape of the steady-state 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.~(\ref{eq:front_shape}).

\subsubsection{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 {\sc rhs} of the coupled PDEs~(\ref{eq:PDE_coupled}). Working in the $d=1$ case, the linearized PDEs read
\begin{align}
\left\{ 
\begin{array}{rl} 
     \partial_\tau \phi +  \partial_X  \phi_1 \!\!\!\! & =  -\phi,
     \\
          \partial_\tau  \phi_1 + \partial_X \phi  \!\!\!\!& = -\phi_1\,.
          \end{array}
\right.
\end{align}
Assuming symmetric initial conditions, \textit{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
\begin{align}
\left\{  \label{eq:inside_light_cone}
\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.
\end{align}
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:
\begin{align}
    \lim\limits_{\tau\to\infty} \phi(\tau,X+\tau) \sim f_+(X) \mbox{  for  } X < 0\, , |X| \gg 1,
\end{align}
where $f_+(X)$ was computed in Eq.~(\ref{eq:front_shape}). Using the late-time solution in Eq.~(\ref{eq:inside_light_cone}), we have
\begin{align}
    \lim\limits_{\tau\to\infty} {\rm e}^{-X-\tau} f_\varphi(X) + {\rm e}^{X+\tau} f_\varphi(-X-2\tau)  \sim f_+(X) \,.
\end{align}
The first term vanishes and we are left with
\begin{align}
    \lim\limits_{\tau\to\infty} {\rm e}^{X+\tau} f_\varphi(-X-2\tau)  \sim f_+(X)\,,
\end{align}
which is solved as 
\begin{align}
    f_\varphi(X) &\sim  {\rm e}^{X/2}\,.
\end{align}
This yields, at large times and far inside from the boundary of the light cone,
\begin{align}
\left\{   \label{eq:phi_cosh}
\begin{array}{rl}  
\phi(\tau,X < \tau) \!\!\!\!& \sim  2{\rm e}^{-\tau/2} \cosh{(X/2)}, \\
\phi_1(\tau,X < \tau ) \!\!\!\! & \sim -2 {\rm e}^{-\tau/2} \sinh{(X/2)}\,.
\end{array}
\right.
\end{align}
Note that these asymptotic solutions have already lost the information about the initial condition. Moreover, they imply the relation $\phi_1(\tau,X < \tau)  = - 2 \partial_X \phi(\tau,X < \tau)$. Once re-injected in the linearized PDEs, this yields the following linearly-driven diffusion equation for $\phi(\tau,X)$ only
\begin{align} \label{eq:PDE_linear_diffusive}
    \partial_\tau \phi - 2 \partial_X^2 \phi = -\phi \,,
\end{align}
where the diffusive constant in front of the $\partial_X^2 \phi$ term is $D := 2 v_{\rm F}^2 \tau_F$ once the original units of time and space are reinstated.
We note that the diffusive {\sc lhs} of Eq.~(\ref{eq:PDE_linear_diffusive}) is similar to the one of the FKPP-like equation that was derived in Ref.~\cite{Aleiner2016} in a similar context.

\section{Discussion and conclusion} \label{sec:discussion}
Starting from the microscopic Hamiltonian of a $(d+1)$-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.

What is perhaps the most striking result of our work is the fact that our traveling front solution scrupulously respects causality: it does not leak outside of the light cone. Relatedly, our solution presents a distinct discontinuity exactly at the boundary of the light cone.
This is different from the traveling front solutions computed in previous works in similar settings, which often exhibited sharp but continuous fronts preceded by exponential tails, therefore surreptitiously breaking causality.
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 and does not explicitly display second-order spatial derivatives, unlike the FKPP equations. 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, number, 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. 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 could smoothen the front at the light-cone boundary.

\section*{Acknowledgements}
The authors are grateful to Adam Nahum and Tianci Zhou for insightful discussions. AM thanks the hospitality of \'Ecole Normale Sup\'erieure (Paris) where this work was initiated.



\paragraph{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. 


\begin{appendix}

\section{Inter-world collision integral}
\label{app:HT_proj}
We recall the collision integral, onshell:
\begin{align}
     \tilde I_{\alpha \beta}({\boldsymbol{k}}) 
      &= 2 \,  {\rm Im} \tilde \Sigma^{R}({\boldsymbol{k}}) \, \tilde F_{\alpha \beta}({\boldsymbol{k}}) + {\rm i} \tilde \Sigma^K_{\alpha\beta}({\boldsymbol{k}}),  \label{eq:coll_inter_onshell}
  \end{align}
  where
  \begin{align}
     2 \,  {\rm Im} \tilde \Sigma^{R}({\boldsymbol{k}}) 
     & =  2
  \sum_{{\boldsymbol{k}}'{\boldsymbol{k}}''} \int \frac{{\rm d} \omega'}{2\pi} \frac{{\rm d} \omega''}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2 
  \nonumber \\
  & \qquad \times   
  \,\mathrm{Im}\, G^R(\omega'',{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}'-{\boldsymbol{k}}) \nonumber \\
  & \qquad \times  \left[ \tanh\left(\frac{\omega''}{2T}\right) +\tanh\left(\frac{\omega'-\omega''}{2T}\right) \right] \left[ \coth\left(\frac{\omega'}{2T}\right) +\tanh\left(\frac{\epsilon_{\boldsymbol{k}}-\omega'}{2T}\right) \right], 
  \end{align}
  and
  \begin{align}
     {\rm i} \tilde \Sigma^K_{\alpha\beta}({\boldsymbol{k}})  & =  2
  \sum_{{\boldsymbol{k}}'{\boldsymbol{k}}''} \int \frac{{\rm d} \omega'}{2\pi} \frac{{\rm d} \omega''}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2 \,\mathrm{Im}\, G^R(\omega'',{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'')
  \nonumber \\
  & \qquad \times   
   \,\mathrm{Im}\, G^R(\omega'-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}'-{\boldsymbol{k}}) 
  \tilde F_{\alpha\beta}({\boldsymbol{k}}'')
  \tilde F_{\alpha\beta}({\boldsymbol{k}}'-{\boldsymbol{k}}'')
  \tilde F_{\beta\alpha}({\boldsymbol{k}}'-{\boldsymbol{k}}).  \label{eq:Sigma_K_inter_onsell}
  \end{align}
Now let us inject the following ansatz in Eq.~(\ref{eq:coll_inter_onshell}),
\begin{align}
    \tilde F^{\rm ansatz}_{\alpha\beta}({\boldsymbol{k}}) &= \left[ \phi  + {\boldsymbol{u}}_{\boldsymbol{k}} \cdot \boldsymbol{\phi}_1\right] \tilde F_{\alpha\beta}^{\rm corr}({\boldsymbol{k}}).
\end{align}
The term in ${\rm Im} \tilde \Sigma^{R}({\boldsymbol{k}})$ is straight-forward:
\begin{align}
     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}}).
  \end{align}
Let us treat the term $ {\rm i} \tilde \Sigma^K_{\alpha\beta}({\boldsymbol{k}}) $ in Eq.~(\ref{eq:Sigma_K_inter_onsell}).
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
\begin{align}
 & 2 \phi^3 
  \sum_{{\boldsymbol{k}}'{\boldsymbol{k}}''} \int \frac{{\rm d} \omega'}{2\pi} \frac{{\rm d} \omega''}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2   
  \,\mathrm{Im}\, G^R(\omega'',{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'')  \nonumber\\
& \qquad\qquad \times 
\,\mathrm{Im}\, G^R(\omega'-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}'-{\boldsymbol{k}})  \tilde F^{\rm corr}_{\alpha\beta}({\boldsymbol{k}}'')
  \tilde F^{\rm corr}_{\alpha\beta}({\boldsymbol{k}}'-{\boldsymbol{k}}'')
  \tilde F^{\rm corr}_{\beta\alpha}({\boldsymbol{k}}'-{\boldsymbol{k}})  \nonumber\\
  & = 2\phi^3\, {\rm Im} \tilde \Sigma^{R}({\boldsymbol{k}}),
  \end{align}
  where we used  the trigonometric relation
\begin{align}\label{corr-id}
& F^{\rm corr}_{\alpha\beta}(\omega'')
   F^{\rm corr}_{\alpha\beta}(\omega'-\omega'')
   F^{\rm corr}_{\beta\alpha}(\omega'-\omega) \nonumber\\
  &+ 
  \left[ \tanh\left(\frac{\omega''}{2T}\right) +\tanh\left(\frac{\omega'-\omega''}{2T}\right) \right] \left[ \coth\left(\frac{\omega'}{2T}\right) +\tanh\left(\frac{\omega-\omega'}{2T}\right) \right]   F^{\rm corr}_{\alpha\beta}(\omega)=0.
\end{align}
Let us now evaluate the term of order $\phi^2 \phi_1$. We have
\begin{align}
&{\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}})\nonumber\\
     &+ 2 \phi^2
  \sum_{{\boldsymbol{k}}'{\boldsymbol{k}}''} \int \frac{{\rm d} \omega'}{2\pi} \frac{{\rm d} \omega''}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2   
  \,\mathrm{Im}\, G^R(\omega'',{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'')\nonumber\\&
  \times  \mathrm{Im}\, G^R(\omega'-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}'-{\boldsymbol{k}}) 
 \biggl\{ \boldsymbol{\phi}_1 \cdot ({\boldsymbol{u}}_{{\boldsymbol{k}}''} +{\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}''} +{\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}}) \biggr\}  F^{\rm corr}_{\alpha\beta}(\omega'')
   F^{\rm corr}_{\alpha\beta}(\omega'-\omega'')
   F^{\rm corr}_{\beta\alpha}(\omega'-\omega). 
\end{align}
The expression can be simplified by use of Eq.~(\ref{corr-id}), yielding
\begin{align}
&{\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}})\nonumber\\
     &-2 \phi^2  \tilde F^{\rm corr}_{\alpha\beta}({\boldsymbol{k}})
  \sum_{{\boldsymbol{k}}'{\boldsymbol{k}}''} \int \frac{{\rm d} \omega'}{2\pi} \frac{{\rm d} \omega''}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2   
  \,\mathrm{Im}\, G^R(\omega'',{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'')  \nonumber\\
&\times \mathrm{Im}\, G^R(\omega'-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}'-{\boldsymbol{k}}) \left[ \tanh\left(\frac{\omega''}{2T}\right) +\tanh\left(\frac{\omega'-\omega''}{2T}\right) \right] \left[ \coth\left(\frac{\omega'}{2T}\right) +\tanh\left(\frac{\epsilon_{\boldsymbol{k}}-\omega'}{2T}\right) \right] \nonumber\\
& \times  
\biggl\{ \boldsymbol{\phi}_1 \cdot ({\boldsymbol{u}}_{{\boldsymbol{k}}''} +{\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}''} +{\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}}) \biggr\}.  \label{eq:Sigma_K_away_from_crit0}
\end{align}

\subsection{Critical case}
Close to criticality, the Cooperon propagator reads
\begin{align}
    D^R(\omega,{\boldsymbol{k}}) \approx \frac{-1/\rho_{\rm F}}{r - {\rm i} a {\omega}/{T}+ \xi^2 k^2 + \ldots},
\end{align}
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\to0$, 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
\newline ${\boldsymbol{u}}_{{\boldsymbol{k}}''} +{\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}''} +{\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}}\approx -  {\boldsymbol{u}}_{\boldsymbol{k}}$. Thus we have
\begin{align}
&{\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}}) \\
      &      +  2 \phi^2 (\boldsymbol{\phi}_1 \cdot {\boldsymbol{u}}_{\boldsymbol{k}}) \tilde F^{\rm corr}_{\alpha\beta}({\boldsymbol{k}}) 
  \sum_{{\boldsymbol{k}}'{\boldsymbol{k}}''} \int \frac{{\rm d} \omega'}{2\pi} \frac{{\rm d} \omega''}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2   
  \,\mathrm{Im}\, G^R(\omega'',{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'') \nonumber\\
  &\times \mathrm{Im}\, G^R(\omega'-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}'-{\boldsymbol{k}}) \left[ \tanh\left(\frac{\omega''}{2T}\right) +\tanh\left(\frac{\omega'-\omega''}{2T}\right) \right]\left[ \coth\left(\frac{\omega'}{2T}\right)+\tanh\left(\frac{\epsilon_{\boldsymbol{k}}-\omega'}{2T}\right) \right].   \nonumber
\end{align}
Performing the sums on $\omega''$ and ${\boldsymbol{k}}''$, we get
\begin{align}
&{\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}})\nonumber\\
     & +  4 \phi^2 (\boldsymbol{\phi}_1 \cdot {\boldsymbol{u}}_{\boldsymbol{k}})
  \sum_{{\boldsymbol{k}}'} \int \frac{{\rm d} \omega'}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2   {\rm Im}\Pi^R(\omega',{\boldsymbol{k}}') \left[ \coth\left(\frac{\omega'}{2T}\right) +\tanh\left(\frac{\epsilon_{\boldsymbol{k}}-\omega'}{2T}\right)  \right]
   \nonumber\\
&\times \mathrm{Im}\, G^R(\omega'-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}'-{\boldsymbol{k}}) \tilde F^{\rm corr}_{\alpha\beta}({\boldsymbol{k}}) .
\end{align}
The above can now be written as
\begin{align}
&{\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}}). 
      \label{eq:Sigma_K_crit}
\end{align}
Altogether,
\begin{align}
    \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}}).
\end{align}

\subsection{Away from criticality}
\label{app:away_from_crit}
In the previous Subsection, we have treated the critical case which yields the following coupled PDEs 
\begin{align} \label{eq:coupled_PDEs_app}
\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},  \\[\medskipamount]
          \partial_t \boldsymbol{\phi}_1 +  v_{\rm F} \boldsymbol{\nabla}_x \phi
          \!\!\!\!& =\, \boldsymbol{\phi}_1   (\phi^2-1) / \tau_{\rm F},
          \end{array}
\right.
\end{align}
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, we separate the expression of ${\rm i} \tilde \Sigma^K_{\alpha\beta}({\boldsymbol{k}})$ in Eq.~(\ref{eq:Sigma_K_away_from_crit0}) into the critical expression computed in Eq.~(\ref{eq:Sigma_K_crit}) and the rest, in particular, the term
\begin{align}
    {\boldsymbol{u}}_{{\boldsymbol{k}}''} +{\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}''} + {\boldsymbol{u}}_{\boldsymbol{k}} + {\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}}.
\end{align}
Explicitly, we have
\begin{align} 
&{\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}})  \nonumber\\
     &-2 \phi^2  \tilde F^{\rm corr}_{\alpha\beta}({\boldsymbol{k}})
  \sum_{{\boldsymbol{k}}'{\boldsymbol{k}}''} \int \frac{{\rm d} \omega'}{2\pi} \frac{{\rm d} \omega''}{2\pi} 
  |D^R(\omega',{\boldsymbol{k}}')|^2   
  \,\mathrm{Im}\, G^R(\omega'',{\boldsymbol{k}}'') \,\mathrm{Im}\, G^R(\omega'-\omega'',{\boldsymbol{k}}'-{\boldsymbol{k}}'')  \nonumber\\
&\times \mathrm{Im}\, G^R(\omega'-\epsilon_{\boldsymbol{k}},{\boldsymbol{k}}'-{\boldsymbol{k}}) \left[ \tanh\left(\frac{\omega''}{2T}\right) +\tanh\left(\frac{\omega'-\omega''}{2T}\right) \right] \left[ \coth\left(\frac{\omega'}{2T}\right) +\tanh\left(\frac{\epsilon_{\boldsymbol{k}}-\omega'}{2T}\right) \right] \nonumber\\
& \times  
\biggl\{ \boldsymbol{\phi}_1 \cdot \left(  {\boldsymbol{u}}_{{\boldsymbol{k}}''} +{\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}''} + {\boldsymbol{u}}_{\boldsymbol{k}} + {\boldsymbol{u}}_{{\boldsymbol{k}}'-{\boldsymbol{k}}}\right) \biggr\} \label{eq:Sigma_K_away_from_crit},
\end{align}
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
\begin{align}
    \int {\rm d} \Omega_{\boldsymbol{k}} \mbox{ last term of Eq.~(\ref{eq:Sigma_K_away_from_crit})}  = 0.
\end{align}
This guarantees that the {\sc rhs} of the first equation in the coupled PDE in ~(\ref{eq:coupled_PDEs_app}) 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, \textit{i.e.}
\begin{align}
    \int {\rm d} \Omega_{\boldsymbol{k}} {\boldsymbol{u}}_{\boldsymbol{k}} \mbox{ last term of Eq.~(\ref{eq:Sigma_K_away_from_crit})}  \neq  0,
\end{align}
and therefore to give an extra contribution to the term in $ \phi^2 \boldsymbol{\phi}_1$ in the {\sc rhs} of the second equation in~(\ref{eq:coupled_PDEs_app}). 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$:
\begin{align}
 \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}}.
\end{align}
This justifies the {\sc rhs} in the second line of the coupled PDEs in (\ref{eq:coupled_PDEs}). $\gamma = 1$ corresponds to the critical case and we expect the near-critical regime to be described by $\gamma \lesssim 1$. 

\section{Discontinuous front for generic $\gamma\neq 1$} \label{app:alpha_neq_1}
In this Section, we compute the late-time solution of the following coupled PDEs
\begin{align}
\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.
\end{align}
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
\begin{align}
\left\{
\begin{array}{rl}
    \phi_+(\tau,X) & \!\!\!\! := \phi(\tau,X+\tau) \stackrel{\tau\to\infty}{\rightarrow} f_+(X), \\
    \phi_{1+}(\tau,X) & \!\!\!\! := \phi_1(\tau,X+\tau) \stackrel{\tau\to\infty}{\rightarrow} f_{1+}(X). 
    \end{array}
    \right.
\end{align}
They obey the equations
\begin{align} \label{2eq}
\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.
\end{align}
As $\tau\to\infty$, this leads to
\begin{align}
 \partial_X (f_{1+}-f_+) = -f_+(1-f_+^2)= f_{1+}(1-\gamma f_+^2).
\end{align}
Hence
\begin{align}
f_{1+}=-f_+\frac{1-f_+^2}{1-\gamma f_+^2},
\label{g1g}
\end{align}
and
\begin{align}
 \partial_X \left(f_+\frac{1-f_+^2}{1-\gamma f_+^2}+f_+\right) = f_+(1-f_+^2).
\end{align}
This is a first-order ODE that can be integrated by means of decomposition into
simple fractions. Implicitly:
\begin{align}
-\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. 
\end{align}
The constant $C$ can be determined by solving $f_+(\lambda^-)$ (\textit{i.e.} the
discontinuity).
We may extract the missing information from ``the other side'', \textit{i.e.} on the left-moving front.
Introduce 
\begin{align}
\phi_{-}  (\tau,X) = \phi(\tau,X-\tau),\qquad
 \phi_{1-}(\tau,X) = \phi_1(\tau,X-\tau).
\end{align}
Using $\partial_{\tau}\phi_- = \partial_{\tau}\phi -\partial_X \phi$, we have
\begin{align}\label{2eqbis}
\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.
\end{align}
Then, integrating the first equation:
\begin{align}
\phi_-(\tau,X)
&= {\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]
\\
&= {\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)
\\
&= {\rm e}^{-\int_0^\tau {\rm d} s\, (1-\phi_-(\tau-s,X)^2)}\phi_0(X) \nonumber \\
& \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).
\end{align}
Add $2\tau$ to $X$, and use $\phi_-(\tau,X+2\tau)=\phi_+(\tau,X)$:
\begin{align}
\phi_+(\tau,X)
&= {\rm e}^{-\int_0^\tau{\rm d} s\, (1-\phi_+(\tau-s,X+2s)^2)}\phi_0(X+2\tau) \nonumber \\
& \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).
\end{align}
When $X > \lambda$, recalling that we postulate a solution such that $\phi_+(\tau,X > \lambda)=1$, the above equation trivially simplifies to $1 = 1 - 0$.
When $X<\lambda$, we take $\tau>\frac{\lambda-X}2$. Then, since $\phi_+(\tau,X > \lambda)=1$, we have:
\begin{align}
\phi_+(\tau,X)
&= {\rm e}^{-\int_0^{\frac{\lambda-X}2}{\rm d} s\, (1-\phi_+(\tau-s,X+2s)^2)} \nonumber \\ 
& \qquad -
\int_0^{\frac{\lambda-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).
\label{12}
\end{align}
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:
\begin{align}
(\partial_X\phi_{1+}+\partial_X\phi_+)(\tau-u,X+2u)
&=\frac12\,\frac{{\rm d}}{{\rm d} u}\Big[(\phi_{1+}+\phi_+)(\tau-u,X+2u)\Big] \nonumber \\
& \qquad 
+\frac12(\partial_\tau\phi_{1+}+\partial_\tau\phi_+)(\tau-u,X+2u) \nonumber
\\&=\frac12\,\frac{{\rm d}}{{\rm d} u}\Big[(\phi_{1+}+\phi_+)(\tau-u,X+2u)\Big] \nonumber \\
& \qquad 
-\frac12\Big[\phi_+(1-\phi_+^2) +\phi_{1+}(1-\gamma\phi_+^2)\Big]
(\tau-u,X+2u), 
\label{13}
\end{align}
using Eq.~(\ref{2eq}). Inserting Eq.~(\ref{13}) into Eq.~(\ref{12}) and integrate by parts:
\begin{align}
\phi_+(\tau,X)
&= {\rm e}^{-\int_0^{\frac{\lambda-X}2}{\rm d} s\, (\cdots)} -
\frac12(\phi_{1+}+\phi_+)\left(\tau-\frac{\lambda-X}2,\lambda\right){\rm e}^{-\int_0^{\frac{\lambda-X}2}{\rm d} s (\cdots)} \nonumber \\
& \qquad\qquad +
\frac12(\phi_{1+}+\phi_+)(\tau,X) 
+\int_0^{\frac{\lambda-X}2}{\rm d} u \, (\cdots).
\end{align}
Crucially, the terms collected in $(\cdots)$ are bounded. Moreover, $\phi_{1+}$ and $\phi_+$ are continuous and $(\phi_{1+}+\phi_+)(\tau,\lambda)=1$ for any finite time $\tau$. Therefore, we can replace $(\phi_{1+}+\phi_+)(\tau-\frac{\lambda-X}2,\lambda)$ by 1 all times. Later sending $\tau\to\infty$, we obtain
\begin{align}
f_+(X)
&= {\rm e}^{-\int_0^{\frac{\lambda-X}2}{\rm d} s\, (\cdots)} -
\frac12{\rm e}^{-\int_0^{\frac{\lambda-X}2}{\rm d} s (\cdots)} +
\frac12(f_{1+}+f_+)(X) 
+\int_0^{\frac{\lambda-X}2}{\rm d} u \, (\cdots),
\end{align}
and send $X \to \lambda^-$
\begin{align}
f_+(\lambda^-) &= \frac12 + \frac12(f_{1+}+f_+)(\lambda^-)
\qquad\mbox{\textit{i.e.}}\qquad f_{1+}(\lambda^-)=f_+(\lambda^-)-1.
\end{align}
Let us call $L:=f_+(\lambda^-)$. Using Eq.~(\ref{g1g}), we have
\begin{align}
-L\frac{1-L^2}{1-\gamma L^2}=L-1.
\end{align}
Simplifying by $L-1$
\begin{align}
(1+\gamma)L^2+L-1=0.
\end{align}
The positive solution is
\begin{align}
L(\gamma)=\frac{\sqrt{5+4\gamma}-1}{2(1+\gamma)}.
\end{align}
This quantity monotonically interpolates from $L=1$ for $\gamma=-1$ to
$L=\frac12$ for $\gamma=1$, going through $L=1/(\text{golden
ratio})$ for $\gamma=0$.


\end{appendix}


