Traveling discontinuity at the quantum butterfly front

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


Introduction
Quantum information scrambling is the mechanism by which localized information in an extended closed quantum many-body system with local interactions flows to non-local degrees of freedom, becoming practically irretrievable. In practice, these information dynamics can be conveniently probed by means of out-of-time-ordered correlators (OTOCs) such as squared commutators of operators inserted at different space-time points, e.g. C(t, x ) := 〈[O(t, x ), O(0, 0)] 2 〉. The effective loss of information is characterized by (i) a ballistic spread of information, often dubbed as the "quantum butterfly effect", (ii) a growth regime, reminiscent of the exponential separation between nearby trajectories in classical chaotic systems, (iii) a purely quantum saturation regime beyond a scrambling time t * .
The ballistic spread of quantum information has been firmly established on the basis of the Lieb-Robinson bound [1]. The causal light-cone structure, with a wavefront traveling at a model-dependent butterfly velocity v B , was confirmed in a wide variety of models, from noninteracting 1d 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: where λ L is the Lyapunov exponent. For truly quantum systems, the rapid growth concentrated near the light cone boundary is understood to be set by model-dependent microscopic scales, and significant efforts were made to compute the particular shape of the butterfly front in a variety of models.
Wavefronts described by power laws, sometimes oscillatory, were found in free and integrable models [2,3]. Sharp wavefronts were found in interacting spin chains [4], nonintegrable systems with diffusive transport [5,6], as well as large-N or semiclassical models [7,8] and holographic models [9,10]. Interestingly, random unitary circuits without conserved quantities, i.e. in the absence of diffusive transport, were found to develop broad fronts controlled by a diffusively growing length scale [11][12][13][14]. Notably, several works [15][16][17][18] have Figure 1: Butterfly effect in the clean interacting metal defined in Eqs. (5) and (22): after a small and localized perturbation at time τ = 0, the space-time dynamics of the inter-world distribution function F du defined in Eq. (15) follows a causal light-cone structure with a front traveling at the maximal butterfly velocity v B = v F / d where v F is the Fermi velocity and d the spatial dimensions. Here, we sketch φ(τ, X ), the component of F du averaged over the Fermi surface, to be introduced in Eq. (31). At early times, the exponential growth regime is governed by the inelastic scattering rate. In the late-time limit, scrambling is governed by shock-wave dynamics with the traveling front developing a distinct discontinuity followed by an exponential decay on the scale of the mean free path ℓ set by the electronic interaction.
pointed towards an effective description of the front in terms of partial differential equations belonging to the class of the Fisher Kolmogorov-Petrovsky-Piskunov (FKPP) equation [19,20] which is known to exhibit traveling wave solutions (see Ref. [21] and references therein).
In this manuscript, we address the dynamics and the geometry of the wavefront for a locally interacting system in the vicinity of a continuous classical phase transition corresponding to the spontaneous breaking of the symmetry associated with a conserved quantity. Practically, we consider a paradigmatic model of interacting electrons where the interactions are due to strong superconducting fluctuations. The long wavelength fluctuations close to criticality produce a separation of scales which we leverage to derive analytic results. How the results get modified on moving away from criticality is transparent in our derivation and our approach can be systematized to address other near-critical quantum many-body systems.
We compute the OTOCs by means of an augmented Keldysh formalism, the so-called manyworld formalism, which was originally proposed in Ref. [22] and recently used to derive kinetic equations for the spreading of quantum information in fermionic interacting systems including electron-phonon scattering, electrons-impurity scattering, as well as electron-electron scattering [15]. The augmented Keldysh formalism has been used in other recent works as well [23][24][25][26]. This formalism can be conveniently harnessed to the standard field-theoretic concepts, tools, and approximation schemes that have been developed over the years in condensed matter theory. Here, we treat the interaction between the electrons and the superconducting fluctuations by means of the random-phase approximation (RPA) in the particleparticle 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.

Summary and main results
The paper is organized as follows. In Sect. 2, we quickly review the many-world formalism which is used to compute OTOCs and access quantum chaotic features of many-body systems. It generalizes the standard Keldysh formalism by studying two replicas of the theory, the socalled worlds, which are only correlated through their initial conditions. In particular, we introduce the inter-world distribution function F αβ (ω, k; t, x ), where α ̸ = β are the world indices, which quantifies the amount of correlation between the two worlds.
In Sect. 3, we introduce a paradigmatic model describing interacting electrons close to a superconducting transition. We avoid the technical challenges of approaching the critical point from within the superconducting phase and work in the normal phase where no long-range order develops. We derive the corresponding kinetic equation for the inter-world distribution which is shown to be of the form with the non-linear collision integral I αβ given in Eq. (30).
In Sect. 4, we propose and numerically validate an ansatz for F αβ leading to a simplified set of non-linear partial differential equations (PDEs) involving two fields φ(t, x ) and φ 1 (t, x ). These are the first terms of the partial-wave expansion in momentum space of F αβ (ω, k; t, x ) evaluated on-shell and at the Fermi surface, i.e. at ω = ε k and k → k F . In terms of dimensionless quantities X for space and τ for time, the PDEs read where the parameter γ effectively encodes the distance from the superconducting phase transition: γ = 1 corresponds to criticality and 0 < γ < 1 corresponds to the off-critical regime in the normal phase.
In Sect. 5, starting from a generic localized initial perturbation, we analytically solve for the dynamics of φ and φ 1 , in any dimension d, at criticality as well as away from criticality. The results are sketched in Fig. 1. The relaxation of the inter-world distribution is found to strictly occur within a light cone growing from the initial perturbation at a constant butterfly velocity v B = v F / d where v F is the Fermi velocity. We work out the early-time dynamics with an exponential growth of scrambling which is controlled by the inelastic scattering rate. Importantly, in the late-time regime, we find that scrambling is governed by shock-wave dynamics, with a traveling wave that develops a discontinuous radial front of the form which extends over a length scale set by the mean free path ℓ related to the scattering of the electrons by superconducting fluctuations. Inside the light cone, f + dies off exponentially away from its boundary (|x | − v B t < 0), and f + = 1 outside the light cone (|x | − v B t > 0). Notably, we find that f + is discontinuous precisely at the boundary (|x | − v B t = 0). We also work out explicitly the exponential falloff governing the approach to the saturation regime within the bulk of the light cone. We conclude in Sect. 6 by discussing the relations of our results to previous works and by giving future directions.

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 Figure 2: Two-world Keldysh contour 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 ψ and ψ † correspond to the OTOC in Eq. (1).
where ρ 0 is the initial density matrix at time t = 0 which is normalized as Tr ρ 0 = 1, H is the Hamiltonian generating the dynamics, and ψ(x ) is a local operator at position x . We have in mind fermionic annihilation/creation operators, but the discussion can be easily adapted to the bosonic case. Here, the four operators ψ(0), ψ(x ), ψ † (0), and ψ † (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 Schwinger-Keldysh formalism defined on the two-fold Baym-Kadanoff contour to a formalism on a four-fold contour, suitable to compute four-point OTOCs. This was originally introduced in Ref. [22] and we refer the reader to Ref. [15] for a detailed presentation which we simply follow. The OTOC in Eq. (1) involves two forward and two backward time-evolution operators. Therefore, following the standard Schwinger-Keldysh construction [27], this yields the four-fold time-ordered contour C depicted in Fig. 2. The forward (backward) branches are labeled by an index a = + (a = −). The first two branches are said to belong to the "up world" and are labeled by the index α = u. The two other branches, posterior on the contour, correspond to the so-called "down world", and are labeled by α = 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. (1) can be rewritten as where T C is the time-ordering operator on the contour C, the operators ψ and ψ † are now written in the Heisenberg picture, and 〈. . .〉 := Tr [. . . ρ 0 ]. In our subsequent study of the quantum butterfly effect, we shall assume equilibrium conditions: the system is initially prepared in thermal equilibrium at temperature T , i.e. with the Gibbs state ρ 0 ∼ e −H/T where we set k 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: with the Keldysh indices a, b = +, −, the world indices α, β = u, d, and where the mixed operatorŜ 0 := ψ − † u (0, 0)ψ − d (0, 0) results from clubbing of the two operators in Eq. (2) that are both evaluated at time t = 0 but at disconnected locations on the time contour C. The OTOC in Eq. (1) is recovered by taking α = d and β = u, while a = ± and b = ± can be chosen arbitrarily. If one is to interpret iG ab αβ (t, x ; t ′ , x ′ ) as a two-point function rather than a four-point function,Ŝ 0 has to be understood as a local modification of the initial condition which, a priori, acts differently on each world: one particle is added to the up world at position x = 0 while one particle is removed from the down world. Other choices of operator content are possible forŜ 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: with the infinitesimal parameter |δφ 0 | ≪ 1. Importantly, an inspection of Eq. (3) forŜ 0 = 1 shows that the intra-world Green's functions (i.e. α = β) correspond to the standard (i.e. single-world) Schwinger-Keldysh two-point correlators: From now on, given that the intra-world quantities are identical in both α = u, d worlds, we simply drop the repeated world indices, e.g. G αα → G, except when this obscures the meaning. Furthermore, when using the indices αβ, we specifically mean α ̸ = β unless specified otherwise.

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 where H 0 is the non-interacting part and the interaction in H int depends on the specific problem at hand, see Eq. (22). The fermionic operator c † kσ creates a electrons with spin σ =↑ or ↓ (σ =↓, ↑) and momentum k in the Brillouin zone (BZ). ε 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 µ can be included via the substitution ε k → ε k − µ. For simplicity, we shall assume that the Fermi surface is spherical, i.e. ε k = 0 when k → k F . Whenever this does not harm the understanding, we shall simply drop the spin indices.

Green's functions in the Keldysh basis
The 16 Green's functions G a b αβ 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 C is such that the inter-world Green's functions (i.e. α ̸ = β) do not depend on the ± basis: On the other hand, the intra-world Green's functions (i.e. α = β) are related via [27] It is therefore customary to perform a rotation from the a = ± basis to the so-called Keldysh basis and work with retarded, advanced, and Keldysh Green's functions: 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 Σ R is the retarded component of the self-energy. It is due to the interaction in H int and can be computed diagrammatically within the standard Schwinger-Keldysh formalism. The last equality is the expression of the fermionic fluctuation-dissipation theorem with F (ω) := tanh (ω/2T ). Concerning inter-world Green's functions, the relation in Eq. (7) together with Eq. (9) immediately implies 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 K ud and G K du which also depend on the choice of the operatorŜ 0 . For the choiceŜ 0 = 1, the two worlds have the same thermal initial conditions, and one may check that G K αβ are space-and time-translational invariant and However, for a generic choice ofŜ 0 , G K αβ is not guaranteed to be space-and time-translational invariant and it can be determined via the Schwinger-Dyson equation reading where Σ K αβ is the Keldysh component of the inter-world self-energy which can be computed diagrammatically in the many-world formalism.

Inter-world kinetic equation
It is useful to work in the Wigner representation and parameterize G K αβ in terms of the real function F αβ : where we introduced the Moyal product ⋆ := exp i 2 ( 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 αβ is dubbed the inter-world distribution function. One has F ud (ω, k; t, x ) ∈ [−2, 0] and F du ∈ [0, 2]. Assuming that the variations of F αβ occur on scales much larger than the microscopic scales involved in G R , we may work in the so-called quasi-classical or gradient approximation which consists in truncating the derivative expansion to its leading terms. The interworld distribution F αβ gives access to the information scrambling as it is related to the original OTOC in Eq. (2). Let us briefly outline this connection. At late times, we expect the decoupling with n 0 := Tr ψ † (0)ψ(0)ρ 0 and where the interworld quantities are computed from the interworld Schwinger-Dyson equations in the presence of a local perturbation to the correlatedworld solution at x = 0 and time t = 0. Massaging Eq. (15) by acting on both sides with the inverse of G R and G A , and using the Dyson equations (10) and (13), one derives the kinetic equation on F αβ which is analogous to the standard (intra-world) kinetic equation: where the LHS corresponds to the non-interacting physics set by H 0 , with the velocity v k :=∇ k ε k , and the right-hand side (RHS) is the so-called collision integral which stems from H int and reads Let us recall that Keldysh components and collision integrals depend on space and time, namely x and t, through the inter-world distribution functions. However, here and from now on, we simplify the notation by dropping the explicit dependence on these objects. As discussed above, the intra-world quantity Σ R cannot depend on F αβ , however Σ K αβ is expected to be a non-linear functional of F αβ . The inter-world kinetic equation in Eq. (17) is therefore a non-linear partial integrodifferential equation. It has a trivial steady-state solution 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Ŝ 0 = 1 in Eq. (12), where both worlds evolve coherently, corresponds to another steady state characterized by and which is dubbed the "correlated-world" solution. As we shall see explicitly later, the correlated-world solution atŜ 0 = 1 is expected to be unstable against small perturbations ofŜ 0 and the only stable steady-state is the uncorrelated-world solution. Figure 3: Superconducting self-energy Σ and polarization bubble Π of the model in Eq. (22) treated in the RPA scheme in the particle-particle channel. The expressions in the Keldysh formalism are given in Eqs. (25) and (24), respectively.
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 where 0 ≤ δφ 0 (x ) ≪ 1 is an initial perturbation to the correlated-world solution localized on a compact support around x = 0. For simplicity, we consider a perturbation that is the same for both F ud and F du .

Superconducting fluctuations 3.1 Model
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. (5) reads where U < 0 is an attractive interaction facilitating superconductivity. In dimensions d ≥ 2, this model exhibits a finite-temperature phase transition towards a superconducting phase associated with the spontaneous breaking of the U(1) symmetry. Here, we consider the nearcritical regime, above the critical temperature where the U(1) symmetry is not broken but the superconducting fluctuations are sizable. The (intra-world) physics of this model is well understood, and we rely on standard and well-tested methods which we extend to the many-world formalism. In practice, we decouple the Hubbard interaction in the Cooper channel and obtain a theory of fermions coupled to bosonic fluctuations. The Cooperon Green's functions within RPA in the particle-particle channel, are given by [28][29][30][31] where Π R is the retarded Cooper bubble and the last equality is the bosonic fluctuationdissipation theorem with P(ω) := coth (ω/2T ). The retarded Cooper bubble and the electronic self-energies within RPA, i.e. limiting ourselves to the one-loop diagrams depicted in and The real parts of the retarded components can be recovered via the Kramers-Kronig relation.
A self-consistent treatment of Green's functions, self-energies, and bubbles ensures that the RPA scheme is a conserving approximation. It is known to be exact in the large N -limit, where N is the number of electronic orbitals [28][29][30][31].
Close to criticality, the Cooperon propagator reads [28][29][30][31][32] where ρ F is the density of states at the Fermi energy, and the parameter r ∝ (T − T c )/T c is the detuning from the critical point. In addition, the remaining parameters are all positive where v F is the Fermi velocity. At criticality r → 0, the Cooperon becomes soft with diverging length scale l ∼ 1/r ν and timescale ∼ l z (here ν = 1/2, z = 2), and the propagator is singular at ω = k = 0.

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. (18). 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: Within the RPA scheme, the inter-world Keldysh components read Altogether, this yields the inter-world kinetic equation (17) with the collision integral We recall that we omitted the local space and time dependence of I αβ and F αβ 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. (19) and the correlated-world solution in Eq. (20). 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 αβ (such as number, energy, momentum conservation) that guarantee sum rules such as dω k I(ω, 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.

Inter-world kinetics in the near-critical regime
In this Section, we propose and numerically validate an ansatz to the inter-world distribution F αβ which allows us to derive a much simpler version of the kinetic equation (17) with its collision integral in Eq. (30). This is done in the vicinity of the superconducting transition where a clear separation of energy scales can be made. The resulting effective description of the information scrambling dynamics consists of a set of coupled PDEs that we solve analytically in Sect. 5.

Partial-wave ansatz
We start by making the quasi-particle approximation: In all practical instances, F αβ (ω, k) appears multiplied by the density of states Im G R (ω, k), see e.g. Eq. (30). When quasi-particles are well defined, with a dispersion relation ε k , the density of states is sharply peaked around ω = ε k and one may seamlessly exchange F αβ (ω, k) with the on-shell quasi-particle distribution functionF αβ (k) := F αβ (ω = ε k , k). From now on, we use the tilde notation to denote the on-shell prescription ω = ε k . Furthermore, given that relaxation is dominated by the electronic states around the Fermi level, at energy scales (e.g. temperature) that are much smaller than the Fermi energy, one may focus on the distribution function close to the Fermi surface by subsequently setting k → k F .
We now propose to simplify drastically the partial integrodifferential kinetic equation in Eq. (17) with the following ansatz: where the unit vector u k := k/k. The two fields φ and φ 1 can be understood as the first terms of a partial-wave expansion [27] ofF αβ , accounting for its isotropic and anisotropic contributions in momentum space: where S d−1 := dΩ k is the surface area of the d − 1-sphere with unit radius and dΩ k is the elementary solid angle in the direction of k. We did not include higher-order terms in the ansatz, e.g. of the form u i k u j k φ i j 2 . Let us give the rationale behind the ansatz proposed in Eq. (31). Firstly, let us note that standard (intra-world) approaches consist in perturbing the distribution function around its equilibrium value, F = F eq + δF , and linearizing the collision integral accordingly. This approach relies on the fact that the equilibrium distribution F eq (ω, k) = tanh(ω/2T ) is a stable steady state of the (intra-world) kinetic equation, guaranteed by the H-theorem. The interworld case is much different as the initial condition set by F corr αβ is unstable and one cannot propose a perturbative ansatz. This explains whyF corr αβ appears multiplicatively in Eq. (31) and why the collision integral cannot be, a priori, linearized. In that regard, it is similar to the ansatz used in Ref. [15].
Secondly, close to the Fermi surface which is assumed to be spherical, the solutionsF corr αβ (k) andF uncorr αβ (k) do not depend on the direction of the momentum k, but only on its norm k ≈ k F . This means that we aim at describing the dynamics ofF αβ from a momentum-space isotropic and real-space homogeneous (unstable) solutioñ to another momentum-space isotropic and real-space homogeneous (stable) solutioñ 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 φ 1 which can be seen as the minimal ingredient to allow for spatial relaxation.
Thirdly, let us note thatF ansatz αβ (t, x ; k) depends on k only throughF corr αβ (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 φ and φ 1 are common to bothF ud andF du . This stems from our choice of initial perturbation in Eq. (21).

Simplified kinetic equation: coupled PDEs
We proceed by injecting the ansatz (31) in the inter-world collision integral in Eq. (30) 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. (30), this means that the term |D R (ω ′ , k ′ )| 2 diverges as ω ′ ≈ 0 and k ′ ≈ 0. The details of the computation are given in Appendix A. We obtaiñ Injecting the ansatz in the kinetic equation (17), dividing byF corr αβ (k), we get We now project on the momentum-space isotropic and first partial-wave contributions by acting with 1 where v 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 k) Note that the temperature dependence of the original problem enters through τ F and γ. γ is a dimensionless parameter that generalizes the computation performed at criticality, for which γ = 1, to near critical regimes for which γ < 1 (see Appendix A.2). After an appropriate rescaling of space and time, together with φ 1 → dφ 1 , the coupled PDEs now only involve dimensionless quantities Importantly, the generic d-dimensional case can be reduced to an effective one-dimensional case. Indeed, assuming spherically-symmetric initial conditions, we may work with the radial coordinate r: φ(τ, r) and φ 1 = φ 1 (τ, r)u r at all times. The coupled PDEs now read The set of PDEs in (37) and the following expressions in Eqs. (40), (41) are one of the main results of this manuscript. Overall, this represents a considerable simplification from the original partial integrodifferential kinetic equation (17) governing the dynamics of the inter-world distribution function F αβ (ω, k; t, x ) with the collision integral in Eq. (30).
To provide a first intuitive understanding of the previous set of PDEs, let us briefly neglect spatial inhomogeneities of φ and φ 1 , and work in d = 1. The equation on φ(t) becomes an autonomous first-order ODE, reading This is a gradient descent in the potential V (φ) = 1 2 φ 2 − 1 4 φ 4 . Reinstating the original units, the rate of escape from the correlated-world solution at the unstable extremum φ = 1, to the uncorrelated-world solution at the global minimum φ = 0 is t * := − 1 2 τ F log δφ 0 where δφ 0 := 1 − φ(0) ≪ 1. At early times, the growth of the perturbation is exponential, while it saturates at late times, Again, the decoupling of φ 1 in such a spatially homogeneous setting is evidence that the momentum anisotropy captured by φ 1 is a minimal ingredient necessary to allow for spatial propagation of the relaxation from the correlated-world solution to the uncorrelated-world solution. This will be the topic of Sect. 5. In the general case, i.e. in the presence of spatial inhomogeneities, it is instructive to compare the inter-world situation to the (standard) intra-world kinetic equations. When an intra-world distribution function F is associated with a conserved quantity (e.g. number of particles or energy), the corresponding hydrodynamic equation is typically expected to display diffusive behavior. Indeed, the timescale associated with the conserved quantity is much slower than the other modes: those can be effectively replaced by their local-equilibrium value in a fixed background of F , typically resulting in a diffusive term of the type ∇ 2 X F . Here, in the inter-world case, the distribution function F αβ is not associated with a conserved quantity (until proven otherwise) and there is no clear separation of timescales in the PDEs (37) governing the dynamics of φ and φ 1 . Consequently, one cannot a priori apply the standard hydrodynamic approach, and one must solve for the dynamics φ and φ 1 on an equal footing.

Validating the ansatz numerically
We perform two independent checks of the ansatz proposed in Eq. (31) by comparing, on the one hand, the solutions of the inter-world kinetic equation in (17) computed with the fullfledged collision integral in Eq. (30) with, on the other hand, the solutions of the coupled PDEs in (37). The first check is performed at early times in a near-critical regime while the second check is performed at later times and at a finite distance from criticality.
Numerically solving the kinetic equation is a formidable task which we simplify as much as possible by working in one dimension, d = 1, with a regular lattice dispersion ε k = − cos(k) for k ∈ [−π, π), and a point-like Fermi surface located at the wave-vector k F = π/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 ≥ 2, with the RPA treatment that we have set up in Section 3 not suited for d = 1. However, the objective here is to put to test the ansatz in conditions that are qualitatively similar to d ≥ 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 and where the Cooperon Green's function D R (ω, k) is computed following Eq. (23). Γ > 0 sets a bare fermionic inverse lifetime. In practice, Γ 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:F αβ (k) := F αβ (ω = ε k , k).

Early times
The first test consists in numerically solving the kinetic equation (17) with the full-fledged collision integral in Eq. (30) for as long as we can ensure numerical convergence of the solutions. In practice, this is challenging and we can only access early times, i.e. on the order of fractions of τ F . Therefore, in order to benchmark the ansatz in both regimes φ ∼ 1 and φ ∼ 0, we work with an initial condition that simultaneously spans those two regimes. We choose an initial condition with a perturbation ofF corr αβ (k) in the shape of a Gaussian droplet of large amplitude δφ 0 ≲ 1 and width R 0 , and localized around X = 0. Explicitly, rescaling time and space according to Eq. (39), we take the following symmetric initial conditioñ and where Θ(X ) is the Heaviside step function. We found such a Gaussian-shaped droplet, defined on the support [−3R 0 , 3R 0 ], to be easier to time-evolve numerically than a semi-circular droplet with sharp edges. The physical parameters are chosen such as to be close to criticality (on the disordered side) and to obey the hierarchy |U|, ε F ≫ T ≫ Γ , where ε F is the Fermi energy.
In Fig. 4, we compare the solutions φ(τ, X ) of the corresponding coupled PDEs with the so-lutionsF du (τ, X ; k) of the full-fledged kinetic equation. The comparison is made by extracting the first partial-wave contributions according to Eq. (32) which, in 1D, simply reads up to the time τ = 0.3. The qualitative agreement is excellent. Notice the splitting of the initial central perturbation into both a left-moving and a right-moving front. We repeated this analysis in a wide range of near-critical parameters and initial conditions and consistently found excellent agreement, even at a finite distance from criticality, in the presence of fast bosonic fluctuations. This validates the ansatz in Eq. (31) at early times.

Late times and partial-wave truncation
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 conditioñ on the kinetic equation side and, correspondingly, on the side of the coupled PDEs which now read Note that these coupled ordinary-differential equations correspond to the discussion around Eq. (42). Notably, Eq. (50) predicts that the relaxation dynamics of φ (but not those of φ 1 ) are independent of the distance to criticality, parameterized by γ. Therefore, complementary to the previous benchmark in Sect. 4.3.1, we test the ansatz at a finite distance from criticality (on the disordered side) by choosing an off-critical set of physical parameters U, T , and Γ . In Fig. 5, we compare the solutions φ(τ) and φ 1 (τ) to the corresponding coupled PDEs with the solutionsF du (τ, k) to the kinetic equation. The comparison is made by extracting the first partial-wave contributions according to Eq. (47). 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 φ 1 (τ) was made by manually adjusting the off-critical value for γ given in the caption. This validates the partial-wave truncation which is made in the ansatz.

Dynamics of information scrambling
In this Section, we solve the set of coupled PDEs (40) 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.

Early-time exponential growth
At early times, the solutions to the inter-world kinetic equation and to the simplified coupled PDEs are expected to be strongly dependent on the system parameters and the initial conditions.
Here, we solve the coupled PDEs in the linear regime around the correlated-world solution, φ ≈ 1 and φ 1 ≈ 0. We expect this linear regime to be all the more valid as the initial perturbation will be small, hence taking a longer time to reach the non-linear regime. To that end, we consider the droplet-shaped initial condition where the perturbation 0 < δφ 0 (X ) ≪ 1 is non-vanishing on a small compact ball of radius R 0 around X = 0. Note that this is a different regime from the numerics presented in Sect. 4.3.1 where the initial condition was probing the non-linear regime. To simplify, we consider the case d = 1. The linearized coupled PDEs on δφ := 1 − φ and φ 1 read This yields the following linearized PDE on δφ(τ, X ) with δφ(τ = 0, X ) = δφ 0 (X ) and ∂ τ δφ(τ = 0, X ) = 2δφ 0 (X ). Integrating the above equation over the whole space, or equivalently the first equation in (52), and introducing the integrated perturbation δM (τ) := dX δφ(τ, X ), we find an exponential growth of the perturbation, echoing the onset of chaos: with the typical timescale to escape from the unstable solution given by τ * := − 1 2 log δM 0 where δM 0 := dX δφ 0 (X ). The corresponding growth rate, λ L = 2/τ F in the original units, can be interpreted as a Lyapunov exponent. Note that the dependence on γ, the parameter quantifying the distance to criticality, has dropped. A more sophisticated calculation restricted to 1 ≪ τ ≪ τ * and |X | ≪ τ 3/4 yields the following solution to Eq. (53) See Appendix B for the details of the computation. The above expression involves an exponential growth and a diffusion kernel. It is similar to the solution obtained in the context of the O(N ) model in Ref. [8] (see Eq. 1.13 therein). Let us provide a quick way to justify this solution. We first absorb the exponential growth of δφ that was found in the solution (54) by working in terms of g(τ, X ) := e −2τ δφ(τ, X ). Using Eq. (53), it obeys with g(τ = 0, X ) = δφ 0 (X ) and ∂ τ g(τ = 0, X ) = 0. Next, given that the solution of the above PDE is expected to vary slowly with time, especially at large times, we may neglect the ∂ 2 τ g term. This yields a simple diffusion equation which, once we switch back to working with δφ, is solved by the expression in (55). Note that the solution in (55) can also be seen as the solution to the diffusion+growth equation: that appears, notably, in the context of branching Brownian motion where it describes the evolution of the expected density of particles [33].

Late-time solutions and discontinuous front
Here, we access the late-time inter-world distribution by analytically solving the coupled PDEs (40) 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.

Traveling front
The solution develops a traveling front located on a sphere of increasing radius. Our goal is to compute its velocity and shape. The first step is to notice that the late-time position of the front is, by definition, far from the origin and we may neglect the 1/r term in the RHS of Eq. (41). 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 γ = 1. However, the generic case for γ ̸ = 1 can also be solved by similar techniques and we refer the reader to Appendix C for the corresponding detailed computation. Introducing the fields the coupled PDEs can be cast as with the initial conditions, The method of characteristics gives the implicit solutions which yield Let us now assume that a right-moving front, traveling at velocity 1 (in units of v F / d), develops at late times, i.e.
pointwise. Naturally, there is also a symmetrical left-moving front. We start from where in the second line we performed a change of dummy variable s → τ−s. In the long-time limit τ → ∞, we get where we used the property φ 0 (X + 2τ) → 1. Consistently with the Eq. (62), we postulate 1 that the front is such that Let us now work with X < R 0 . Hence, lim which implies that the first term in Eq. (66) is zero. After the change of variables X + 2s → u, we are left with At X X <R 0 −→ R 0 , this yields f + (R 0 ) = 1/2. Given that we assume f + (X > R 0 ) = 1, this signals a discontinuity in f + (X ) at X = R 0 . Moreover, Eq. (68) implies that f + (X ) obeys which can be solved by separation of variables, yielding the discontinuous right-moving front with the shape Note that, here, the information on the precise shape of the initial perturbation is lost. Except for a trivial spatial offset in the shape of the front, R 0 drops out of the problem. This is expected for generic initial conditions defined on a compact support |X | < R 0 . Indeed, one can show that the first non-vanishing derivative ∂ (n) X φ 0 (X = R − 0 ) > 0 is responsible for the generation of a first-order derivative f + (X = R − 0 ) > 0 which grows exponentially with time, therefore creating a discontinuity at large times. This independence of the steady-state solution with respect to the initial condition is the hallmark of a universal solution. Notably, as previously discussed in Sec. 5.1, R 0 and the precise shape of the initial perturbation are however controlling the timescale for the traveling front to form. Barring this point, the universal features of the steady state may be safely accessed by sending R 0 → 0 + after τ → ∞. We obtain the right-moving traveling front with The last equation above follows from φ 1 = ϕ − ψ/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 with the butterfly velocity given by v B := v F / d and the radial front shape f + given by Eq. (72) that spans over a length scale controlled by the mean free path ℓ := v F τ F . We recall that v F is the Fermi velocity and the scattering time τ F was defined in Eq. (38). This solution can also be generalized to cases away from criticality, i.e. for a generic γ ̸ = 1. The qualitative features are found to be very similar to the critical case γ = 1. In the Appendix C, we show that the discontinuity of the front is now from f + (X = 0 − ) = L(γ) to f + (X > 0) = 1 with This quantity monotonically interpolates from L = 1/2 for γ = 1 to L = 1/(golden ratio) for γ = 0. The discontinuous traveling wavefronts in Eqs. (72), (73), and their generalization to noncritical regimes in Eq. (74) 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 γ = 1 (critical regime). In Fig. 6, we display the numerical solution starting from the following droplet-shaped initial perturbation of the correlated-world solution where δφ 0 sets the amplitude of the droplet, and R 0 sets its radius. This illustrates unambiguously the spatial growth of the loss of quantum coherence as time goes on. The light cone

Saturation inside the light cone
Inside the light cone and far enough from its boundaries, we have φ ≪ 1 and we can neglect the non-linearities in the RHS of the coupled PDEs (59). Working in the d = 1 case, the linearized PDEs read Assuming symmetric initial conditions, i.e., 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 where the function f ϕ (X ) should in principle be determined by solving the early-tointermediate time problem. In practice, we can determine the asymptotic behavior of the function f ϕ (X ) at large X by requiring a matching to the left side of the late-time front computed previously: where f + (X ) was computed in Eq. (72). Using the late-time solution in Eq. (77), we have The first term vanishes and we are left with which is solved as This yields, at large times and far inside the boundary of the light cone, i.e. for τ − |X | ≫ 1, These solutions can also be obtained directly from injecting the expression (72) for the traveling front in Eq. (61), using the spatial symmetry lim Note that these asymptotic solutions have already lost the information about the initial condition. Moreover, they imply the relation φ 1 (τ, X ) = −2∂ X φ(τ, X ) for τ − |X | ≫ 1. Once re-injected in the linearized PDEs, this yields the following linearly-driven diffusion equation for φ(τ, X ) only where the diffusive constant in front of the ∂ 2 X φ term is D := 2v 2 F τ F once the original units of time and space are reinstated. We note that the diffusive LHS of Eq. (83) is similar to the one of the FKPP-like equation that was derived in a similar context in Ref. [15], see Eq. (76) therein.

Discussion and conclusion
Starting from the microscopic Hamiltonian of a d-dimensional quantum many-body system of interacting electrons close to a superconducting phase transition, we carefully derived the corresponding dynamics of quantum information scrambling.
Quite expectedly, we found a ballistic spread of information governed by a non-universal butterfly velocity v B . We presented analytical solutions in the different regimes relevant to quantum chaotic dynamics: the early exponential growth, the geometry of the late-time front, and the saturation within the light cone.
Perhaps the most striking result of our work is the fact that the scrambling of quantum information at late times is governed by shock-wave dynamics. Scrambling propagates at the maximum velocity allowed by causality and develops a distinct discontinuity exactly at the boundary of the light cone. Notably, these dynamics scrupulously respect causality: scrambling does not leak outside of the light cone. This is different from the findings of previous works in similar settings, which often exhibited sharp but continuous fronts preceded by exponential tails. At a formal level, this difference arises from the fact that our effective dynamics are governed by a set of coupled PDEs that does not belong to the reaction-diffusion class.
While our full solution does share strong similarities with the diffusive FKPP solutions in the linearized regimes (at early times or deep in the bulk of the light cone), we attribute this explicit absence of a diffusive term to the fact that information dynamics is not directly associated with a conserved quantity, unlike usual transport which is associated with, say, charge, energy, or momentum conservation.
Addressing the robustness of the finite spacetime discontinuity in the traveling front is perhaps one of the most pressing questions raised by our results. Let us first note that the discontinuity has to be understood on the scale of the mean free path set by the inelastic scattering, ℓ := v F τ F , and not on the microscopic scale 1/k F where the kinetic equation breaks down. The discontinuity is unlikely an artifact of the one-loop RPA scheme (exact in the limit where the number of electronic orbitals is sent to infinity), nor of our truncated partial-wave expansion in momentum space (key to the ensuing strict causal structure, and likely to be a very good approximation in some geometries). However, higher-order terms in the Moyal product expansion that were neglected in the quasi-classical approximation, or any source of noise or disorder [34], e.g. caused by elastic scattering on random impurities, could smoothen the front at the light-cone boundary and possibly decrease the butterfly velocity. Numerical confirmation of the shock-wave dynamics starting from a microscopic Hamiltonian is expected to be difficult as the scrambling front propagates fast and the discontinuity only develops at late times. This means that one should simulate large systems of linear size L ∼ 10 ℓ until long times t max ∼ 10 τ F . Tensor network approaches in 1d may be suited to meet the challenge [4].
Now let us inject the following ansatz in Eq. (A.1), The term in ImΣ R (k) is straight-forward: Let us treat the term iΣ K αβ (k) in Eq. (A.3). It produces terms in φ 3 , φ 2 φ 1 , φφ 2 1 , and φ 3 1 . In practice, consistently with our choice of ansatz which consists of tracking only the two first multipolar contributions toF (k), we discard the terms of order φ 2 1 and φ 3 1 which yield higherorder multipolar contributions to the collision integral. The term in φ 3 reads 2φ 3 where we used the trigonometric relation Let us now evaluate the term of order φ 2 φ 1 . We have The expression can be simplified by use of Eq. (A.7), yielding

A.1 Critical case
Close to criticality, the Cooperon propagator reads At criticality r → 0, the Cooperon becomes soft with a diverging length scale l ∼ 1/r ν (here ν = 1/2), and the propagator is singular at ω = k = 0. In this case, we can approximate the term Performing the sums on ω ′′ and k ′′ , we get The above can now be written as Altogether, we obtaiñ ImΣ R (k) . (A.14)

A.2 Away from criticality
In the Subsection A.1 above, we have treated the critical case which yields the following coupled PDEs where v F is the Fermi velocity and we defined the timescale τ F as 1/τ F := −2 ImΣ R (k F ). Away from criticality, assuming 1/τ F > 0, we separate the expression of iΣ K αβ (k) in Eq. (A.9) into the critical expression computed in Eq. (A.13) and the rest by simply writing u k ′′ + u k ′ −k ′′ + u k = −u k + u k ′′ + u k ′ −k ′′ + u k + u k ′ −k .
(A. 16) Explicitly, we have iΣ K αβ (k) = −2φ 2 φ + (φ 1 · u k ) F corr αβ (k)ImΣ R (k) − 2φ 2F corr αβ (k) where the first term is the critical expression while the second term collects the rest. Assuming that ε k = ε −k , one may check that this second term is odd under k → −k and therefore cannot contribute to the projection on the momentum-space isotropic contribution and therefore to give an extra contribution to the term in φ 2 φ 1 in the RHS of the second equation in (A. 15). 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 γ > 0: This justifies the RHS in the second line of the coupled PDEs in (37). γ = 1 corresponds to the critical case and we expect the near-critical regime to be described by γ < 1.