# Universal out-of-equilibrium dynamics of 1D critical quantum systems perturbed by noise coupled to energy

Alexios Christopoulos, Pierre Le Doussal, Denis Bernard, Andrea de Luca

► To cite this version:

Alexios Christopoulos, Pierre Le Doussal, Denis Bernard, Andrea de Luca. Universal out-of-equilibrium dynamics of 1D critical quantum systems perturbed by noise coupled to energy. Physical Review X, 2023, 13 (1), pp.011043. 10.1103/PhysRevX.13.011043. hal-04245132

HAL Id: hal-04245132

https://hal.science/hal-04245132v1

Submitted on 11 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Universal Out-of-Equilibrium Dynamics of 1D Critical Quantum Systems Perturbed by Noise Coupled to Energy

Alexios Christopoulos,¹ Pierre Le Doussal,² Denis Bernard,² and Andrea De Luca¹

¹Laboratoire de Physique Théorique et Modélisation, CY Cergy Paris Université,
CNRS, F-95302 Cergy-Pontoise, France

²Laboratoire de Physique de l'École Normale Supérieure, CNRS, ENS and PSL University,
Sorbonne Université, Université de Paris, 75005 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 3 January 2022; revised 3 November 2022; accepted 8 December 2022; published 22 March 2023)

We consider critical one-dimensional quantum systems initially prepared in their ground state and perturbed by a smooth noise coupled to the energy density. By using conformal field theory, we deduce a universal description of the out-of-equilibrium dynamics. In particular, the full time-dependent distribution of any two-point chiral correlation function can be obtained from solving two coupled ordinary stochastic differential equations. In contrast to the general expectation of heating, we demonstrate that, over the noise realizations, the system reaches a nontrivial and universal stationary distribution of states characterized by broad tails of physical quantities. As an example, we analyze the entanglement entropy production associated to a given interval of size $\ell$. The corresponding stationary distribution has a 3/2 right tail for all $\ell$ and converges to a one-sided Levy stable for large $\ell$. We obtain a similar result for the local energy density: While its first moment diverges exponentially fast in time, the stationary distribution, which we derive analytically, is symmetric around a negative median and exhibits a fat tail with 3/2 decay exponent. We show that this stationary distribution for the energy density emerges even if the initial state is prepared at finite temperature. Our results are benchmarked via analytical and numerical calculations for a chain of noninteracting spinless fermions with excellent agreement.

DOI: 10.1103/PhysRevX.13.011043

Subject Areas: Condensed Matter Physics
Statistical Physics
Strongly Correlated Materials

## I. INTRODUCTION

The coherent dynamics of macroscopic quantum systems has attracted a lot of interest in the past years, both for its fundamental importance and for its relevance to experimental setups [1–4]. From the theoretical point of view, recent advances have led to a much deeper level of understanding about the mechanism of equilibration and thermalization of isolated many-body quantum systems brought out of equilibrium. In most cases, the eigenstate thermalization hypothesis ensures relaxation to the canonical Gibbs ensemble and the emergence of standard thermodynamics. In recent years, the quest for intriguing out-of-equilibrium phases which escape thermalization has pinpointed a few phenomena, such as many-body localization (MBL) [5], equilibration toward generalized Gibbs ensembles due to integrability [6], quantum scars [7–9], and Hilbert space fragmentation [10,11]. In particular, for integrable systems, anomalous transport has

been observed beyond the expected ballistic, with super-diffusive behavior [12–15]. Nevertheless, weak integrability breakings have been shown to eventually lead to thermalization [16–18].

Recently, the possibility of exploring periodically driven systems has led to a larger class of setups, which culminated in the discovery of discrete time crystalline order [19,20]. The stability of these phases is based either on MBL, which prevents heating to infinite temperature, or on high-frequency expansion, which leads to long-lived prethermal behavior [21].

Stochastic unitary dynamics has recently also attracted a lot of interest, also for the possibility of realizing it in concrete experiments, including cold-atom platforms [22,23] and trapped ions [24–26]. From the theory side, several exact results have recently emerged with discrete time evolution involving random unitary circuits [27–29]. Although finite time dynamics exhibits interesting connections to growth processes, at large time the system relaxes to a trivial infinite temperature ensemble. Other results are obtained in the context of stochastic dynamics in continuous time [30–32], on the Fredrickson-Andersen model [33], and the quantum simple symmetric exclusion process (QSSEP) [34,35]. In these cases, by looking at the average dynamics over noise

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
realizations, one obtains an effective Lindbladian description, which once again admits only the infinite temperature state as a stationary point. Physically, this is the consequence of ergodicity in the unitary group, while, from a mathematical perspective, it follows from the difficulty of having, in general, any nontrivial steady state for a Lindbladian combining a unitary dynamics and jump operators [36]. It is not obvious, however, that the average dynamics is always representative of the typical one. Recent studies, thus, explore the statistical behavior beyond the average dynamics [31], including large deviations [37–40]. In particular, for the finite QSSEP with periodic boundary conditions, it was shown that the stationary state [41] is uniformly distributed among Gaussian states with the same occupancy of the initial state [38]. In the presence of appropriate open boundary conditions, a nontrivial steady state can be attained in the QSSEP [35,42]. Another remarkable mechanism involves the repeated monitoring of local degrees of freedom which compete with the inner unitary dynamics of the system to produce the measurement-induced phase transition, with nontrivial stationary states visible in the statistics of entanglement [43–45]. In general, an important question for quantum stochastic dynamics is whether nontrivial stationary distribution can emerge when the thermodynamic limit is taken before the large-time one. In this direction, the finite time fluctuations are studied for QSSEP in Ref. [46]. From the experimental perspective, the observation of the statistics beyond the average has recently been explored. In general, it is based on the physical realization of multiple copies of the same system undergoing the same stochastic evolution. This is certainly possible in those setups where the noise has a quenched physical origin (e.g., laser speckles in cold atoms [47]). It is more problematic for the stochasticity induced by quantum measurements, where two identical realizations can only be the consequence of postselection [48]. Nevertheless, experimental indications of the measurement-induced phase transition have been recently reported in small systems [26].

In this paper, we consider critical one-dimensional quantum systems initially prepared in their ground state. In practice, the starting Hamiltonian is homogeneous and gapless so that scale invariance is present. In this case, its low-energy spectrum is independent on the microscopic details and is well described by a conformal field theory (CFT). The behavior of quantum systems perturbed by different sources of noises has attracted great interest in the recent years [49–51], in particular, investigating their stability properties under $1/f$ noise [52,53]. Here, we introduce a spatially smooth, white noise in time, coupled at $t > 0$ to the energy density, and bring the system out of equilibrium by evolving it under the corresponding unitary dynamics. By using conformal field theory, we deduce a universal description of the out-of-equilibrium dynamics. We show that the full distribution of correlation functions reaches a nontrivial stationary limit, not visible at the level

of noise averages which instead exhibit apparent heating. Before going into the details, we summarize the model and the results in the next section.

## II. MODEL AND RESULTS

We consider a 1D model at a second-order quantum phase transition, initially prepared in the ground state $|\Psi_0\rangle$ of its gapless Hamiltonian $\hat{H}_0$. To simplify the notation, we indicate the ground state averages as $\langle\dots\rangle = \langle\Psi_0|\dots|\Psi_0\rangle$. Also, we assume continuous space, but the treatment can be readily extended to lattice systems in the scaling limit. At time $t = 0$, a perturbation $\hat{H}_1$ is turned on, by coupling a space-dependent noise term with the system energy density. The total Hamiltonian then takes the form

$$ \hat{H} = \hat{H}_0 + \hat{H}_1 = \int dx[1 + \eta(x, t)]\hat{h}(x), \quad (1) $$

where $\hat{h}(x)$ is the Hamiltonian density and $\eta(x, t)$ is the noise perturbation. Examples of this kind of Hamiltonian with an experimental counterpart can be obtained, for instance, by considering tight-binding noninteracting fermions with a noisy hopping (see Sec. VIII) or antiferromagnetic Heisenberg chains with a stochastic modulation of the coupling, whose out-of-equilibrium properties can be investigated by means of neutron scattering [54]. A related perturbation has been recently considered in Ref. [55], where, however, the noise $\eta(x, t)$ in Eq. (1) is restricted to a single-wavelength modulation, so that the Hamiltonian is spanned by the $SL_2$ finite subgroup of the Virasoro algebra. Here, instead, we choose to characterize it by the space-time correlation $\overline{\eta(x, t)\eta(x', t')}$ = $\delta(t - t')f(x - x')$. The function $f(x)$ controls the noise correlation and has the dimension of a (turnover) time; it is even and has positive Fourier transform, and we choose it smooth, monotonically decreasing for $x > 0$, with a fast decay at infinity when $x \gg 1$. In the following, we use a hat to denote quantum operators $\hat{O}$ and indicate as $\hat{Q}$ the average of any quantity $Q$ over the noise realizations. It is well known that the low-energy behavior of $\hat{H}_0$ can be described in terms of a CFT. In the absence of the noise, the Heisenberg evolution of chiral primary fields under $H_0$ reduces to a translation in time at the light velocity $v$: $\hat{\phi}^{\pm}(y, t) = \hat{\phi}^{\pm}(y \pm vt)$. As shown in Ref. [56], the coupling in Eq. (1) is such that the time evolution of chiral primary fields is still equivalent to a flow but along a random stochastic trajectory [see Eq. (7)]. In Ref. [56], the focus is on quench protocols resulting in an initial state with short spatial correlation length, while here the initial state is gapless with quasi-long-range order. This implies that, for primary fields $\hat{\Phi}(x, t) = \hat{\phi}_+(x, t) \times \hat{\phi}_-(x, t)$ with a scaling dimension $\Delta = \Delta_+ + \Delta_-$, the time evolution of the two-point correlation function can be parametrized in terms of $\kappa^{\pm}$, two random functions of the spatial separation $\ell$ and time $t$:
$$C(\ell, t) = a_0^{2\Delta} \langle \hat{\Phi}(0, t) \hat{\Phi}(\ell, t) \rangle = C(\ell, t=0) e^{-\Delta_+ \kappa^+ - \Delta_- \kappa^-}, \quad (2)$$

where $a_0$ is an ultraviolet cutoff. As we discuss below in Sec. IV, the statistics of the random variable $\kappa^\pm$ can be related to the joint distribution of four stochastic trajectories in the same noisy environment. By studying in detail the resulting stochastic process, we show that the time-dependent distribution of each $\kappa = \kappa^\pm$ can actually be reduced to the solution of two coupled ordinary stochastic differential equations in Eq. (16) below. Quite surprisingly, this leads to a stationary state $t \to \infty$ characterized by broad distributions of the variable $\kappa^\pm$ for all values of the separation $\ell$. Such a distribution takes a simple universal form in the limit of small and large $\ell$, compared to the correlation length of the noise, here chosen of the order of unity.

At large $\ell \gg 1$, we set $\kappa = \theta\ell^2\chi/f(0) + O(\ell)$ and we find that the random variable $\chi$ is distributed according to the stable one-sided Levy distribution of index 1/2:

$$\mathcal{L}(\chi) \equiv \frac{1}{\sqrt{2\pi}} \frac{e^{-1/2\chi}}{\chi^{3/2}} \Theta(\chi). \quad (3)$$

For small $\ell \ll 1$, we set $\kappa^\pm = \ell^2 \tilde{\kappa}_0[(\omega/\omega_0) - 1]$, where $\tilde{\kappa}_0$ and $\omega_0$ are constants expressed in terms of derivatives of $f(x)$ around $x = 0$ [see Eq. (24)]. Then, we obtain the following stationary distribution for the random variable $\omega$:

$$\mathcal{B}(\omega) \equiv \frac{\Gamma(3/4)}{\Gamma(1/4)} \frac{1}{\sqrt{\pi}(1 + \omega^2)^{3/4}}. \quad (4)$$

As explained in Sec. VII, in the small $\ell$ limit, one has the proportionality $\langle T^\pm \rangle \propto \lim_{\ell \to 0} \kappa^\pm / \ell^2$, where $T^\pm$ are the chiral components of the stress-energy tensors. We can, thus, infer the stationary distribution of the local energy density, since $h(x, t) = \langle \hat{h}(x, t) \rangle = v(\langle \hat{T}^+ \rangle + \langle \hat{T}^- \rangle)$. This leads to

$$\lim_{t \to \infty} h(x, t)^{\text{in law}} \sim \frac{v c \tilde{\kappa}_0}{4\pi} (\Omega/\omega_0 - 2), \quad (5)$$

where $c$ is the central charge of the underlying CFT and $\Omega = \omega^+ + \omega^-$. At large times, one expects the two chiral components to be only weakly correlated, as they depend from more and more distant regions. Thus, $\omega^+$ and $\omega^-$ can be assumed to be two independent random variables both distributed according to $\mathcal{B}(\omega)$. Thus, the stationary distribution of the local energy density reaches a universal form, still displaying a fat tail with 3/2 decay exponent and no finite integer moments. This is consistent with our finding in Eq. (40), where we compute the first moment at all time and show that, as a consequence of conformal invariance, it diverges exponentially fast in time. This is a clear

manifestation that the average over noise is not indicative of the typical behavior of the full distribution. Interestingly enough, we are able to extend such a result for the steady distribution of the local energy density even if the initial state is prepared at a finite temperature. Remarkably, the initial temperature affects only the transient dynamics, while the stationary distribution remains unchanged. As a consequence, the 3/2 tail is robust even at finite temperature.

Equations (3) and (4) characterize asymptotic behavior of the stationary distribution for $\kappa$ at small and large $\ell$, respectively. For intermediate values of $\ell$, a stationary distribution for $\kappa^\pm$ is still reached, but its explicit form depends on the function $f(x)$. Nonetheless, we prove that the 3/2 exponent for the right tail is always present for any $\ell$ and a sufficiently smooth $f(x)$ (see Appendix E 4). In Appendix E 7, we provide a general method to find the stationary distribution, together with analytical formula for some specific solvable choices for $f(x)$.

A direct manifestation of these results can be seen in bipartite entanglement Rényi entropies. Indeed, one can express entanglement entropies in terms of correlators of twist fields [57], which leads to the relation (see Sec. VI)

$$S_t^{(n)} = \frac{1}{1-n} \ln \text{Tr}[\rho_{\ell,t}^n] = S_{t=0}^{(n)} + \frac{(n+1)c}{24n}(\kappa^+ + \kappa^-), \quad (6)$$

where $\rho_{\ell,t}$ is the reduced density matrix at time $t$ of an interval of length $\ell$ and $S_{t=0}^{(n)}$ is the initial entropy of the ground state. Assuming that at large times $\kappa^+$ and $\kappa^-$ are only weakly correlated, the corresponding stationary distribution can be extracted by convolution. In general, the 3/2 right tail emerges for all sizes of the interval $\ell$.

To test the validity of the field theory description, in Sec. VIII, we study analytically and numerically a chain of noninteracting spinless fermions at half filling with a noisy hopping. At low energy, they are well described by Dirac fermions corresponding to a $c = 1$ CFT. We identify a scaling limit where the noise correlation length on the lattice $\xi$ diverges and the CFT predictions are exactly recovered, as confirmed numerically by computing the local energy and entanglement entropy on the lattice. We make use of this lattice model also to estimate the timescales at which the CFT framework breaks down. In practice, in the fermionic formulation, the initial ground state corresponds to a filled Fermi sea up to the Fermi momentum $|k| < k_F$. Because of the noise, the sharp edges start to effectively broaden and smoothen. We can, thus, estimate the regime of validity of the CFT prediction as the timescale at which the scale of broadening of each Fermi point $\delta k_F$ becomes comparable to their relative distance $\delta k_F \sim k_F$. As shown in Sec. VIII C, the precise calculation leads to the timescale $\tau^* \sim \xi \ln \xi$, such that, at times much larger than $\tau^*$, significant deviations from CFT are to be expected. Such a prediction is qualitatively confirmed in our numerics.
### III. SOLUTION OF THE DYNAMICS VIA CFT

In this section, we summarize some results of Ref. [56] which allow computing dynamical correlation functions in the presence of the noise term (1) within the framework of CFT. Beyond its application to equilibrium systems, conformal symmetry [58] has been successfully employed even in quantum quenches and general out-of-equilibrium dynamics [59,60]. It implies that all local operators can be split into two chiral [+(−) for right (left) moving] components organized into families [61]. All operators within each chiral family descend from a primary field $\hat{\phi}^{\pm}(x)$ with given conformal dimension $\Delta^{\pm}$. In particular, the Hamiltonian density can be represented as $\hat{h}(x) \sim v[\hat{T}^{+}(x) + \hat{T}^{-}(x)]$, where $v$ is the light velocity and $\hat{T}^{+}$ and $\hat{T}^{-}$ are the two components of the stress-energy tensor, with $\hat{T}^{+}(x) - \hat{T}^{-}(x)$ the momentum density. This implies that, under $\hat{H}_0$, chiral primary fields simply translate in time $\hat{\phi}^{\pm}(x, t) = \hat{\phi}(x \mp vt)$. Although the Hamiltonian (1) manifestly breaks Lorentz invariance, it is still possible to use the machinery of CFT to investigate the corresponding dynamics. In Ref. [56], it is shown that the time evolution of primary fields under $\hat{H}$ can be interpreted as a conformal transformation. In practice, one introduces the stochastic trajectories $q^{\pm}(s)$ as a solution of the Langevin equation

$$ \frac{dq^{\pm}(s)}{ds} = \pm v \{1 + \eta[q^{\pm}(s), s]\}, \quad (7) $$

where the Ito convention is assumed. Then, we define the functions $X_t^{\pm}(y)$ as the initial condition for Eq. (7) at $t=0$ [i.e., $q^{\pm}(0) = X_t^{\pm}(y)$] such that $q^{\pm}(t) = y$. It follows that the evolution of a primary field is simply given by

$$ \hat{\phi}^{\pm}(y, t) = [X_t^{\pm}(y)]^{\pm} \hat{\phi}(X_t^{\pm}(y), 0), \quad (8) $$

and arbitrary correlation functions at time $t$ can be reduced to those in the initial state via

$$ \left\langle \prod_{i=1}^{n} \hat{\phi}_{i}^{+}(y_{i}, t) \right\rangle = \mathcal{J}(y_{1}, \dots, y_{n}) \left\langle \prod_{i=1}^{n} \hat{\phi}_{i}^{+}(X_{t}^{+}(y_{i})) \right\rangle, \quad (9) $$

where the factor $\mathcal{J}(y_1, \dots, y_n) = \prod_i [X_t^{\pm}(y_i)]^{\Delta_i^{\pm}}$ accounts for the Jacobian of the conformal transformation. An analogous equation can be written for the other chiral component with $X_t^{\pm} \to X_t^{\mp}$ and $\Delta_i^{\pm} \to \Delta_i^{\mp}$.

To study the correlation functions in Eq. (9) and their sample to sample fluctuations, we thus need the joint probability distribution function (JPDF) of the set of $2n$ random variables $X_t^{\pm}(y_j)$, $j=1, \dots, n$. Let us first focus on the JPDF of $X_t^{\pm}(y_j) = x_j$ for a fixed chirality, choosing either $\pm$, which we denote $P_t^{\pm}(\mathbf{x}|\mathbf{y})$. Given $n$ trajectories $q_j^{\pm}(s)$ satisfying Eq. (7) with end points $q_j^{\pm}(0) = x_j$ and

$q_j^{\pm}(t) = y_j$, $P_t^{\pm}(\mathbf{x}|\mathbf{y})$ is, thus, the JPDF of the initial positions $\mathbf{x} = (x_1, \dots, x_n)$ of these $n$ trajectories conditioned to the positions of their final point $\mathbf{y} = (y_1, \dots, y_n)$. As shown in Appendix A, it satisfies the Fokker-Planck (FP) equation also studied in the context of turbulence and passive scalar [62–64]:

$$ \partial_t P_t^{\pm}(\mathbf{x}|\mathbf{y}) = \left( \pm v \sum_{i=1}^{n} \partial_i + \frac{v^2}{2} \sum_{i,j=1}^{n} \partial_i \partial_j f(x_i - x_j) \right) P_t^{\pm}(\mathbf{x}|\mathbf{y}), \quad (10) $$

where $\partial_i = \partial/\partial x_i$. Since all trajectories are evolving according to Eq. (7) within the same realization of the noise, they are correlated, which appears as an interaction in Eq. (10). The martingale property from the initial time implies $\bar{x}_i = y_i \mp vt$; additionally, the trajectories cannot cross one another, so that the coordinates $\mathbf{y}$ and $\mathbf{x}$ are always ordered in the same way.

### IV. TWO-POINT CORRELATIONS

Consider a primary field $\hat{\Phi}(x, t) = \hat{\phi}_{+}(x, t) \times \hat{\phi}_{-}(x, t)$ with a scaling dimension $\Delta = \Delta_{+} + \Delta_{-}$. We are interested in the two-point correlator $C(y_1, y_2, t) = a_0^{2\Delta} \langle \hat{\Phi}(y_1, t) \hat{\Phi}(y_2, t) \rangle$. At $t=0$, conformal invariance implies

$$ C(x_1, x_2, t=0) = \left( \frac{a_0}{x_1 - x_2} \right)^{2(\Delta_+ + \Delta_-)}. \quad (11) $$

Making use of Eq. (9), we can relate the correlator at time $t$ with the one at initial time, which leads to

$$ C(y_1, y_2, t) = C(y_1, y_2, t=0) e^{-\Delta_+ \kappa^+ - \Delta_- \kappa^-}, \quad (12) $$

where we assume $y_1 > y_2$ and we set

$$ \kappa^{\pm}(y_1, y_2, t) = \ln \left| \frac{[X_t^{\pm}(y_1) - X_t^{\pm}(y_2)]^2}{(y_1 - y_2)^2 X_t^{\pm}(y_1) X_t^{\pm}(y_2)} \right|. \quad (13) $$

This expression gives access to the full statistics of the correlation function. We first focus on either $\kappa^+$ or $\kappa^-$. Indeed, the trajectories corresponding to the two chiral components are typically separated by a distance of approximately $2vt$. Therefore, at large time for $2vt \gg 1$, the noises they feel become uncorrelated, and it is reasonable to expect the two components to have little statistical correlation. Let us define $\ell \equiv y_1 - y_2 > 0$ and the dimensionless ratio $r \equiv [X_t^{\pm}(y_1) - X_t^{\pm}(y_2)]/\ell$. Using spatial homogeneity, the one-point PDF of $\kappa = \kappa^{\pm}(y_1, y_2, t)$ can be only a function of $\ell$ and $t$ and must be independent on the chirality. Although this PDF does not obey a closed equation, we can derive a FP equation for the JPDF $P_t(r, \kappa)$ of $\kappa$ and $r$. Here, we sketch the procedure, leaving the
technical details to Appendix B. We first perform a point splitting in the derivatives in Eq. (13):

$$X_{t}^{\pm}(y) = \lim_{\epsilon \to 0} \frac{X_{t}^{\pm}(y+\epsilon) - X_{t}^{\pm}(y)}{\epsilon}. \quad (14)$$

Therefore, in the limit $\epsilon \to 0$, both $\kappa^{\pm}(y_1, y_2, t)$ and $r$ depend only on the initial points $x_1, x_2, x_3$, and $x_4$ of four trajectories ending, respectively, at $y_1, y_2, y_1 + \epsilon$, and $y_2 + \epsilon$, whose JPDF is given by Eq. (10) for $n=4$ (see Fig. 1 for a sketch). Because of homogeneity, the center of mass $x_1 + x_2$ decouples and one obtains, in the limit $\epsilon \to 0$, a FP equation for the JPDF involving only the three variables

$$r = r(y_1, y_2, t) \equiv \frac{X_t(y_1) - X_t(y_2)}{y_1 - y_2} = \frac{x_1 - x_2}{\ell}, \quad (15)$$

$x'_1 = (x_3 - x_1)/\epsilon$, and $x'_2 = (x_4 - x_2)/\epsilon$. Remarkably, performing another change of variable, one finds that the random variables $r$ and $\kappa = \log[r^2/(x'_1x'_2)]$ satisfy a stand-alone FP equation. Instead of working directly with this equation, it is more convenient to reformulate the problem in terms of an equivalent stochastic differential equation (SDE) in Ito's form (see Appendix B)

$$dr = v dW_1(t), \quad d\kappa = v^2 g(r) dt + v dW_2(t), \quad (16)$$

where $v^2g(r)$ is a drift term and $dW_1(t)$ and $dW_2(t)$ are two centered Gaussian white noises in time of $r$-dependent variances $\overline{dW_1(t)^2} = 2A(r)dt$ and $\overline{dW_2(t)^2} = 2C(r)dt$ and cross-correlation $\overline{dW_1(t)dW_2(t)} = B(r)dt$ [65]. To simplify the notation, we introduce

<figure><img src="image_3.png" /><figcaption>FIG. 1. Sketch of the four trajectories required to access the full PDF of two-point correlation functions in the limit $\epsilon \to 0$.</figcaption></figure>

$$\begin{align}
A_{\ell}(r) &= \frac{f(0) - f(\ell r)}{\ell^2}, \nonumber \\
B_{\ell}(r) &= \frac{2f'(\ell r)}{\ell} + \frac{4[f(0) - f(\ell r)]}{\ell^2 r}, \nonumber \\
C_{\ell}(r) &= \frac{4[f(0) - f(\ell r) + \ell r f'(\ell r)]}{\ell^2 r^2} - f''(0) - f''(\ell r), \nonumber \\
g_{\ell}(r) &= -f''(0) - \frac{2[f(0) - f(\ell r)]}{\ell^2 r^2}. \tag{17}
\end{align}$$

Equations (16) must be solved with the initial condition $r=1$ and $\kappa=0$ at $t=0$. A remark is in order: By construction, Eqs. (16) [and Eq. (10)] are derived by adding a time slice $dt$ in the backward time direction (Fig. 1). This poses no problems if one is interested in the distribution of $\kappa$ in Eq. (12) at a fixed time $t$, and this is the focus in this work. However, the backward and forward evolutions lead to different stochastic processes in $t$ [66]. In particular, correlations in time of the forward process (13), e.g., $\kappa^{\pm}(y_1, y_2, t)\kappa^{\pm}(y'_1, y'_2, t')$, cannot be obtained from Eqs. (16).

Since the equation for $r$ does not involve $\kappa$, one may first solve for $r(t)$ and then insert the solution for $r(t)$ in the equation for $\kappa(t)$. For finite $\ell$, Eq. (16) cannot be solved explicitly for an arbitrary $f(x)$. Nevertheless, one can understand its behavior at finite time in two regimes $\ell = y_1 - y_2 \gg 1$ and $\ell = y_1 - y_2 \ll 1$, as well as in the large-time limit.

### A. Small separation $\ell \ll 1$

For small $\ell \ll 1$, we can Taylor expand the function $f$ in Eq. (17). One finds the leading behavior at small $\ell$ of each function as

$$A_{\ell}(r) \approx -\frac{1}{2}f^{(2)}(0)r^2, \quad B_{\ell}(r) \approx \frac{\ell^2}{6}f^{(4)}(0)r^3, \quad (18)$$

$$C_{\ell}(r) \approx -\frac{\ell^4}{72} f^{(6)}(0) r^4, \quad g_{\ell}(r) \approx \frac{\ell^2}{12} r^2 f^{(4)}(0). \quad (19)$$

The positivity of the Fourier transform of $f(x)$ implies alternating signs of the even derivatives $\operatorname{sgn}[f^{(2n)}(0)] = (-1)^n$, and, thus, all variables above are positive. We see that Eq. (16) can be rewritten, after a redefinition of the noises $dW_1(t) = r dB_1(t)$ and $dW_2(t) = -(\ell^2/6)r^2 dB_2(t)$, in the form

$$dr = r v dB_1(t), \quad d\kappa = \frac{\ell^2 r^2 v}{6} \left( \frac{v f^{(4)}(0) dt}{2} - dB_2(t) \right), \quad (20)$$

where now $dB_1(t)$ and $dB_2(t)$ are $r$-independent Gaussian white noises with fixed correlation matrix $\overline{dB_1(t)^2} = -f''(0)dt$, $\overline{dB_2(t)^2} = -f^{(6)}(0)dt$, and $\overline{dB_1(t)dB_2(t)} = -f^{(4)}(0)dt$.
Let us first discuss the marginal distribution $P_t(r)$ of $r$ in this regime $\ell \ll 1$. One can solve the stochastic equation for $r$ and obtain after an application of Ito's lemma:

$$r(t) = e^{-\theta t + v B_1(t)}. \qquad (21)$$

We define $\theta = -v^2 f''(0)/2 > 0$, which is the inverse of a characteristic time. Alternatively, one can change the variable to $\rho = \log r$, which obeys the stochastic equation $d\rho = -\theta dt + v dB_1(t)$, implying that $\rho(t)$ is a Brownian motion with a drift. Hence, $P_t(r)$ is a log-normal distribution for the variable $r$, with average and variance

$$\overline{\ln r} = -\theta t, \quad \mathrm{Var}[\ln r] = 2\theta t. \qquad (22)$$

Since $\theta > 0$, this shows, interestingly, that the trajectories $X_t(y_1)$ and $X_t(y_2)$ tend to get closer as time grows, a manifestation of the phenomenon of sticky particles observed in turbulent fluids [68–71]. On the other hand, $\bar{r} = 1$ holds independently of time, which shows that, although the typical value $r_{\text{typ}} = e^{\overline{\log r}}$ decreases to zero, the distribution of $r$ is necessarily broadening with time. Hence, higher moments grow with time, and, within the small $\ell$ approximation (21), one has $\bar{r}^{\overline{\ell}} \sim e^{n(n-1)\theta t}$. Using this result, it is easy to calculate the noise average of $\kappa = \kappa^{\pm}$ by simply averaging Eq. (20) using that $\overline{dB_2} = 0$. Integrating over time, one finds

$$\bar{\kappa} = \frac{v^2 f^{(4)}(0) \ell^2}{24\theta} (e^{2\theta t} - 1) + O(\ell^4). \qquad (23)$$

This result allows one to obtain the noise average of $\ln C(y_1, y_2, t)$, irrespective of the possible correlations between $\kappa^+$ and $\kappa^-$, by averaging the logarithm of Eq. (12). It is possible to calculate the higher integer moments $\kappa^{\overline{\ell}}$, and one finds (see Appendix F) that they all grow exponentially with time within the validity of the small $\ell$ regime. However, upon averaging Eq. (16) over the noise, we observe that $\bar{\kappa} \le 2\theta t$ is an exact bound, since $g(r) \le g(\infty) = -f''(0)$. Thus, while the moments are still diverging at large time, the fast exponential growth is valid only when $\ell^2 e^{2\theta t} \lesssim 1$.

Nevertheless, as we now show, the PDF of $\kappa$ remarkably converges to a stationary distribution at large time. This distribution is very broad and consistently does not possess any integer moments for $n \ge 1$. To obtain the PDF of $\kappa$, we proceed in two steps: First, we plug the solution (21) into the equation (20) for $\kappa$; second, we reparametrize the Wiener processes $B_1(t)$ and $B_2(t)$ by looking at them from their final point, i.e., setting $\tilde{B}_{1,2}(s) = B_{1,2}(t) - B_{1,2}(t-s)$ (see Appendix D for details). This allows us to recast the resulting stochastic equation in the standard form studied by Bougerol [72,73]. This equation is best expressed with a change of variable $\kappa = \ell^2 \tilde{\kappa}_0(\omega/\omega_0 - 1)$ with

$$\begin{aligned} \tilde{\kappa}_0 &= -\frac{1}{12} \frac{f^{(4)}(0)}{f''(0)}, \\ \omega_0 &= \frac{f^{(4)}(0)}{\sqrt{f^{(6)}(0)f^{(2)}(0) - f^{(4)}(0)^2}}, \end{aligned} \qquad (24)$$

which leads to the SDE for the variable $\omega$ in the form

$$d\omega = 2\theta\omega dt + \sqrt{8\theta} \sqrt{1 + \omega^2} d\tilde{B}, \quad \omega(t=0) = \omega_0, \qquad (25)$$

where $\tilde{B}(t)$ is a single standard Wiener process $[d\tilde{B}(t)^2 = dt]$ expressed as a linear combination of $\tilde{B}_{1,2}(t)$. Note that $\omega_0$ is always real and positive because of Cauchy-Schwarz inequality [see Eq. (D12)]. In the large-time limit, the variable $\omega$ attains the stationary distribution Eq. (4). To see this explicitly, we perform another change of variable $\omega = \sinh(Y)$, so that the random variable $Y$ satisfies

$$dY = -2\theta \tanh Y dt + \sqrt{8\theta} d\tilde{B}, \quad \sinh[Y(0)] = \omega_0. \qquad (26)$$

Equation (26) describes the Langevin motion of a particle in a confining potential $U(y) = 2\theta \log \cosh y \approx 2\theta|y|$ at temperature $4\theta$. Hence, it reaches an (equilibrium) stationary measure at large time, $P_{\text{stat}}(Y) = C/\sqrt{\cosh(Y)}$ with $C = \sqrt{2\pi}/\Gamma(1/4)^2$, which corresponds to the stationary distribution (4) for the variable $\omega$. Consistently, the power-law tail $\propto |\omega|^{-3/2}$ implies that $\omega$ (and $\kappa$ as well) does not have finite integer moments.

We remark that, because of the reparametrization $B_j(t) \to \tilde{B}_j(s)$, the stochastic process (25) for $\omega$ and the original one for $\kappa$ in Eq. (20) are not equivalent: The former is ergodic in time, while the latter has a finite (random) limit $\kappa(t \to \infty)$. Nevertheless, they are equivalent in law at fixed $t$ (see also the discussion in Appendix D).

## B. Large interval $\ell \gg 1$ at short times

The leading behavior at large $\ell$ of each function up to $O(1/\ell^2)$ reads

$$\begin{aligned} A_{\ell}(r) &\approx \frac{f(0)}{\ell^2}, & B_{\ell}(r) &\approx 4 \frac{f(0)}{\ell^2 r}, \\ C_{\ell}(r) &\approx -f''(0) + 4 \frac{f(0)}{\ell^2 r^2}, & g_{\ell}(r) &\approx -f''(0) - \frac{2f(0)}{\ell^2 r^2}, \end{aligned} \qquad (27)$$

where we assume that $f(x)$ decays faster than a power law. At leading order, one can set $B_{\ell}(r) \approx 0$, which implies that the equation (16) for $\kappa$ becomes independent of $r$. Using $g_{\ell}(r) \approx -f''(0)$ and $C_{\ell}(r) \approx -f''(0)$ to leading order, one obtains

$$\kappa = 2\theta t + \sqrt{4\theta W(t)}, \quad r = 1 + \frac{\sqrt{2f(0)v}}{\ell} \tilde{W}(t), \qquad (28)$$
where $W(t)$ and $\tilde{W}(t)$ are two uncorrelated standard Brownian motions with $\overline{dW(t)^2} = \overline{d\tilde{W}(t)^2} = dt$. Note that the growth of $\bar{\kappa}$ saturates the exact bound $\bar{\kappa} < 2\theta t$. To obtain this result, we assume not only that $\ell \gg 1$, but also that $\ell r \gg 1$. Although this condition holds for finite time, since $r$ is undergoing diffusion, we see from Eq. (28) that, for $t \sim \ell^2/f(0)$, $r(t)$ may become close to zero and the condition is violated. So, this large $\ell$ expansion holds only at short times. In order to access the large-time limit, we need a different approach that we present in the next section.

## V. LARGE-TIME LIMIT FOR ARBITRARY $\ell$

To search for a stationary distribution for any $\ell$, we derive in Appendix E an evolution equation for the characteristic function of $\kappa$, $Q_k(r_0, t) = e^{-ik\bar{\kappa}r_0}$, where the superscript $r_0 = r(t=0)$ indicates the initial condition for the variable $r$ in Eq. (16), ultimately setting $r_0 = 1$. Looking for a time-dependent solution in the large-time limit $Q_k(r_0, t) \to Q_k(r_0)$, it is useful to perform a shift by extracting an exponential factor $Q_k(r) = e^{ik\kappa_0(\ell r)} G_k(r)$, with

$$ \kappa_0(x) = -\log\left(\frac{2[f(x) - f(0)]}{x^2 f''(0)}\right). \qquad (29) $$

In this way, the function $G_k(r)$ satisfies the Schrödinger-like equation for $r \ge 0$

$$ -G_k''(r) - k(k + i)V(r)G_k(r) = 0 \qquad (30) $$

with boundary conditions $G_k(0) = 1$ and $\lim_{r \to +\infty} G_k(r) = 0$. The potential takes the form

$$ V(r) = -\frac{d^2}{dr^2}\log[f(0) - f(\ell r)] + \frac{\ell^2 f''(0)}{f(0) - f(\ell r)}. \qquad (31) $$

Studying this equation (see Appendix E), we find that the stationary distribution $P_{\text{stat}}(\kappa)$ [obtained by Fourier inversion of $Q_k(r=1)$] depends nontrivially on $\ell$ and $f(x)$. Additionally, the fact that $G_k(r)$ depends on $k$ only via the combination $k(k+i)$ implies the exact symmetry relation

$$ \frac{P_{\text{stat}}[-\kappa_0(\ell) + u]}{P_{\text{stat}}[-\kappa_0(\ell) - u]} = e^u, \qquad (32) $$

reminiscent of a fluctuation theorem [74]. In the limit of small $\ell$, one has $\kappa_0(\ell) \to \ell^2 \bar{\kappa}_0$, and as a consistency check, in Appendix E 5, we reobtain the asymptotic solution at small $\ell$ (4) derived in Sec. IVA in a completely different way. Note that in this limit Eq. (32) reduces to $B(\omega) = B(-\omega)$ in Eq. (4).

In the opposite limit of large $\ell$, one observes that the potential approaches a constant value $V(r) \to \ell^2 f''(0)/f(0)$, which leads to the simple solution of Eq. (30) in the form $G_k(r) \sim e^{-\sqrt{2k(k+i)\theta/f(0)}r\ell}$. As a consequence, setting $\kappa = \theta\ell^2\chi/f(0)$, one finds that, in the limit $\ell \to \infty$, the variable $\chi$ is distributed according to

$$ \mathcal{L}(\chi) \equiv \frac{1}{\sqrt{2\pi}} \frac{e^{-1/2\chi}}{\chi^{3/2}} \Theta(\chi), \qquad (33) $$

i.e., to the stable one-sided Levy distribution of index $1/2$. Note that in this limit the shift $\kappa_0(\ell) = O[\log(\ell)]$ and is negligible with respect to $\kappa = O(\ell^2)$. Additionally, the symmetry Eq. (32) is consistent with $\mathcal{L}(\chi < 0) = 0$.

Equations (4) and (33) characterize asymptotic behavior of the stationary distribution for $\kappa$ at small and large $\ell$, respectively. For intermediate values of $\ell$, an explicit expression is not available for generic $f(x)$. However, we provide some cases which are analytically solvable for any $\ell$, e.g., for $f(x) = (1/\cosh x)^n$ with $n=1, 2$.

In Appendix E 4, we prove that, for sufficiently smooth $f(x)$, $Q_k(r) = 1 + O(\sqrt{k})$ at small $k$, irrespectively of $\ell$. This translates onto a $-3/2$ power-law tail on the positive $\kappa$ side. For the left tail $\kappa \to -\infty$, it follows from Eq. (32) that $P_{\text{stat}}(\kappa)$ decays exponentially for any finite $\ell > 0$. An intuition on why the $3/2$ exponent appears is as follows: The random variable $r$ is attracted to $r=0$ since, for $r < 1$, $\log r$ is approximately Brownian with a negative drift [see Eq. (22)]. In this regime, $\kappa$ remains almost constant as $dk \propto r^2 dt$. However, $r$ starts from 1 and has a finite probability to move right toward $r > 1$; in this case, $\kappa$ increases by $\sim \Delta t$ [see Eq. (28)], where $\Delta t$ is the time interval spent before $r(t)$ hits again 1. Since in this regime $r(t)$ is approximately an unbiased Brownian, $\Delta t$ is the corresponding first-passage time, distributed as $1/(\Delta t)^{3/2}$.

Finally, let us recall that for time $t \gg 1/(2v)$ the two chiral components $\kappa^+$ and $\kappa^-$ are expected to decorrelate; hence, their joint distribution reaches in the large-time limit the factorized form $P_{\text{stat}}(\kappa^+)P_{\text{stat}}(\kappa^-)$. We expect that $\kappa^\pm(t)$ defined in Eq. (13) as a stochastic process in time has $P_{\text{stat}}(\kappa^\pm)$ as the one-time ergodic measure. As we already clarified, this is not in contradiction with the fact that $\kappa$ defined in the auxiliary stochastic system (16) has a (random) finite limit $\kappa^\pm(t \to \infty)$ (see also the remark in Appendix D).

## VI. ENTANGLEMENT PRODUCTION

An interesting application of the above results is the calculation of the entanglement entropies. Let us define as $\rho_{A,t}$ the reduced density matrix for the interval $A = [y_1, y_2]$ at time $t$. Then, the Rényi entropies are defined as $S_t^{(n)} = 1/(1-n) \ln \text{Tr} \rho_{A,t}^n$. Introducing the twist fields $\hat{\Phi}_n(y,t)$ [57], one can identify
$\mathrm{Tr} \rho_{A,t}^n \propto \langle \hat{\Phi}_n(y_1, t) \hat{\Phi}_n(y_2, t) \rangle$. As anticipated in Eq. (6), we can, thus, express the entropy production $S_t^{(n)} = S_t^{(n)} - S_{t=0}^{(n)}$ using Eq. (12) as

$$S_t^{(n)} = \frac{(n+1)c}{24n}(\kappa^+ + \kappa^-), \quad (34)$$

where we use that for twist fields $\Delta^\pm = c(n-1/n)/24$. The von Neumann entropy corresponds to $n=1$, and we denote it simply as $S_t$. Equation (34) shows that all the Rényi entropies are controlled by the same random variable and that the details of the model enter only in the prefactor via the central charge. From Ref. [75], it implies that the entanglement spectrum retains the same form in each noise realization. In other words, denoting as $\lambda_{i,t}$ the eigenvalues of $\rho_{A,t}$ and as $\lambda_{\max,t}$ their maximum, the density of $\log \lambda/\log \lambda_{\max,t}$ is independent of the noise, while $\log \lambda_{\max,t} = -c/24(\kappa^+ + \kappa^-) + O(1)$. Using that the noise average $\overline{S_t^{(n)}} = [(n+1)c/12n]\bar{\kappa}$, we see that at early times two different growth regimes exist: For large intervals ($\ell \gg 1$) the average entropy production grows linearly with time as in Eq. (28), while for small intervals it grows as in Eq. (23). Finally, at large time, we find that the entropy production (34) reaches a stationary distribution given up to a scale by the convolution $P_{\mathrm{stat}} * P_{\mathrm{stat}}$ determined above, still with a $-3/2$ power-law tail and no finite integer moments.

## VII. DISTRIBUTION OF THE ENERGY DENSITY

An interesting quantity to look at is the dynamics of the energy density, which, in the CFT mapping, is encoded in the stress-energy tensor $\hat{T}^+(x)$ and $\hat{T}^-(x)$. Their time evolution can be obtained by direct calculation of the commutators $[\hat{H}, \hat{T}^+(x)]$ in the Heisenberg equation (see Appendix F). Alternatively, one can use the fact that time evolution can be seen as a conformal mapping, and the corresponding transformation of the stress energy reads [61] (for simplicity, we omit the $\pm$ superscript)

$$\hat{T}(y, t) = |X_t'(y)|^2 \hat{T}(X_t(y), t=0) - \frac{c}{24\pi} (S \cdot X_t)(y), \quad (35)$$

where the second term is proportional to the central change $c$ and the Schwarzian derivative

$$(S \cdot X_t)(y) = \frac{2X'''(y)X'(y) - 3X''(y)^2}{2X'(y)^2}. \quad (36)$$

Since in the initial ground state $\langle \hat{T}(y, 0) \rangle = 0$, one has

$$\langle \hat{T}(y, t) \rangle = - \frac{c}{24\pi} (S \cdot X_t)(y). \quad (37)$$

The time evolution of this quantity can be simply obtained from the above analysis of $\kappa^\pm$ in the regime $\ell \to 0$. Indeed, expanding Eq. (13) in small $y_1 - y_2$, we see that

$$\lim_{\ell \to 0} \frac{1}{\ell^2} \kappa(y_2 + \ell, y_2, t) = -\frac{1}{6} (S \cdot X_t)(y_2). \quad (38)$$

Equation (38) is well known in the CFT context: It reflects the occurrence of the conformal anomaly in the transport of the stress tensor in a Gaussian free field theory; see Ref. [61]. Here, it implies that the distribution of $\langle \hat{T}^\pm(y, t) \rangle$ can be obtained from the one of $\kappa^\pm$ in the limit of small $\ell$. Following the derivation in Sec. IVA, one can express the full distribution over the noise as

$$\langle \hat{T}(y, t) \rangle \stackrel{\text{in law}}{=} \frac{c \bar{\kappa}_0}{4\pi} \left( \frac{\omega}{\omega_0} - 1 \right), \quad (39)$$

where the random variable $\omega$ is a solution of Eq. (25). We observe a remarkable identification between the entanglement of an infinitesimal interval and the local energy density. Indeed, denoting the time-dependent local energy density as the quantum expectation $h(x, t) = \langle \hat{h}(x, t) \rangle = v(\langle \hat{T}^+ \rangle + \langle \hat{T}^- \rangle)$, one has $h(x, t) = \lim_{\ell \to 0} v S_t / (\pi \ell^2)$. Additionally, defining the noise average (which is space independent) as $\mathbf{e}(t) = \overline{h(x, t)}$, one has from Eq. (23)

$$\mathbf{e}(t) = \frac{c v^3 f^{(4)}(0)}{48\pi\theta} (e^{2\theta t} - 1). \quad (40)$$

For $t \gg 1/(2v)$, we expect $\langle \hat{T}^\pm \rangle$ to decorrelate; thus, from Eq. (24), we extract the stationary distribution of the one-point energy density:

$$\lim_{t \to \infty} h(x, t) \stackrel{\text{in law}}{\sim} \frac{v c \bar{\kappa}_0}{4\pi} (\Omega/\omega_0 - 2), \quad (41)$$

where $\Omega = \omega^+ + \omega^-$ and $\omega^\pm$ are independent random variables both distributed according to $\mathcal{B}(\omega)$ in Eq. (4).

Surprisingly enough, these considerations can be extended to an initial state prepared at a finite inverse temperature $\beta$. Let us define $\langle \hat{O} \rangle_\beta = \mathrm{Tr}[\hat{O} e^{-\beta H_0}]/Z_\beta$ with the partition function $Z_\beta = \mathrm{Tr}[e^{-\beta H_0}]$. In this case, the quantum expectation of the stress-energy tensor in the initial state takes the form $\langle \hat{T}(y, 0) \rangle_\beta = \pi c / (12\beta^2 v^2)$ [76,77]. Starting from Eq. (35), we can employ Eqs. (38), (14), and (15) to express the random variable $T_\beta = \langle \hat{T}(y, t) \rangle$ in terms of $r$ and $\kappa$ in the limit of small $\ell$ as

$$T_\beta = \lim_{\ell \to 0} \left[ \frac{\pi c r^2}{12\beta^2 v^2} + \frac{c}{4\pi} \frac{\kappa(y, y + \ell, t)}{\ell^2} \right]. \quad (42)$$

Following a procedure similar to Sec. IVA, in Appendix H, we show that Eq. (39) retains its validity even for a finite temperature initial state, where the random variable $\omega$ is still a solution of Eq. (25) but the inverse temperature $\beta$ modifies the initial condition as
$$ \omega_{\beta}(t=0) = \omega_0 + \frac{\pi^2 \omega_0}{3\beta^2 \bar{\kappa}_0 v^2}. \quad (43) $$

It thus follows that, after an initial transient, the distribution of the stress-energy tensor reaches the same stationary limit irrespectively of the inverse temperature $\beta$.

## VIII. COMPARISON WITH FREE FERMIONS

### A. Tight-binding model

We consider a model of noninteracting spinless fermions, where we denote $\tau$ the time in the lattice model with Hamiltonian

$$ \hat{H}_F = \sum_i [1 + \eta_i(\tau)] \hat{h}_i, \\ \hat{h}_i = -J(\hat{a}_i^\dagger \hat{a}_{i+1} + \hat{a}_{i+1}^\dagger \hat{a}_i). \quad (44) $$

Note that the local energy density is defined only up to a total derivative. Such an ambiguity could, in principle, spoil the identification with the continuum limit. A good way to fix it is to impose local parity invariance, which is verified in Eq. (44) [78–80]. On the contrary, in the presence of a finite chemical potential, a slightly different definition is required (see Appendix G). In the absence of noise, this corresponds to a dispersion relation $\epsilon(k) = -2J \cos(k)$. For the noise, we choose the correlation

$$ \overline{\eta_i(\tau)\eta_j(\tau')} = \tau_0 \delta(\tau - \tau') F(i-j), \quad F(j) = f(j/\xi), \quad (45) $$

where $\xi$ is the characteristic correlation length of the discrete model and in the numerics we choose $f(x) = 1/\cosh(x)$, as it corresponds to an analytically solvable case in the CFT limit. The dynamics induced by Eq. (44) is better studied in terms of the noise-averaged Wigner function $n_{\tau}(k) = \sum_{j'} \langle \hat{a}_{j+j'}^{\dagger} \hat{a}_{j} \rangle_{\tau} e^{ikj'}$, where $\langle \cdots \rangle_{\tau}$ denotes the quantum average at time $\tau$ under the evolution $\hat{H}_F$ in Eq. (44). We choose the initial state as the ground state so that $n_{\tau}(k)$ does not depend on the lattice site $j$. Also, $n_{\tau=0}(k) = \Theta(k + k_F) - \Theta(k - k_F)$, where $\Theta(z)$ is the Heaviside function and $k_F$ is such that $\epsilon(k_F) = 0$ and $k_F = \pi/2$, which corresponds to half filling. The system is critical and can be described with a conformal field theory with central charge $c = 1$ [81].

### B. Scaling limit

Using the Wigner function, we can express the noise-averaged energy density with respect to the ground state:

$$ \mathbf{e}_F(\tau) \equiv \overline{\langle \hat{h}_i \rangle_{\tau} - \langle \hat{h}_i \rangle_0} = \int \frac{dk}{2\pi} \epsilon(k) [n_{\tau}(k) - n_0(k)]. \quad (46) $$

One can derive (see Appendix G) an exact evolution equation for $n_{\tau}(k)$, which reads

$$ \partial_{\tau} n_{\tau}(k) = \tau_0 \int \frac{dk'}{2\pi} \tilde{F}(k') \epsilon(k + k'/2)^2 [n_{\tau}(k + k') - n_{\tau}(k)]. \quad (47) $$

By studying this equation (see Appendix G), we show that, around the Fermi points, $n_{\tau}(k)$ takes the scaling form

$$ n_{\tau}(k) \approx \mathbf{n}(\xi(k_F - k), \tau/\xi) + \mathbf{n}(\xi(k + k_F), \tau/\xi). \quad (48) $$

In turns, this leads to the result for fixed $t$ as $\xi \to +\infty$:

$$ \lim_{\xi \to \infty} \xi^2 \mathbf{e}_F(t\xi) = \tilde{\mathbf{e}}^+(t) + \tilde{\mathbf{e}}^-(t), \quad (49) $$

where we split the contribution of the energy from the two Fermi points expanding $k = \pm(k_F - p/\xi)$. In particular, introducing the Fermi velocity as $\epsilon'(\pm k_F) = \pm v$, we have for both chiral components

$$ \tilde{\mathbf{e}}^{\pm}(t) = -v \int_{-\infty}^{\infty} \frac{dp}{2\pi} p [\mathbf{n}(p, t) - \mathbf{n}(p, 0)] = \mathbf{e}(t)/2. \quad (50) $$

The last equality is proven in Appendix G, with $\mathbf{e}(t)$ given in Eq. (40) at $c = 1$. Hence, the CFT predicts, upon rescaling, the mean energy for the fermion system, thus confirming the exponential growth from a first-principles lattice calculation. This result suggests that, in the scaling limit of large $\xi$, the noisy dynamics in Eq. (44) is fully captured by the universal description provided by the CFT, upon rescaling space and time as $j = x\xi$, $\tau = t\xi$ and setting $\tau_0 = \xi$.

In order to validate this hypothesis, we compute numerically the two-point correlation matrix $C_{ij}(\tau) \equiv \langle a_i^\dagger a_j \rangle_{\tau}$. Since the model in Eq. (44) is noninteracting and the initial ground state is Gaussian, for each realization of the noise, all quantities can be expressed via the Wick theorem in terms of the coefficients $C_{ij}$. Nonetheless, we stress that, despite the Gaussianity of the quantum state, the distribution over different realizations of the noise of quantum expectation remains hard to access analytically. In Fig. 2 (top), we show the convergence for $\xi \to \infty$ of the noise-averaged energy density $\mathbf{e}_F(\tau)$ to its CFT prediction, consistently with Eqs. (49) and (50). The correlation matrix $C(\tau)$ can also be used to compute explicitly the Rényi entropies for any interval $I$ of size $\ell_F$. Indeed, setting $C_I$ as the $\ell_F \times \ell_F$ submatrix obtained restricting the indexes of $C(\tau)$ to $I$, we have for the von Neumann entropy of the interval $I: S_F(\tau) \equiv -\mathrm{Tr}[C_I \ln C_I + (1 - C_I) \ln(1 - C_I)]$. In Fig. 3 (left), we show the noise-averaged entanglement entropy production $S_F(\tau) \equiv S_F(\tau) - S_F(0)$ for the fermion system for intervals of various sizes $\ell_F$ on the lattice. Our prediction is that it should equal at large $\xi$ the CFT value $S_I^{(1)}$ (without any prefactor) with $\ell = \ell_F/\xi$.

We also predict that the CFT describes the distribution over the noise of these quantities. We show in Fig. 3
<figure><img src="image_4.png" /><figcaption>FIG. 2. Top: the average energy $\mathcal{E}_F(\tau)$ vs $t = \tau/\xi$, the scaled time, for different values of the noise correlation length $\xi$. Continuous lines are obtained from the numerical solution of the Wigner function equation (47), while the markers correspond to the exact dynamics of Eq. (44) for $L = 2048$. The dot-dashed line is the CFT result (40), which from Eq. (49) is predicted to hold for large $\xi$. Bottom: the median of the distribution of $\xi^2(\langle\hat{h}_i\rangle_\tau - \langle\hat{h}_i\rangle_0)$ vs the scaled time $\tau/\xi$. In the limit of $\xi \to \infty$, the median is expected to decrease toward the negative asymptotic value predicted by CFT (dot-dashed horizontal line). For finite $\xi$, the median starts to grow at large times, suggesting that heating may eventually dominate on the lattice.</figcaption></figure>

(middle) the one-point PDF for the local energy density $\xi^2\langle\hat{h}_i\rangle_\tau$ at the largest time $\tau$ available. As one sees, it compares reasonably well, with no free parameters, with the prediction from the CFT, i.e., the convolution $P_{\text{stat}} * P_{\text{stat}}$ where $P_{\text{stat}}$ is obtained in Eq. (4). This confirms that with the chosen $f(x)$, at this observation time, $2v_F\tau/\xi$ is large enough so that the two chiral components are only weakly correlated. As we see in Fig. 2 (top), the average energy grows with time, consistent with $P_{\text{stat}}$ having an infinite first moment. The median of the energy distribution, shown in Fig. 2 (bottom), is thus a better probe of the typical behavior. Remarkably, it is found to *decrease* with time, approaching at large $\xi$ a stationary value compatible with the CFT prediction $e_{\text{stat}}^{\text{median}} = -(c/2\pi)\tilde{\kappa}_0 < 0$; see Eq. (24). Finally, in Fig. 3 (right), as a representative of the finite-$\ell$ behavior, we compare the distribution of the entanglement entropy at different times for intervals of size $\ell_F = \ell\xi$ and $\ell = 1/2$, with the analytic prediction obtained from CFT.

### C. Deviation from the CFT prediction

For any finite $\xi$, we expect lattice effects to eventually break the CFT description. Both the average energy [Fig. 2 (top)] and the entanglement entropy [Fig. 3 (left)] indicate the emergence of a characteristic timescale $\tau^*(\xi)$, diverging with $\xi$, after which the CFT description becomes inaccurate. It is expected that the ultraviolet cutoff induces *heating* and stationarity does not hold anymore, as hinted by the rebound in the median observed at larger times [Fig. 2 (bottom)]. The solvability of the free-fermion model gives us the possibility to estimate this timescale. In principle, the full scaling function $n(p; t)$ in Eq. (48) should be universal, i.e., independent of the microscopic details of the discrete

<figure><img src="image_5.png" /><figcaption>FIG. 3. Left: noise-averaged von Neumann entanglement entropy production $\overline{S^{(1)}(\tau)}$ vs $t = \tau/\xi$ evolving in time under Eq. (44), for increasing values of $\xi$ and fixed ratio $\ell_F/\xi = \ell = 1/2$. The dotted line is obtained from the numerical solution of Eqs. (16) and using Eq. (34) for $n \to 1$, i.e., $\overline{S_F(\tau)} = c\bar{\kappa}/6$ with $c = 1$. Middle: distribution of the scaled energy density $\xi^2\langle\hat{h}_i(\tau)\rangle$ at $\xi = 64$. For the analytical prediction, we use that in the scaling limit $\xi^2\langle\hat{h}_j\rangle \to h(x, t)$ and at large time $h(x, t)$ is distributed as Eq. (5), which in the present case reduces to $h^{\text{in law}} = (3\Omega - 5)/(24\pi)$. In the inset, the right tail of the distribution is shown in log-log scale, showing the predicted $\propto h^{-3/2}$ tail. Right: distribution of the entanglement entropy $S_F(\tau)$ at $\xi = 64$ for an interval of size $\ell_F = 32$. For the analytical prediction, we use that in the stationary limit $\lim_{\tau\to\infty} S_F(\tau)^{\text{in law}} = c(\kappa^+ + \kappa^-)/12$, with $\kappa^\pm$ independently distributed according to $P_{\text{stat}}(\kappa)$. The stationary distribution $P_{\text{stat}}(\kappa)$ is obtained by numerically inverting the Fourier transform $Q_k(1)$ as a function of $k$, for $\ell = 1/2$, defined in Eqs. (E9) and (E39) in Appendix E. All simulations are performed on systems of total length $L = 2048$ and are repeated for $N_{\text{sample}} = 800$ samples.</figcaption></figure>
model, but depending on $f(x)$. Using that $n'(p; t) = \partial_p n(p; t)$ can be formally interpreted as a probability distribution (although not strictly positive), its time evolution can be characterized in terms of the moments

$$M_n(t) = \int_{-\infty}^{\infty} \frac{dp}{2\pi} p^n n'(p, t). \quad (51)$$

Since $n(p, t)$ ranges from 0 to 1 when $p$ increases, one has $M_0(t) = 1$. Also, integrating by part and using Eq. (50), $M_2(t) = 2\bar{\xi}^+(t)/v$. For higher $n$, one can show that these moments satisfy a hierarchy of differential equations which connects each moment with the previous ones having the same parity, i.e.,

$$\partial_t M_n(t) = \sum_{k=1}^{[n/2]} \alpha_k M_{n-2k}(t). \quad (52)$$

Because of the conservation of the density, one has $M_1(t) = 0$, and similarly all odd moments $M_{2k+1}(t)$ vanish. One can, thus, interpret $M_2(t)$ as the width of the distribution $n'(p; t)$. We use this result to estimate the timescale for the breakdown of the CFT description at least at the level of average quantities. Indeed, the Fermi points of $k \sim \pm k_F$ independently broaden with time, up to a time $\tau^*(\xi)$, where their width $\sim \xi\sqrt{M_2(t)}$ is comparable with their initial separation $\sim 2k_F$. From Eq. (50), $M_2(t) = \mathbf{e}(t)/v$, and Eq. (40) leads to

$$\tau^* \sim \frac{\xi}{v^2|f''(0)|} \ln \left( \frac{48\pi|f''(0)|k_F^2\xi^2}{f^{(4)}(0)} \right) \propto \xi \ln \xi. \quad (53)$$

Although we do not attempt a systematic verification of such a scaling, it is in qualitative agreement with our numerics, where the time at which the deviation from CFT occurs in the rescaled variable $t = \tau/\xi$ grows linearly in $\ln \xi$. Even thought the growth is only logarithmic, this shows that, consistently with the sizes of our numerics, a time window to observe the CFT phenomenology exists for systems involving a few thousands of lattice sites.

## IX. CONCLUSION

We have identified an out-of-equilibrium protocol which leads to a nontrivial stationary state [41] for a generic gapless one-dimensional system.

Several questions and directions remain open. First, it would be exciting to see the fingerprints of our predictions in a concrete experimental setup. As already mentioned, the study of the statistics of trajectories and unravelings has recently received a lot of attention in particular in the context of the measurement-induced phase transition. In contrast with these studies, our protocol is not cursed by postselection and is, thus, a promising platform for concrete observability. In order to facilitate experimental

observation, the robustness of the observed phenomenology has to be addressed, for instance, analyzing the role of thermal fluctuations. Our preliminary investigations suggest that a small initial temperature does not modify the power-law behavior of the energy density distribution. Second, it would be of great interest to obtain the full space-time statistics of the local energy, or of any other local operator, beyond the one-point distribution, especially since the latter exhibit heavy tails. It also remains a challenge to extend the present method to four-point (and higher) quantum correlation functions, thus providing a full characterization of the quantum state. For instance, the mutual information between two intervals requires a four-point correlator of the twist fields and reflects finer details of the specific CFT than the mere central charge. Although we focused on the infinite system, coupling between the chiral components becomes relevant at finite volume and can modify the behavior of the system at large times.

From a more concrete perspective, it would be interesting to test the theory at other quantum critical points beyond noninteracting systems, in particular, to observe the role played by interactions in controlling the deviations from the CFT predictions. We expect a connection between the recently derived large-scale description known as generalized hydrodynamics [82,83] to be helpful, in particular, in its zero-temperature extension toward quantum fluctuations [84].

Finally, the surprising existence of a stationary distribution in our setup raises the question about the fundamental ingredients to observe similar phenomenology in other quantum stochastic systems. It is possible that the stationary state that we identified within the continuous field theory description corresponds to a long-lived prethermal state that delays heating in the corresponding lattice system, similarly to what has been observed in the context of many-body quantum scars [85]. More generally, it remains an open question whether lattice effects or the presence of finite correlation length are compatible with the emergence of nontrivial steady states in the thermodynamic limit.

## ACKNOWLEDGMENTS

P. L. D. acknowledges support from the ANR Grant No. ANR-17-CE30-0027-01 RaMaTraF. D. B. acknowledges support from the project “ESQuisses,” ANR-20-CE47-0014-01. A. D. L. acknowledges support by the ANR JCJC Grant No. ANR-21-CE47-0003 (TamEnt).

## APPENDIX A: FOKKER-PLANCK EQUATION

In this section, we derive the Fokker-Planck equation (10) of the text, for the joint PDF of the backward stochastic trajectories $x_1 = X_t^+(y_1), \dots, x_n = X_t^+(y_n)$ associated to the Langevin equation (7). We thus consider
only a given chirality—here, we choose +, but the same Fokker-Planck equation holds for the chirality −, with $v \to -v$. We do not consider here the joint PDF of both chiralities. So here we denote simply $X_t^+ \to X_t$.

One can show that Eq. (7) translates into a stochastic equation for the variable $X_t(y_i)$ as a function of $t$ which takes the form [see Eqs. (58)and (59) in Supplemental Material of Ref. [56]]

$$dX_t(y_i) = \frac{v^2}{2} X_t''(y_i) f(0) dt - X_t'(y_i) v [dt + dW_t(y_i)], \quad (A1)$$

where $dX_t(y_i) = X_{t+dt}(y_i) - X_t(y_i)$. The $W_t(y_i)$ are mutually correlated Wiener processes in time $t$, which relates to the noise in Eqs. (1) and (7) via

$$W_t(y) = \int_0^t ds \eta(y, s), \quad \overline{dW_t(y)dW_t(y')} = dt f(y - y'). \quad (A2)$$

Here, $dW_t(y) = W_{t+dt}(y) - W_t(y)$, and $f(y)$ is the noise correlation function defined in the text.

Consider an arbitrary smooth function of $n$ variables $G(x_1, ..., x_n)$. In the case of these variables being the

backward stochastic trajectories $x_i = X_t(y_i)$ of Eq. (7), we define

$$g_t(y_1...y_n) = G[X_t(y_1)...X_t(y_n)]. \quad (A3)$$

By Ito calculus, the time variation $dg_t = g_{t+dt} - g_t$ of this observable is obtained by expanding up to second order:

$$\begin{aligned} dg_t = & \sum_{j=1}^{n} \partial_j G[X_t(y_1), ..., X_t(y_n)] dX_t(y_j) \\ & + \frac{1}{2} \sum_{j,m=1}^{n} \partial_j \partial_m G[X_t(y_1), ..., X_t(y_n)] dX_t(y_j) dX_t(y_m). \end{aligned} \quad (A4)$$

Here, we shorten the notation by setting $\partial_{x_j}...\partial_{x_m} G(x_1...x_n) = \partial_j...\partial_m G(x_1...x_n)$. Using Eqs. (A1) and (A2), we derive

$$dX_t(y_i)dX_t(y_j) = v^2 X_t'(y_i)X_t'(y_j)f(y_i - y_j)dt + O(dt^{3/2}). \quad (A5)$$

Averaging over the noise, we obtain

$$\begin{aligned} \frac{d\overline{g_t}}{dt} = & -v \sum_{j=1}^{n} \overline{\partial_j G[X_t(y_1), ..., X_t(y_n)] X_t'(y_j)} + \frac{v^2 f(0)}{2} \sum_{j=1}^{n} \overline{\partial_j G[X_t(y_1), ..., X_t(y_n)] X_t''(y_j)} \\ & + \frac{v^2}{2} \sum_{j,m=1}^{n} f(y_j - y_m) \overline{\partial_j \partial_m G[X_t(y_1), ..., X_t(y_n)] X_t'(y_j) X_t'(y_m)}. \end{aligned} \quad (A6)$$

From the chain rule for the derivation with respect to the variables {$y_i$}, it is easy to check that

$$\partial_{y_j} \partial_{y_m} g_t = X_t'(y_j) X_t'(y_m) \partial_j \partial_m G + \delta_{j,m} X_t''(y_j) \partial_j G, \quad (A7)$$

which finally leads to

$$\frac{d\overline{g_t}}{dt} = \left[ -v \sum_{j=1}^{n} \partial_{y_j} + \frac{v^2}{2} \sum_{j,m=1}^{n} f(y_j - y_m) \partial_{y_j} \partial_{y_m} \right] \overline{g_t}. \quad (A8)$$

It is useful to reexpress Eq. (A8) in terms of the Fokker-Planck Hamiltonian (and its Hermitian adjoint). In order to do so, we introduce the operators $q_j$ and $p_j$, defined by their action on any smooth function $\omega(\mathbf{y})$ as $q_j \cdot \omega(\mathbf{y}) = y_j \omega(\mathbf{y})$ and $p_j \cdot \omega(\mathbf{y}) = -i \partial_j \omega(\mathbf{y})$, respectively. For these conjugate variables, the canonical quantization holds $[q_i, p_j] = i \delta_{i,j}$. We can then define

$$\mathcal{H}_{\text{FP}} \equiv -iv \sum_{i=1}^{n} p_i + \frac{v^2}{2} \sum_{ij} p_i p_j f(q_i - q_j), \quad (A9)$$

so that we can rewrite

$$\begin{aligned} \frac{d\overline{g_t}}{dt} = & -\mathcal{H}_{\text{FP}}^{\dagger} \cdot \overline{g_t}, \\ \mathcal{H}_{\text{FP}}^{\dagger} = & iv \sum_{i=1}^{n} p_i + \frac{v^2}{2} \sum_{ij} f(q_i - q_j) p_i p_j. \end{aligned} \quad (A10)$$

To deduce the Fokker-Planck equation for $P(\mathbf{x}|\mathbf{y})$, defined in the text, we need one more step. From the definition of $P(\mathbf{x}|\mathbf{y})$, we have

$$\overline{g_t(\mathbf{y})} = \int d\mathbf{x}' G(\mathbf{x}') P_t(\mathbf{x}'|\mathbf{y}). \quad (A11)$$

We can choose $G(\mathbf{x}') = \delta_{\epsilon}(\mathbf{x}' - \mathbf{x})$, where $\delta_{\epsilon}(\mathbf{x})$ is a mollifier of the Dirac delta function. In the limit $\epsilon \to 0$, we recover $\overline{g_t(\mathbf{y})} \to P_t(\mathbf{x}|\mathbf{y})$ the JPDF for the initial points $\mathbf{x} = (x_1, ..., x_n)$ of the stochastic trajectories and finally deduce from Eq. (A10)
$$ \partial_t P_t(\mathbf{x}|\mathbf{y}) = -\mathcal{H}_{\text{FP}}^T[\mathbf{y}] \cdot P_t(\mathbf{x}|\mathbf{y}). \quad (A12) $$

This equation can be formally solved by starting from the initial condition at $t = 0$ $P_t(\mathbf{x}|\mathbf{y}) = \delta(\mathbf{x} - \mathbf{y})$. It is useful to employ the bra $\langle \mathbf{y} |$ and ket $| \mathbf{x} \rangle$ notation for the eigenstates of the position operators $\hat{q}_j$. Then, we can represent the probability distribution as

$$ P_t(\mathbf{x}|\mathbf{y}) = \langle \mathbf{y} | e^{-r\mathcal{H}_{\text{FP}}^T} | \mathbf{x} \rangle = \langle \mathbf{x} | e^{-r\mathcal{H}_{\text{FP}}} | \mathbf{y} \rangle, \quad (A13) $$

where the last equality follows from the Hermitian conjugation. Therefore, the following equation must also hold:

$$ \partial_t P_t(\mathbf{x}|\mathbf{y}) = -\mathcal{H}_{\text{FP}}[\mathbf{x}] \cdot P_t(\mathbf{x}|\mathbf{y}), \quad (A14) $$

where the action of the differential operator $\mathcal{H}_{\text{FP}}[\mathbf{x}]$ is now over the variables $\mathbf{x}$. More explicitly, using Eq. (A9), we arrive at

$$ \partial_t P_t(\mathbf{x}|\mathbf{y}) = v \left( \sum_{i=1}^{n} \partial_i + \frac{v}{2} \sum_{i,j=1}^{n} \partial_i \partial_j f(x_i - x_j) \right) P_t(\mathbf{x}|\mathbf{y}), \quad (A15) $$

which coincides with Eq. (10) given in the main text.

## APPENDIX B: JOINT EVOLUTION OF $r$ AND $\kappa$

In this section, we give two equivalent methods to derive the joint evolution of the variables $r$ and $\kappa$ defined in the text. The first one uses the Fokker-Planck equation, and the second one is a direct derivation using stochastic equations.

### 1. Derivation of the evolution equation for $P_t(r, x'_1, x'_2)$

We first derive the evolution equation for the JPDF $P_t(x_1, x_2, x'_1, x'_2)$ of the random variables $x_1, x_2, x'_1$, and $x'_2$ defined in the text. The starting point is the application of the Fokker-Planck equation (A15) in the case of four trajectories $(x_1, x_2, x_3, x_4)$. To fix the variables $y_i$, we remind that $y_1$ and $y_2$ are kept fixed while two additional variables are taken infinitesimally away from $y_1$ and $y_2$. More explicitly,

$$ \begin{aligned} x_1 &= X_t^{\pm}(y_1), & x_2 &= X_t^{\pm}(y_2), \\ x_3 &= X_t^{\pm}(y_1 + \epsilon), & x_4 &= X_t^{\pm}(y_2 + \epsilon). \end{aligned} \quad (B1) $$

We then perform the change of variables $x_1, x_2, x'_1 = (x_3 - x_1)/\epsilon, x'_2 = (x_4 - x_2)/\epsilon$ together with the limit $\epsilon \to 0$. By applying this change of variables in the Fokker-Planck equation above and taking the limit $\epsilon \to 0$, one obtains the following equation for $P_t = P_t(x_1, x_2, x'_1, x'_2)$:

$$ \begin{aligned} \partial_t P_t = & v \left( (\partial_{x_1} + \partial_{x_2}) + \frac{v}{2} f(0) (\partial_{x_1}^2 + \partial_{x_2}^2) \right. \\ & \left. \qquad + v \partial_{x_1} \partial_{x_2} f(x_1 - x_2) - \frac{v}{2} F''(0) [\partial_{x'_1}^2 (x'_1)^2 + \partial_{x'_2}^2 (x'_2)^2] \right. \\ & \left. \qquad + v \partial_{x_2} \partial_{x'_1} x'_1 f'(x_1 - x_2) - v \partial_{x_1} \partial_{x'_2} x'_2 f'(x_1 - x_2) \right. \\ & \left. \qquad - v \partial_{x'_1} \partial_{x'_2} x'_1 x'_2 f''(x_1 - x_2) \right) P_t. \end{aligned} \quad (B2) $$

Then, we exploit the invariance under translation by making the change of variables from $(x_1, x_2, x'_1, x'_2)$ to $[R = (x_1 + x_2/2), r = (x_1 - x_2)/\ell, x'_1, x'_2]$. Finally, integrating out the center of mass variable $R$, we obtain

$$ \begin{aligned} \partial_t P = & v^2 \left( \frac{1}{\ell^2} \partial_r^2 [f(0) - f(r\ell)] - \frac{f''(0)}{2} [\partial_{x'_1}^2 (x'_2)^2 + \partial_{x'_1}^2 (x'_2)^2] \right. \\ & \left. \qquad - f''(r\ell) \partial_{x'_1} \partial_{x'_2} x'_1 x'_2 - \frac{1}{\ell} \partial_r f(r\ell) (\partial_{x'_1} x'_1 + \partial_{x'_2} x'_2) \right) P. \end{aligned} \quad (B3) $$

From this result, we now obtain the Fokker-Planck equation for the joint distribution $P_t(\kappa, r)$ of the variable $\kappa$ and the variable $r$ defined as

$$ P_t(\kappa, r) = \int dx'_1 dx'_2 P(r, x'_1, x'_2) \delta \left[ \kappa - \ln \left( \frac{r^2}{x'_1 x'_2} \right) \right]. \quad (B4) $$

Interestingly, $P_t(\kappa, r)$ satisfies a closed evolution equation. To obtain it, we compute the time derivative of Eq. (B4) using Eq. (B3), and then we integrate by part and obtain

$$ \begin{aligned} \frac{\ell^2}{v^2} \partial_t P = & \left[ \partial_r^2 [f(0) - f(r\ell)] + 4 \frac{f(0) - f(r\ell) + r\ell f'(r\ell)}{r^2} \partial_\kappa^2 \right. \\ & \left. \qquad + \ell^2 f''(0) (\partial_\kappa - \partial_\kappa^2) - \ell^2 f''(r\ell) \partial_\kappa^2 \right. \\ & \left. \qquad + 2\ell \partial_r \partial_\kappa f'(r\ell) + 4\partial_r [f(0) - f(r\ell)] \frac{1}{r} \partial_\kappa \right. \\ & \left. \qquad + 2 \frac{f(0) - f(r\ell)}{r^2} \partial_\kappa \right] P. \end{aligned} \quad (B5) $$

Then, it is easy to see that the FP equation in Eq. (B5) is equivalent to the stochastic equations for $r$ and $\kappa$ (in Ito convention) given in Eq. (16) of the main text.

### 2. Direct derivation of the stochastic equations for $r$ and $\kappa$

Instead of working with the Fokker-Planck equation (4), we consider now an equivalent system of stochastic equations. It is equivalent in the sense that the Fokker-Planck equation associated with this system coincides with Eq. (4). It provides an alternative way to derive Eq. (16) in the main text. This stochastic system involves the set of trajectories $x_i(t)$, with $dx_i(t) = x_i(t + dt) - x_i(t)$, and reads
$$ \frac{dx_i(t)}{dW_i(t)} = \mp vdt + dW_i(t), \\ \left. \frac{dW_i(t)dW_j(t)}{dW_i(t)} \right|_t = f[x_i(t) - x_j(t)]dt, \quad (B6) $$

where the $W_i(t)$'s are mutually correlated Wiener processes in time and should not be confused with $W_t(y)$ introduced in Eq. (A2). Formally, the expectation value $\bar{O}|_t$ used here is conditioned to the position of the trajectories $x_i(t)$ at time $t$. The $\mp$ refers to $\pm$ chiralities, but the drift term plays no role in the following and so is omitted in the following. Here, the initial condition $x_i(t=0) = y_i$ is assumed. Note that this system of equation can be also seen as the time reversal of Eq. (7).

Let us set $\ell = 1$ and restore $\ell$ later. One performs the linear change of variable $x_3 = x_1 + \epsilon x'_1$, $x_4 = x_2 + \epsilon x'_2$, and $x_1 - x_2 = r$, which leads to (we keep the time dependence implicit for convenience)

$$ \begin{aligned} dr &= dW_1 - dW_2, & dx'_1 &= \frac{1}{\epsilon}(dW_3 - dW_1), \\ dx'_2 &= \frac{1}{\epsilon}(dW_4 - dW_2). \end{aligned} \quad (B7) $$

The correlations can be expressed in terms of $r$, $x'_1$, and $x'_2$ (the center of mass decouples), and performing the limit $\epsilon \to 0$ one finds with $j = 1, 2$

$$ \begin{aligned} \overline{dr^2} &= 2[f(0) - f(r)]dt, & \overline{drdx'_j} &= -f'(r)x'_jdt, \\ \overline{(dx'_j)^2} &= -f''(0)(x'_j)^2dt, & \overline{dx'_1dx'_2} &= -f''(r)x'_1x'_2dt. \end{aligned} \quad (B8) $$

Now we have for $\kappa = \log(r^2/x'_1x'_2)$ using Ito's rule

$$ d\kappa = 2\frac{dr}{r} - \frac{dx'_1}{x'_1} - \frac{dx'_2}{x'_2} + g(r)dt, \quad (B9) $$

where $g(r)dt$ is the Ito drift which can be computed using Eq. (B8):

$$ \begin{aligned} & -\frac{1}{r^2}\overline{dr^2} + \frac{1}{2}\overline{(dx'_1)^2} + \frac{1}{2}\overline{(dx'_2)^2} = g(r)dt, \\ & g(r) = \frac{2}{r^2}[f(0) - f(r)] - f''(0), \end{aligned} \quad (B10) $$

and depends only on $r$. We also need the correlations of the noise part

$$ \begin{aligned} \overline{d\kappa dr} &= 2\frac{dr^2}{r} - \frac{\overline{drdx'_1}}{x'_1} - \frac{\overline{drdx'_2}}{x'_2} = B(r)dt, \\ B(r) &= \frac{4}{r}[f(0) - f(r)] + 2f'(r), \end{aligned} \quad (B11) $$

$$ \begin{aligned} \overline{d\kappa^2} &= 4\frac{\overline{dr^2}}{r^2} + \frac{(\overline{dx'_1})^2}{(x'_1)^2} + \frac{(\overline{dx'_2})^2}{(x'_2)^2} - 4\frac{\overline{drdx'_1}}{rx'_1} \\ &\quad - 4\frac{\overline{drdx'_2}}{rx'_2} + 2\frac{\overline{dx'_1dx'_2}}{x'_1x'_2} \\ &= 2C(r)dt, \\ C(r) &= 4\frac{f(0) - f(r) + rf'(r)}{r^2} - f''(0) - f''(r), \end{aligned} \quad (B12) $$

to which we must add $\overline{dr^2} = 2A(r)$ with $A(r) = f(0) - f(r)$. Restoring $\ell$, one finds the equations in the text.

## APPENDIX C: NUMERICAL SOLUTION OF EQ. (16)

In this section, we provide some details about how to numerically solve the stochastic equations for $\kappa$ and $r$. The main difficulty arises because $r$ can become typically very small [see Eq. (22)]. It is, thus, useful to rewrite such a system of SDE in terms of another variable:

$$ \rho = \ln r, \quad r = e^{\rho}. \quad (C1) $$

Then, using Ito's lemma, the stochastic equation for $\rho$ takes the form

$$ d\rho = e^{-\rho}vdW_1 - e^{-2\rho}v^2A[r(\rho)]dt. \quad (C2) $$

To solve the equation, we then discretize time $dt \to \Delta t$ and define

$$ \begin{aligned} \Delta W_1 &= \sqrt{2A(r)}\beta_1, \\ \Delta W_2 &= \sqrt{\frac{B(r)^2}{2A(r)}\beta_1 + \sqrt{\left(2C(r) - \frac{B(r)^2}{2A(r)}\right)\beta_2}}, \end{aligned} \quad (C3) $$

where $\beta_1$ and $\beta_2$ are independently Gaussian random variables with zero average and $\Delta t$ variance. To rewrite the equation in a more numerically stable way, we extract from the quantities $A(r)$ and $B(r)$ their leading behavior at small $r$:

$$ \tilde{A}(r) = A(r)/r^2, \quad \tilde{B}(r) = B(r)/r^3, \quad (C4) $$

so that $\tilde{A}(r)$ and $\tilde{B}(r)$ are both finite in the limit $r \to 0$. So we finally have the discrete evolution equations

$$ \rho(t + \Delta t) = \rho(t) + v\sqrt{2\tilde{A}(r)\Delta t\beta_1} - v^2\tilde{A}(r)\Delta t, \quad (C5) $$
$$ \kappa(t + \Delta t) = \kappa(t) + v^2 g(r) \Delta t + vr^2 \sqrt{\frac{\tilde{B}(r)^2}{2\tilde{A}(r)}} \Delta t \beta_1 + \sqrt{\left(2C(r) - \frac{r^4 \tilde{B}(r)^2}{2\tilde{A}(r)}\right)} \Delta t \beta_2. \quad (C6) $$

## APPENDIX D: ANALYSIS OF THE SHORT-DISTANCE REGIME FOR $\kappa$

In the limit $\ell \ll 1$, one has Eq. (20) in the main text that we report for convenience:

$$ \begin{aligned} dr &= rvdB_1, \\ d\kappa &= -\frac{\ell^2}{6} r^2 v \left( -\frac{v}{2} f^{(4)}(0) dt + dB_2 \right), \end{aligned} \quad (D1) $$

$$ \begin{aligned} \overline{dB_1 dB_1} &= -f''(0)dt, & \overline{dB_2 dB_2} &= -f^{(6)}(0)dt, \\ \overline{dB_1 dB_2} &= -f^{(4)}(0)dt. \end{aligned} \quad (D2) $$

We now show from these equations that $\kappa$ satisfies a closed SDE which leads to the stationary measure given in the

main text. We first of all solve the equation for the variable $r$, which takes the form

$$ r(t) = \exp \left[ vB_1(t) + \frac{1}{2} v^2 f''(0)t \right]. \quad (D3) $$

Injecting this solution in the equation for $\kappa$ and integrating in time, we arrive at

$$ \kappa(t) = -\frac{\ell^2}{6} \int_0^t e^{2vB_1(s)+v^2f''(0)s} v \left( -\frac{v}{2}f^{(4)}(0)ds + dB_2(s) \right). \quad (D4) $$

This equation gives already a closed representation for $\kappa(t)$.

We can further simplify Eq. (D4) by making use of the following reparametrization: For each $t$, we define, for $i = 1, 2$, $\tilde{B}_i(s') = B_i(t) - B_i(t-s')$, $s' \in [0, t]$, which is, thus, an equivalent Brownian process which measures the deviation from the final point $B_i(t)$ (which is kept fixed) of the original one. One clearly has $\tilde{B}_i(0) = B_i(0) = 0$, $\tilde{B}_i(t) = B_i(t)$, and $d\tilde{B}_i(s') = dB_i(t-s')$. We want now to rewrite Eq. (D4) in terms of the processes $\tilde{B}_i$. A little bit of care is needed for the stochastic integral. Indeed, writing explicitly the Ito integral, we have

$$ \begin{aligned} \int_0^t e^{2vB_1(s)+v^2f''(0)s} dB_2(s) &= \lim_{n \to \infty} \sum_{i=0}^{n-1} e^{2vB_1(s_i)+v^2f''(0)s_i} [B_2(s_{i+1}) - B_2(s_i)] \\ &= e^{2v\tilde{B}_1(t)+v^2f''(0)t} \lim_{n \to \infty} \sum_{j=1}^n e^{-2v\tilde{B}_1(\tilde{s}_j)-v^2f''(0)\tilde{s}_j} [\tilde{B}_2(\tilde{s}_j) - \tilde{B}_2(\tilde{s}_{j-1})], \end{aligned} \quad (D5) $$

where the $s_i$'s are a partition of $n$ elements of $[0, t]$, with $s_0 = 0$ and $s_n = t$. We define $\tilde{s}_j = t - s_i$ with $j = n - i$, which is an equivalent partition. Clearly, the last expression does not converge to a stochastic integral in the Ito form. We, thus, rewrite in the last term as

$$ e^{-2v\tilde{B}_1(\tilde{s}_j)} = e^{-2v\tilde{B}_1(\tilde{s}_{j-1})} e^{-2v[\tilde{B}_1(\tilde{s}_j) - \tilde{B}_1(\tilde{s}_{j-1})]} \quad (D6) $$

and then expand the second exponential, using $[\tilde{B}_1(\tilde{s}_j) - \tilde{B}_1(\tilde{s}_{j-1})][\tilde{B}_2(\tilde{s}_j) - \tilde{B}_2(\tilde{s}_{j-1})] = -f^{(4)}(0)(\tilde{s}_j - \tilde{s}_{j-1})$. We then arrive at

$$ \int_0^t e^{2vB_1(s)+v^2f''(0)s} dB_2(s) = e^{2v\tilde{B}_1(t)+v^2f''(0)t} \left[ \int_0^t e^{-2v\tilde{B}_1(s)-v^2f''(0)s} dB_2(s) + 2vf^{(4)}(0) \int_0^t e^{-2v\tilde{B}_1(s)-v^2f''(0)s} ds \right]. \quad (D7) $$

Applying these transformations to Eq. (D4), we obtain

$$ \kappa(t) = -\underbrace{\frac{\ell^2}{6} e^{2v\tilde{B}_1(t)+v^2f''(0)t}}_{\kappa_1} \underbrace{\int_0^t e^{-2v\tilde{B}_1(s)-v^2f''(0)s} v \left( \frac{3v}{2}f^{(4)}(0)ds + dB_2(s) \right)}_{\kappa_2}. \quad (D8) $$
Using Ito's lemma, $d\kappa = \kappa_2 d\kappa_1 + \kappa_1 d\kappa_2 + d\kappa_1 d\kappa_2$, which, after collecting different contributions, leads to

$$d\kappa = v^2 \left( \frac{\ell^2}{12} f^{(4)}(0) - f''(0)\kappa \right) dt - \frac{\ell^2 v}{6} d\tilde{B}_2(t) + 2vd\tilde{B}_1(t)\kappa. \quad (D9)$$

The correlations of the $d\tilde{B}_j(t)$ are the same as the ones of the $dB_j(t)$ in Eq. (D2). However, $\kappa(t)$ defined by Eq. (D8) is a different process in $t$ than $\kappa(t)$ defined by Eq. (D4). At fixed $t$ the two random variables have the same law, but as $t$ is varied the trajectories are different (since the relation between $\tilde{B}_i$ and $B_i$ involves $t$ explicitly). As a consequence, the two stochastic equations (D9) and (D4) are inequivalent, although they lead to the same single-time distribution for $\kappa(t)$. One illustration of that is that, while the second process converges, i.e., $\kappa(t) \to \kappa_\infty$, where the distribution of $\kappa_\infty$ is given below, the first process is ergodic (with the same law). A simpler example which allows one to understand better this point (using more explicit notations) is worked out in the remark below.

We can recast Eq. (D9) as an equation with a single Brownian process $dB$ [we are using that $a_1 dB_1 + a_2 dB_2 = \sqrt{-a_1^2 f''(0) - a_2^2 f^{(6)}(0) - 2a_1 a_2 f^{(4)}(0)} d\tilde{B}$, where $d\tilde{B}$ is a new Wiener process with standard normalization, $d\tilde{B}^2 = dt$]:

$$d\kappa = v^2 \left( \frac{\ell^2}{12} f^{(4)}(0) - f''(0)\kappa \right) dt + \frac{vd\tilde{B}}{6} \sqrt{-[\ell^4 f^{(6)}(0) - 24\ell^2 f^{(4)}(0)\kappa + 144f''(0)\kappa^2]}. \quad (D10)$$

With the change of variable, $\kappa = \kappa_0(\omega/\omega_0 - 1)$, where we define as in the main text (24)

$$\kappa_0 = -\frac{\ell^2 f^{(4)}(0)}{12 f''(0)}, \quad \omega_0 = \frac{1}{\sqrt{\frac{f^{(6)}(0)f''(0)}{f^{(4)}(0)^2} - 1}}. \quad (D11)$$

Note that $\omega_0$ is real and positive, as it is guaranteed by the positivity of the Fourier transform $\hat{f}(k) > 0$ of $f(x)$. Indeed,

$$\begin{aligned} & f^{(6)}(0)f''(0) - f^{(4)}(0)^2 \\ &= \int k^2 \hat{f}(k) \int k^6 \hat{f}(k) - \left( \int k^4 \hat{f}(k) \right)^2 > 0, \end{aligned} \quad (D12)$$

which is a consequence of the Cauchy-Schwarz inequality.

The SDE for $\omega$ becomes Eq. (25), which is solved by

$$\omega(t) = e^{\sqrt{8\theta}B(t) - 2\theta t} \left( \omega_0 + \sqrt{8\theta} \int_0^t e^{-\sqrt{8\theta}B(s) + 2\theta s} d\gamma_s \right). \quad (D13)$$

One recognizes a Bougerol variable with drifted Brownian motion in the exponent [73]. It is useful to do the change of variable (25):

$$\omega = \sinh Y, \quad Y = Y(\omega) = \operatorname{argsinh}\omega. \quad (D14)$$

One has

$$Y'(\omega) = \frac{1}{\sqrt{1 + \omega^2}} = \frac{1}{\cosh Y}, \quad Y''(\omega) = -\frac{\omega}{(1 + \omega^2)^{3/2}} \quad (D15)$$

and from Ito's rule follows Eq. (26) in the main text.

**Remark.**—We present a method to transform the random variable Eq. (D4) expressed as an integral into the solution of the stochastic process (D10). As mentioned above, the two describe random variables having the same single-time distribution, but they differ as stochastic processes, i.e., in the way the realization at time $t$ and $t+dt$ are connected. In particular, Eq. (D4) almost surely has a fixed limit $\kappa(t \to \infty)$ for each realization, while the latter is ergodic [being equivalent to a Langevin equation (25)]. Nevertheless, the distribution of $\kappa(t \to \infty)$ over the noise realizations is the same in the two cases. To further clarify the aspect, we include a self-contained minimal example of this mechanism where we make the notation more explicit. Consider the process

$$Z_t = \int_0^t ds e^{B(s)-s/2}, \quad (D16)$$

where $B(s)$ is a standard Wiener process with $B(0) = 0$. Since $B(s)$ is almost surely subleading with respect to $s/2$, for each realization of $B(s)$ the limit $Z_\infty = \lim_{t \to +\infty} Z_t$ exists almost surely. The value of $Z_\infty$ changes from realization to realization and is, thus, a random variable. On the other end, as in the previous discussion, we define at fixed $t \tilde{B}_t(s) = B(t) - B(t-s)$, which is also a standard Wiener process. Here, $t$ denotes the time around which the Brownian has been reflected and is indicated for clarity as a subscript. One can now write

$$\begin{aligned} Z(t) &= \int_0^t ds' e^{B(t-s')-t/2+s'/2} \\ &= e^{-t/2+\tilde{B}_t(t)} \int_0^t ds' e^{-\tilde{B}_t(s')+s'/2} = \tilde{Z}_t(s) \Big|_{s=t}, \end{aligned} \quad (D17)$$

where we define
$$ \tilde{Z}_t(s) = e^{-s/2 + \tilde{B}_t(s)} \int_0^s ds' e^{-\tilde{B}_t(s') + s'/2}, \quad (\text{D18}) $$

This is a new process for $s \in [0, t]$, which obeys the following stochastic equation:

$$ d\tilde{Z}_t(s) := \tilde{Z}_t(s + ds) - \tilde{Z}_t(s) = ds + \tilde{Z}_t(s) d\tilde{B}_t(s), \quad (\text{D19}) $$

where $\tilde{Z}_t(0) = 0$. By running this equation until $s = t$, we can recover $Z(t)$ via Eq. (D17). We can solve this stochastic equation by setting $Q_t(s) = \ln Z_t(s)$. One has

$$ \begin{aligned} dQ_t(s) &= d\tilde{B}_t(s) + \left(e^{-Q_t(s)} - \frac{1}{2}\right)dt \\ &= -V'[Q_t(s)]ds + d\tilde{B}_t(s) \end{aligned} \quad (\text{D20}) $$

with $V(Q) = e^{-Q} + Q/2$. This last equation is of Langevin type, and it implies a stationary distribution at large $s = t \to \infty$ for $Q_\infty$ and $Z_\infty$:

$$ P(Q_\infty) = C e^{-2V(Q_\infty)} \Rightarrow P(Z_\infty) = \frac{C}{Z^2} e^{-2/Z}, \quad (\text{D21}) $$

which leads to the known inverse gamma distribution for $Z_\infty$. Note that we choose $t$ sufficiently large so that the stationary measure is reached when $s = t$.

# APPENDIX E: STATIONARY MEASURE FOR $\kappa$ FOR ANY FINITE $\ell$

## 1. Backward method

We now apply the backward method to the full stochastic equation for $r(t)$ and $\kappa(t)$ in Eq. (16). As explained in the main text, we define

$$ Q_k(r_0, t) \equiv e^{-ik\kappa(t)^{r_0}}, \quad (\text{E1}) $$

where the superscript $r_0 = r(t=0)$ indicates the initial condition for the variable $r$. In the end, we set $r_0 = 1$ as it is required in our case, but it is useful to keep it free. One has

$$ Q_k(r_0, t + dt) = \overline{e^{-ik[v^2 g(r_0)dt + v dW_2(0)]} Q[r_0 + v dW_1(0), t]^{r_0}}; \quad (\text{E2}) $$

expanding with Ito's lemma and averaging, we arrive at

$$ \partial_t Q_k = v^2 [A_\ell(r) \partial_r^2 - ikB_\ell(r) \partial_r - k^2 C_\ell(r) - ikg_\ell(r)] \times Q_k(r, t) \quad (\text{E3}) $$

with the boundary conditions

$$ \begin{aligned} Q_k(r_0, 0) &= 1, & Q_k(r_0 = 0, t) &= 1, \\ \lim_{r_0 \to +\infty} Q_k(r_0, t) &= e^{-2\theta(k^2+ik)t}, \end{aligned} \quad (\text{E4}) $$

and we recall that

$$ \begin{aligned} A_\ell(r) &= \frac{f(0) - f(\ell r)}{\ell^2}, & B_\ell(r) &= 2 \frac{f'(\ell r)}{\ell} + 4 \frac{f(0) - f(\ell r)}{\ell^2 r}, \\ C_\ell(r) &= 4 \frac{f(0) - f(\ell r) + \ell r f'(\ell r)}{\ell^2 r^2} - f''(0) - f''(\ell r), \\ g_\ell(r) &= -f''(0) - 2 \frac{f(0) - f(\ell r)}{\ell^2 r^2}. \end{aligned} \quad (\text{E5}) $$

The second condition in Eq. (E4) comes from the fact that $d\kappa = 0$ and $dr = 0$ for $r = 0$, since $g(0) = 0$ and $C(0) = 0$. The third condition is obtained using the fact that the dynamics of $\kappa$ is pure diffusion at large $r_0$.

Let us denote $Q_k(r)$ the stationary solution of Eq. (E3). It, thus, satisfies

$$ A_\ell(r) Q_k''(r) - ikB_\ell(r) Q_k'(r) - [k^2 C_\ell(r) + ikg_\ell(r)] Q_k(r) = 0 \quad (\text{E6}) $$

with the boundary conditions

$$ Q_k(r_0 = 0) = 1, \quad Q_k(r_0 \to +\infty) = 0. \quad (\text{E7}) $$

From this stationary solution, one obtains the stationary measure $P_{\text{stat}}(\kappa)$ for $\kappa$ by Fourier inversion:

$$ P_{\text{stat}}(\kappa) = \int_{-\infty}^{+\infty} \frac{dk}{2\pi} e^{ik\kappa} Q_k(r=1). \quad (\text{E8}) $$

## 2. Schrödinger equation for the stationary measure

We can further simplify Eqs. (E4) and (E6) by removing the first derivative term. This can be achieved by setting

$$ Q_k(r) = \phi_k(r) G_k(r). \quad (\text{E9}) $$

We choose $\phi_k(r)$ so that

$$ \begin{aligned} \frac{\phi_k'(r)}{\phi_k(r)} &= ik \frac{B(r)}{2A(r)}, \\ \phi_k(r) &= e^{ik \int_0^r dr' [B(r')/2A(r')]}, \\ &= \left( \frac{\ell^2 r^2 f''(0)}{2[f(\ell r) - f(0)]} \right)^{ik} = e^{ik\kappa_0(\ell r)}, \end{aligned} \quad (\text{E10}) $$

where $\kappa_0(x)$ is defined in Eq. (29) of the main text. Then, one finds that $G_k(r)$ satisfies the Schrödinger equation

$$ -G_k''(r) - k(k+i)V(r)G_k(r) = 0 \quad (\text{E11}) $$

with the potential
$$ 
\begin{aligned}
V(r) &= \frac{\ell^2\{f'(\ell r)^2 + [f(0) - f(\ell r)][f''(\ell r) + f''(0)]\}}{[f(0) - f(\ell r)]^2} \\
&= -\frac{d^2}{dr^2}\log[f(0) - f(\ell r)] + \frac{\ell^2 f''(0)}{f(0) - f(\ell r)}
\end{aligned}
\quad (E12) $$

on the positive half-space $r \ge 0$ with the boundary condition inherited from Eq. (E7):

$$ G_k(0) = 1, \quad \lim_{r \to \infty} G_k(r) = 0. \quad (E13) $$

### 3. Symmetry of the stationary distribution function

Before solving various cases, let us indicate a nice symmetry property. One notes that the dependence in $k$ in Eq. (E11) is only via the prefactor $\gamma \equiv k(k+i)$ of the potential: $G_k(r) \equiv G(\gamma; r)$. This implies that

$$ \int d\kappa P_{\text{stat}}(\kappa) e^{-ik\kappa} = \phi_k(1)G_k(1) = e^{ik\kappa_0}G(\gamma; 1), \quad (E14) $$

where $\kappa_0$ is explicitly

$$ \kappa_0 = \frac{1}{ik} \log \phi_k(1) = -\log \left( \frac{2[f(\ell) - f(0)]}{\ell^2 f''(0)} \right). \quad (E15) $$

Hence, from Eq. (E8),

$$ 
\begin{aligned}
P_{\text{stat}}(\kappa) &= \int \frac{dk}{2\pi} e^{ik(\kappa+\kappa_0)} \mathcal{G}[k(k+i); 1] \\
&= e^{(\kappa+\kappa_0)/2} \int \frac{du}{2\pi} e^{iu(\kappa+\kappa_0)} \mathcal{G}\left(u^2 + \frac{1}{4}; 1\right),
\end{aligned}
\quad (E16) $$

where in the last equality we change variable $k = u - i/2$. It implies, in particular, that

$$ 
\begin{aligned}
P_{\text{stat}}(\kappa) &= e^{(\kappa+\kappa_0)/2} \tilde{P}(|\kappa + \kappa_0|), \\
\tilde{P}(y) &= \int \frac{du}{2\pi} e^{iu y} \mathcal{G}\left(u^2 + \frac{1}{4}; 1\right)
\end{aligned}
\quad (E17) $$

or, equivalently, the symmetry (reminiscent of a Nishimori condition or a Galavotti-Cohen theorem)

$$ \frac{P_{\text{stat}}(-\kappa_0 + u)}{P_{\text{stat}}(-\kappa_0 - u)} = e^u. \quad (E18) $$

### 4. Proof of the right 3/2 tail

Here, we show that, for a smooth noise correlation $f(x)$, the right tail of the stationary distribution is always a power-law $P_{\text{stat}}(\kappa) \propto \kappa^{-3/2}$ for $\kappa \to +\infty$, independently of $\ell$. Thanks to Eq. (E1), it is enough to prove the following expansion at small $k$ for its Fourier transform:

$$ Q_k(r=1) = 1 + C\sqrt{k} + O(k), \quad (E19) $$

where $C$ is a constant (see below). To prove Eq. (E19), we proceed as follows. First of all, since $\phi_k(r)$ is analytic in $k$, using Eq. (E9), we can focus on the small $k$ expansion of $G_k(r)$. The small $k$ behavior of the solution of Eq. (E11) can be obtained by setting $x = r\sqrt{\gamma}$, with $\gamma = k(k+i)$. Then, setting $G_k(r) = g(r\sqrt{\gamma})$, we can rewrite Eq. (E11) in the limit $k \to 0$ as

$$ -g'(x) - V_{\infty}g(x) = 0 \Rightarrow g(x) = e^{-\sqrt{-V_{\infty}x}}, \quad (E20) $$

where we set $V_{\infty} = \lim_{r \to \infty} V(r)$ and enforce the boundary conditions (E13). This implies

$$ G_k(r) \sim g(r\sqrt{\gamma}) = e^{-\sqrt{-V_{\infty}\gamma}r}, \quad \forall r = O(\gamma^{-1/2}). \quad (E21) $$

This is still not enough, because we require the expansion (E19) for $r=1$. However, fixing $\delta > 0$ and $r$, we can write

$$ 
\begin{aligned}
|G'_{k}(r) - G'_{k}(\delta/\sqrt{\gamma})| &= \left| \int_{\delta/\sqrt{\gamma}}^{r} dr' G''_{k}(r') \right| \\
&\le \gamma \int_{\delta/\sqrt{\gamma}}^{r} dr' |V(r')G_{k}(r')| \\
&\le K\gamma(r - \delta/\sqrt{\gamma}),
\end{aligned}
\quad (E22) $$

where we set $K = \sup_r |V(r)G_k(r)|$, which is finite for a sufficiently smooth $f(x) \in C^6$ fast decaying at infinity. As a consequence, at small $\gamma$,

$$ 
\begin{aligned}
G'_{k}(r) &= G'_{k}(\delta/\sqrt{\gamma}) + O(\delta\sqrt{\gamma}) \\
&= -\sqrt{-V_{\infty}\gamma}[g(\delta) + O(\delta)] \\
&\xrightarrow{\delta \to 0} -\sqrt{-V_{\infty}\gamma} + O(\gamma).
\end{aligned}
\quad (E23) $$

Finally, integrating over $r$,

$$ G_k(1) = 1 + \int_0^1 dr G'_k(r) = 1 - \sqrt{-V_\infty \gamma} + O(\gamma), \quad (E24) $$

which proves Eq. (E19) with $C = \sqrt{-V_{\infty}}$.

### 5. Small $\ell$ limit

At small $\ell$, one finds that the potential is a harmonic oscillator:

$$ V(r) = \frac{\ell^4 r^2 [f^{(4)}(0)^2 - f^{(6)}(0)f^{\prime\prime}(0)]}{36 f^{\prime\prime}(0)^2} + O(\ell^6 r^4) \quad (E25) $$

and
$$ \kappa_0 = \tilde{\kappa}_0 \ell^2 + O(\ell^4), \qquad \tilde{\kappa}_0 = - \frac{f^{(4)}(0)}{12 f''(0)}. \qquad (E26) $$

We see that $\kappa_0$ is $O(\ell^2)$ as $\ell \to 0$, so the random variable $\kappa = O(\ell^2)$ in this limit. Therefore, we can obtain an $\ell$-independent limit by scaling $k = \tilde{k}/\ell^2$. In terms of this variable, the potential term in the Schrödinger equation has a finite limit:

$$ -k(k+i)V(r) \approx \frac{4\tilde{k}^2 r^2 \tilde{\kappa}_0^2}{\omega_0^2}, \qquad (E27) $$

where we define as in the text

$$ \omega_0 = \frac{1}{\sqrt{\frac{f''(0)f^{(6)}(0)}{f^{(4)}(0)^2} - 1}}. \qquad (E28) $$

A solution of Eq. (30) with the potential (E25) and the boundary conditions (E13) can be expressed in terms of the Bessel function as

$$ G_k(r) = \frac{2^{3/4}}{\Gamma(1/4)} \left( \frac{\tilde{\kappa}_0 |\tilde{k}|}{\omega_0} \right)^{1/4} \sqrt{r} K_{1/4} \left( \frac{|\tilde{k}| r^2 \tilde{\kappa}_0}{\omega_0} \right), \qquad (E29) $$

and the prefactor is fixed, imposing that $G_k(r=0) = 1$. This leads to

$$ Q_k(r=1) = \frac{2^{3/4}}{\Gamma(1/4)} \left( \frac{\tilde{\kappa}_0 |\tilde{k}|}{\omega_0} \right)^{1/4} e^{i\tilde{k}\tilde{\kappa}_0} K_{1/4} \left( \frac{|\tilde{k}|\tilde{\kappa}_0}{\omega_0} \right), \qquad (E30) $$

which allows one to determine $P_{\text{stat}}(\kappa)$ by Fourier inversion from Eq. (E8).

One can check that this coincides with the result in the text, identifying $\tilde{\kappa} = \tilde{\kappa}_0[(\omega/\omega_0) - 1]$. Equivalently, we obtain the scaling form

$$ P_{\text{stat}}(\kappa) \stackrel{\ell \ll 1}{\approx} \frac{1}{\ell^2} \tilde{P}\left(\frac{\kappa}{\ell^2}\right), \\ \tilde{P}(\tilde{\kappa}) \equiv \frac{C\omega_0}{\tilde{\kappa}_0} \left[ 1 + \omega_0^2 \left( \frac{\tilde{\kappa} + \tilde{\kappa}_0}{\tilde{\kappa}_0} \right)^2 \right]^{-3/4}, \qquad (E31) $$

and one can check that the Fourier transform of $\tilde{P}$

$$ \int \frac{d\tilde{k}}{2\pi} e^{i\tilde{k}\tilde{\kappa}} \tilde{P}(\tilde{\kappa}) = Q_k(r=1) \qquad (E32) $$

as given in Eq. (E30). This can be seen by using the identity

$$ \int dx e^{ikx} \frac{1}{(1+x^2)^{3/4}} = \sqrt{2\pi} \frac{(2|k|)^{1/4} K_{1/4}(|k|)}{\Gamma(\frac{3}{4})}. \qquad (E33) $$

## 6. Large $\ell$ limit

At large $\ell$, under the hypothesis that $f(x)$ and its derivatives decay at infinity, the potential term reaches a constant value

$$ -k(k+i)V(r) \approx -k(k+i)\ell^2 \frac{f''(0)}{f(0)} \sim \frac{2i\tilde{k}\theta}{f(0)}, \qquad (E34) $$

where we use once again the scaling $k = \tilde{k}/\ell^2$, which implies again $\kappa = O(\ell^2)$, but with $\ell \to \infty$ in this case. From this potential, we immediately derive the solution respecting the boundary conditions (E13) in the form

$$ G_k(r) = e^{-\sqrt{2i\tilde{k}\theta/f(0)}r}. \qquad (E35) $$

Note that at large $\ell$

$$ \kappa_0(\ell) \stackrel{\ell \to \infty}{\approx} \log\{\ell^2[-f''(0)]/2f(0)\} \qquad (E36) $$

so that $k\kappa_0 \stackrel{\ell \to \infty}{\to} 0$ and, therefore, $Q_k(r) = G_k(r)$. Inverting the Fourier transform (E8), we obtain once again the stationary distribution

$$ P_{\text{stat}}(\kappa) \stackrel{\ell \gg 1}{\approx} \frac{1}{\ell^2} \tilde{P}\left(\frac{\kappa}{\ell^2}\right), \\ \tilde{P}(\tilde{\kappa}) \equiv \sqrt{\frac{\theta}{2\pi f(0)}} \frac{e^{-\theta/2f(0)\tilde{\kappa}}}{\tilde{\kappa}^{3/2}} \Theta(\tilde{\kappa}). \qquad (E37) $$

Equivalently, denoting $\kappa = \theta\ell^2\chi/f(0)$, one finds that $\chi$ is distributed according to $\mathcal{L}(\chi)$ in Eq. (33) in the text, i.e., the stable one-sided Levy distribution of index $1/2$.

## 7. Solvable cases for $f(x)$

For some particular choice of the noise correlation function $f(x)$, the potential $V_k(r)$ takes a form which is explicitly integrable. In Table I, we list a few interesting cases. Here, we focus on the case

$$ f(x) = 1/\cosh(x), \qquad (E38) $$

which is analytic and fast decaying. Setting $\gamma = k(k+i)$, the solution $G_k(r)$ respecting the boundary conditions (E13) can be expressed in terms of hypergeometric function

<table><caption>TABLE I. A few examples of noise correlation functions $f(x)$ leading to Schrödinger equations with a solvable potential $V(r)$.</caption><thead><tr><th>f(x)</th><th>V(r)/ℓ<sup>2</sup></th><th>ŷ(k)</th><th>Smoothness</th></tr></thead><tbody><tr><td>e<sup>-|x|</sup></td><td><sup>1</sup>&frasl;<sub>(1-e<sup>-rℓ</sup>)<sup>2</sup></sub></td><td><sup>2</sup>&frasl;<sub>1+k<sup>2</sup></sub></td><td>C<sup>0</sup></td></tr><tr><td>2e<sup>-|x|</sup> - e<sup>-2|x|</sup></td><td>-<sup>2</sup>&frasl;<sub>1-e<sup>-rℓ</sup></sub></td><td><sup>12</sup>&frasl;<sub>k<sup>4</sup> + 5k<sup>2</sup> + 4</sub></td><td>C<sup>2</sup></td></tr><tr><td><sup>1</sup>&frasl;<sub>cosh x</sub></td><td>- tanh(ℓr)<sup>2</sup></td><td><sup>π</sup>&frasl;<sub>cosh(kπ/2)</sub></td><td>C<sup>∞</sup></td></tr><tr><td><sup>1</sup>&frasl;<sub>(cosh x)<sup>2</sup></sub></td><td>-2 tanh(ℓr)<sup>2</sup></td><td><sup>πk</sup>&frasl;<sub>sinh(kπ/2)</sub></td><td>C<sup>∞</sup></td></tr></tbody></table>
$$G_k(r) = e^{-\sqrt{\gamma}\ell r} [1 + \tanh(\ell r)] \sqrt{\gamma} \frac{2 F_1\left(\frac{\sqrt{\gamma}}{2} - \frac{1}{2}\sqrt{\gamma + \frac{1}{4}} + \frac{1}{4}, \frac{\sqrt{\gamma}}{2} + \frac{1}{2}\sqrt{\gamma + \frac{1}{4}} + \frac{1}{4}; \sqrt{\gamma} + 1; \frac{1}{\cosh(\ell r)^2}\right)}{2 F_1\left(\frac{\sqrt{\gamma}}{2} - \frac{1}{2}\sqrt{\gamma + \frac{1}{4}} + \frac{1}{4}, \frac{\sqrt{\gamma}}{2} + \frac{1}{2}\sqrt{\gamma + \frac{1}{4}} + \frac{1}{4}; \sqrt{\gamma} + 1; 1\right)}. \quad (E39)$$

Equivalently, this expression can be represented in terms of generalized Legendre functions

$$G_k(r) = \frac{P_{(1/2)(\sqrt{1+4\gamma}-1)}^{-\sqrt{\gamma}}[\tanh(\ell r)]}{P_{(1/2)(\sqrt{1+4\gamma}-1)}^{\sqrt{\gamma}}(0)}. \quad (E40)$$

As a first check, we verify that the solution Eq. (E39) reproduces the known solutions in the small and large $\ell$ limits.

### a. Asymptotic limits $\ell \to \infty$

At large $\ell$, we simply have

$$\lim_{\ell \to \infty} Q_{k'/\ell^2}(r) = \lim_{\ell \to \infty} G_{k'/\ell^2}(r) = e^{-r\sqrt{i\ell'}} \quad (E41)$$

in agreement with Eq. (E35) [$\theta = 1/2$ for Eq. (E38)]. The limit in Eq. (E41) can be easily obtained using Eq. (E39) using that

$$\lim_{\gamma \to 0^2} F_1 \left( \frac{\sqrt{\gamma}}{2} - \frac{1}{2} \sqrt{\gamma + \frac{1}{4}} + \frac{1}{4}, \frac{\sqrt{\gamma}}{2} + \frac{1}{2} \sqrt{\gamma + \frac{1}{4}} + \frac{1}{4}; \sqrt{\gamma} + 1; x \right) = {}_2F_1 \left( 0, \frac{1}{2}, 1; x \right) = 1, \quad (E42)$$

irrespectively of $x$.

### b. Small $\ell$ check

In this case, the limit is less trivial as $k = k'/\ell^2$ becomes large in the limit of small $\ell$, so that simultaneously the parameters of the hypergeometric are diverging, while its argument is going to 1. Thus, we first apply the transformation between hypergeometric functions

$${}_2F_1(a, b, c; z) = \frac{(1-z)^{-a-b+c}\Gamma(c)\Gamma(a+b-c) {}_2F_1(c-a, c-b; -a-b+c+1; 1-z)}{\Gamma(a)\Gamma(b)} + \frac{\Gamma(c)\Gamma(-a-b+c) {}_2F_1(a, b; a+b-c+1; 1-z)}{\Gamma(c-a)\Gamma(c-b)}. \quad (E43)$$

Then, we use that

$$\lim_{\ell \to 0} {}_2F_1 \left[ \frac{1}{4} \left( 2\sqrt{\gamma} - \sqrt{1+4\gamma} + 1 \right), \frac{1}{4} \left( 2\sqrt{\gamma} + \sqrt{1+4\gamma} + 1 \right); \frac{1}{2}; \tanh(\ell^2) \right] = \frac{e^{k'/2} (k')^{1/4} \Gamma(\frac{3}{4}) I_{-1/4}(k'/2)}{\sqrt{2}} \quad (E44)$$

to recover, after some manipulations, Eq. (E30).

## APPENDIX F: DISTRIBUTION OF THE STRESS-ENERGY TENSOR

As we see in the main text, the distribution of the stress-energy tensor can be deduced from the one of $\kappa$ in the limit of small $\ell$, simply by using Eq. (38). An alternative approach is to compute explicitly the time evolution of the stress-energy tensor directly from the expression of the Hamiltonian.

### 1. Direct derivation

Consider the explicit expression of the Hamiltonian (1) in terms of the stress-energy tensor

$$\hat{H} = v \int dx [1 + \eta(x, t)] [\hat{T}^{+}(x) + \hat{T}^{-}(x)]. \quad (F1)$$

To derive the time evolution of an operator, we introduce the infinitesimal generator of time evolution

$$d\hat{H} = v \int dx[dt + dW_t(x)][\hat{T}^+(x) + \hat{T}^-(x)] \quad (F2)$$

and define the time evolution $U_t$ operator up to time $t$ with the equation

$$U_{t+dt} = e^{-id\hat{H}} U_t. \quad (F3)$$

The evolution of an operator $\hat{O}(t) \equiv U_t^\dagger \hat{O} U_t$ takes the form

$$\begin{aligned} d\hat{O}(t) &= U_t^\dagger e^{id\hat{H}} \hat{O} e^{-id\hat{H}} U_t - \hat{O}(t) \\ &= U_t^\dagger \left( t[d\hat{H}, \hat{O}] - \frac{1}{2}[d\hat{H}, [d\hat{H}, \hat{O}]] + \dots \right) U_t. \end{aligned} \quad (F4)$$
Note that because of the Ito's convention we need to keep terms up to the double commutators.

Consider the particular case of $\hat{O} = \hat{T}^{\pm}(y)$, the stress-energy tensor at the point $y$. We have for the equal-time commutator

$$[\hat{T}^{+}(x), \hat{T}^{+}(y)] = -i \left( 2\delta'(x-y)\hat{T}^{+}(y) - \delta(x-y)\hat{T}^{+'}(y) - \frac{c}{24\pi}\delta'''(x-y) \right), \quad (F5)$$

$$[\hat{T}^{-}(x), \hat{T}^{-}(y)] = i \left( 2\delta'(x-y)\hat{T}^{-}(y) - \delta(x-y)\hat{T}^{-'}(y) - \frac{c}{24\pi}\delta'''(x-y) \right). \quad (F6)$$

From this, we deduce

$$[d\hat{H}, \hat{T}^{\pm}(y)] = \pm 2i \hat{T}^{\pm}(y) v dW_{t}^{\prime}(y) \pm i [dt + dW_{t}(y)] v \hat{T}^{\pm'}(y) \mp i \frac{c}{24\pi} v dW_{t}^{\prime\prime}(y), \quad (F7)$$

where $dW_t'(x) = \partial_x dW_t(x)$ (the noise is smooth in space) and similarly for the higher derivatives. This implies the evolution equation for the operators $\hat{T}^{\pm}(y, t)$, which reads

$$d\hat{T}^{\pm}(y, t) = \mp \left( 2\hat{T}^{\pm}(y, t) v dW_{t}^{\prime}(y) + [dt + dW_{t}(y)] v \hat{T}^{\pm'}(y, t) - \frac{c}{24\pi} v dW_{t}^{\prime\prime}(y) \right) \\ + \frac{cv^{2}f^{(4)}(0)}{48\pi} dt - f''(0) v^{2} \hat{T}^{\pm}(y, t) dt + \frac{1}{2} f(0) v^{2} \partial_{y}^{2} \hat{T}^{\pm}(y, t) dt. \quad (F8)$$

Although the two chiral components do not couple at the CFT level, their evolutions are not statistically independent, since they feel the same noise. There are several quantities that one can study from there. One is the noise average $\overline{\hat{T}^{\pm}(y, t)}$, which is still a quantum operator. Its evolution is obtained by taking the noise average of Eq. (F8) and reads

$$\partial_t \overline{\hat{T}^{\pm}(y, t)} = \frac{cv^2 f^{(4)}(0)}{48\pi} - f''(0) v^2 \overline{\hat{T}^{\pm}(y, t)} + \frac{1}{2} f(0) v^2 \partial_y^2 \overline{\hat{T}^{\pm}(y, t)} \mp v \partial_y \overline{\hat{T}(y, t)}. \quad (F9)$$

Another observable is the quantum expectation $\langle \hat{T}^{\pm}(y, t) \rangle = \langle \Psi_0 | \hat{T}^{\pm}(y, t) | \Psi \rangle$ on any translational invariant state $|\Psi\rangle$ in a given noise realization. It satisfies a stochastic differential equations obtained by taking the quantum expectation of Eq. (F8) [which leads to a similar equation as Eq. (F8) but now for a scalar $\langle \hat{T}^{\pm}(y, t) \rangle$]. Note that it describes the coupled stochastic evolution of the two fields $\langle \hat{T}^{\pm}(y, t) \rangle$, a complicated problem. Here, we focus only on (i) noise moments, which have a solvable dynamics (ii) the one-point PDF of $\langle \hat{T}^{\pm}(y, t) \rangle$. This one-point distribution being independent of $y$ at $t=0$, it remains independent of $y$ for all times. In addition, it does not depend on the chirality; thus, we omit the $\pm$ superscript and denote $T(t) = \langle \hat{T}^{\pm}(y, t) \rangle$ the corresponding random variable. Then, one can show that its PDF can be obtained from the stochastic equation (F10):

$$dT = 2T v dB_1(t) - \frac{c}{24\pi} v dB_2(t) + \frac{cv^2 f^{(4)}(0)}{48\pi} dt \\ - f''(0) v^2 T dt \quad (F10)$$

with $dB_1(t) = \mp dW'(y, t)$ and $dB_2(t) = \mp dW'''(y, t)$. This can be seen intuitively: Indeed, we expect that the spatial derivative terms in Eq. (F8) are irrelevant, since the one-point PDF of $T = \langle \hat{T}^{\pm}(y, t) \rangle$ does not depend on $y$. More formally, one can prove that Eqs. (F10) and (F8) lead to the same evolution equation for the noise average $Z(\bar{T})$ for any smooth function $Z$. Although Eq. (F10) is obtained here by a completely different method, one can check that Eq. (F10) is equivalent to Eq. (D9) with the correspondence

$$T = \lim_{\ell \to 0} \frac{c}{4\pi \ell^2} \kappa, \quad (F11)$$

noting also that the noise satisfies Eq. (D2).

For the translationally invariant initial state chosen here, it is easy to obtain from Eq. (F10) the recursive equation for the moments (over the noise) of the quantum expectation of $T$ as

$$\partial_t \overline{T^n} = - \frac{v^2 c^2 f^{(6)}(0) n(n-1)}{1152 \pi^2} \overline{T^{n-2}} + \frac{c v^2 f^{(4)}(0) (4n-3) n}{48\pi} \overline{T^{n-1}} - n(2n-1) v^2 f''(0) \overline{T^n}. \quad (F12)$$
These can be solved with $T(t=0) = 0$, which corresponds to $|\Psi\rangle = |\Psi_0\rangle$ being the ground state. This leads to the first two moments:

$$ \bar{T} = \frac{c f^4(0)}{48 \pi f''(0)} (1 - e^{-v^2 f''(0)t}), \quad (\text{F13}) $$

$$ \bar{T}^2 = \frac{c^2}{3456\pi^2} \left[ \frac{1}{2} \left( \frac{f^{(4)}(0)}{f''(0)} \right)^2 (e^{-6v^2 f''(0)t} - 6e^{-v^2 f''(0)t} + 5) \right. \\ \left. + \frac{f^{(6)}(0)}{f''(0)} (e^{-6v^2 f''(0)t} - 1) \right]. \quad (\text{F14}) $$

# APPENDIX G: FREE FERMION DYNAMICS

## 1. Evolution equation for the Wigner function

Here, we derive the evolution equation Eq. (47) for the Wigner function. We consider a model of spinless non-interacting fermions in one dimension. Let us first consider the Hamiltonian in the absence of noise:

$$ \hat{H}_0 = \sum_i \hat{h}_i, \quad \hat{h}_i = -J(\hat{a}_{i+1}^\dagger \hat{a}_i + \hat{a}_i^\dagger \hat{a}_{i+1} - \mu \hat{a}_i^\dagger \hat{a}_i), \quad (\text{G1}) $$

where we introduce a chemical potential $\mu$ (chosen to be zero in the text).

In general, it is useful to represent densities which are conserved under the $H_0$ evolution with support around the site j in the following form:

$$ \hat{z}_j = \sum_{j'} \zeta_{j'} \hat{a}_{j+j'}^{\dagger} \hat{a}_j = \int_{-\pi}^{\pi} \frac{dk}{2\pi} \sum_{j'} e^{-ikj'} \zeta(k) \hat{a}_{j+j'}^{\dagger} \hat{a}_j, \\ \left[ \sum_j \hat{z}_j, \hat{H}_0 \right] = 0, \quad (\text{G2}) $$

where $\zeta(k) = \sum_j \zeta_j e^{ikj}$ is the Fourier transform of $\zeta_j$, and consider the coupling

$$ \hat{H} = \hat{H}_0 + \sum_j \eta_j(\tau)(\hat{z}_j + \hat{z}_j^\dagger), \quad (\text{G3}) $$

where the correlation of the noise $\eta_j(\tau)$ is defined in the main text. We keep $\zeta(k)$ arbitrary and specify its value only at the end.

One can introduce the Wigner function which completely characterizes all correlation functions in any Gaussian state. Here, we focus on the noise average Wigner function, which is defined as

$$ n_{\tau}(k) \equiv \sum_{j'} e^{ikj'} \overline{\langle \hat{a}_{j+j'}^{\dagger} \hat{a}_{j} \rangle}_{\tau} = \frac{1}{L} \sum_{jj'} e^{ikj'} \text{Tr}[\hat{a}_{j+j'}^{\dagger} \hat{a}_{j} \overline{\hat{Q}_{\tau}}], \quad (\text{G4}) $$

where $\langle \cdots \rangle_{\tau} = \langle \Psi_0(\tau) | \cdots | \Psi_0(\tau) \rangle$ denotes the quantum average at time $\tau$ and we introduce the density matrix

$\overline{Q}_{\tau} = |\Psi_0(\tau)\rangle\langle\Psi_0(\tau)|$. Although we are interested in the ground state of $\hat{H}_0$, these considerations apply to any translational invariant Gaussian initial state $|\Psi_0\rangle$. In this case, because of the noise average, Eq. (G4) is independent of the position j.

We now consider the quantum evolution of $\overline{Q}_{\tau}$. It is easy to verify that, after the noise average, one obtains the Lindblad form

$$ \frac{d\overline{Q}_{\tau}}{d\tau} = -t[H_0, \overline{Q}_{\tau}] - \frac{\tau_0}{2} \sum_{j,j'} F(j-j') [\hat{z}_j + \hat{z}_{j'}^{\dagger}, [\hat{z}_{j'} + \hat{z}_{j}^{\dagger}, \overline{Q}_{\tau}]]. \quad (\text{G5}) $$

We can use Eqs. (G4) and (G5) to obtain an evolution equation for $n_{\tau}(k)$. Note that the first term in Eq. (G5) does not contribute because of translational invariance and Eq. (G2). Using Eq. (G2), we obtain explicitly

$$ \partial_{\tau} n_{\tau}(k) = -\frac{\tau_0}{2} \sum_{n,n',m,j,j'} F(n-n') \times e^{ikm} \left( \zeta_j \zeta_{j'} \overline{\langle [a_{n+j}^{\dagger} \hat{a}_n, [a_{n'+j}^{\dagger} \hat{a}_{n'}, a_m^{\dagger} \hat{a}_0]] \rangle}_{\tau} \right. \\ \left. + \zeta_j^* \zeta_{j'}^* \overline{\langle [a_n^{\dagger} \hat{a}_{n+j}, [a_{n'}^{\dagger} \hat{a}_{n'+j}, a_m^{\dagger} \hat{a}_0]] \rangle}_{\tau} \right. \\ \left. + \zeta_j^* \zeta_{j'} \overline{\langle [a_n^{\dagger} \hat{a}_{n+j}, [a_{n'+j}^{\dagger} \hat{a}_{n'}, a_m^{\dagger} \hat{a}_0]] \rangle}_{\tau} \right. \\ \left. + \zeta_j \zeta_{j'}^* \overline{\langle [a_{n+j}^{\dagger} \hat{a}_n, [a_{n'}^{\dagger} \hat{a}_{n'+j}, a_m^{\dagger} \hat{a}_0]] \rangle}_{\tau} \right). \quad (\text{G6}) $$

We can expand all the commutators applying twice

$$ [\hat{a}_{j_1}^{\dagger} \hat{a}_{j_2}, \hat{a}_{j_3}^{\dagger} \hat{a}_{j_4}] = \delta_{j_2 j_3} \hat{a}_{j_1}^{\dagger} \hat{a}_{j_4} - \delta_{j_1 j_4} \hat{a}_{j_3}^{\dagger} \hat{a}_{j_2} \quad (\text{G7}) $$

from the anticommutation rules $\{\hat{a}_i^\dagger, \hat{a}_j^\dagger\} = 0$ and $\{\hat{a}_i^\dagger, \hat{a}_j\} = \delta_{ij}$, and one gets

$$ \begin{aligned} [\hat{a}_{j_1}^{\dagger} \hat{a}_{j_2}, [a_{j_3}^{\dagger} \hat{a}_{j_4}, a_{j_5}^{\dagger} \hat{a}_{j_6}]] = & \delta_{j_2 j_3} \delta_{j_4 j_5} a_{j_1}^{\dagger} \hat{a}_{j_6} - \delta_{j_1 j_6} \delta_{j_4 j_5} a_{j_3}^{\dagger} \hat{a}_{j_2} \\ & - \delta_{j_3 j_6} \delta_{j_2 j_5} a_{j_1}^{\dagger} \hat{a}_{j_4} + \delta_{j_3 j_6} \delta_{j_1 j_4} a_{j_5}^{\dagger} \hat{a}_{j_2}. \end{aligned} \quad (\text{G8}) $$

As a result, a single correlator survives. This is expected, since the full Hamiltonian (G3) is quadratic, so the correlation matrix $\langle \hat{a}_j^\dagger \hat{a}_{j'} \rangle$ must satisfy a closed and linear evolution equation. Next, one obtains, using the parity of the noise correlation $F(x) = F(-x)$ and translational invariance,

$$ \partial_{\tau} n_{\tau}(k) = -\frac{\tau_0}{2} \sum_{m,j,j'} e^{ikm} \{2(\zeta_j \zeta_{j'} + \zeta_{-j}^* \zeta_{-j'}^*)[F(j) - F(m+j)] \\ + (\zeta_j \zeta_{-j}^* + \zeta_{-j} \zeta_{j'}^*)[F(0) + F(j+j')] \\ -F(m) - F(m+j+j')\} \overline{\langle \hat{a}_{m+j+j'}^{\dagger} \hat{a}_{0} \rangle}_{\tau}. \quad (\text{G9}) $$
Using that $\zeta(k)^* = \sum_j e^{ikj} \zeta_{-j}^*$, we finally obtain

$$ \partial_{\tau} n_{\tau}(k) = \tau_0 \int_q \tilde{F}(q) |\zeta(-k) + \zeta(-q-k)^*|^2 [n(k+q) - n(k)], \qquad (G10) $$

where

$$ \tilde{F}(q) = \sum_j e^{iqj} F(j). \qquad (G11) $$

Consider first half filling, i.e., $\mu = 0$, as in the main text. To recover the coupling with the noise defined there, we must choose $\zeta(k) = -Je^{ik}$. With this choice,

$$ |\zeta(-k) + \zeta(-q-k)^*|^2 = \epsilon(k+q/2)^2, \qquad (G12) $$

and one recovers Eq. (47) given in the main text.

For finite $-2 < \mu < 2$, the Hamiltonian in Eq. (G1) remains critical and described by a $c = 1$ CFT. In order to preserve parity and time-reversal symmetry, a simple choice is $\zeta(k) = \epsilon(k)/2$ with $\epsilon(k) = -J(2 \cos k - \mu)$, which corresponds to the coupling to the noise of the form

$$ \hat{H} = \hat{H}_0 + \frac{1}{2} \sum_j \eta_j(\tau) (\hat{h}_j + \hat{h}_{j-1}). \qquad (G13) $$

In this case, we obtain after replacing $\zeta(k) \to \epsilon(k)/2$ in Eq. (G10)

$$ \partial_{\tau} n_{\tau}(k) = \frac{\tau_0}{4} \int_q \tilde{F}(q) [\epsilon(k) + \epsilon(q+k)]^2 [n(k+q) - n(k)]. \qquad (G14) $$

## 2. Scaling limit

Here, we study the evolution equation for the Wigner function (G14). In the scaling limit of large $\xi$ described in the main text, we look for a solution which has the following scaling form around the two Fermi points:

$$ n_{\tau}(k) \approx \mathbf{n}[\xi(k_F - k), \tau/\xi] + \mathbf{n}[\xi(k + k_F), \tau/\xi]. \qquad (G15) $$

The initial condition which corresponds to the ground state reads

$$ \begin{aligned} n_{\tau=0}(k) &= \Theta(k + k_F) - \Theta(k - k_F) \\ &= -1 + \Theta(k + k_F) + \Theta(k_F - k), \end{aligned} \qquad (G16) $$

which gives for the initial scaling function $\mathbf{n}(p, t=0) = -1/2 + \Theta(p)$. We can now derive directly an evolution equation for the scaling function $\mathbf{n}(p, t)$. In order to do this, we replace in Eq. (G14) $\hat{\xi}(k)$ with the energy dispersion relation $\epsilon(k)$, which gives Eq. (47) in the text. Next, we inject the scaling form Eq. (G15) around each Fermi point (choosing $+k_F$) and replace $k = k_F - p/\xi$, $k' = -p'/\xi$. From Eqs. (45) and (G11), we obtain that at large $\xi$, $\tilde{F}(k) \approx \xi \hat{f}(\xi k)$, which leads to (using $\tau_0 = \xi$)

$$ \begin{aligned} \partial_{\tau} \mathbf{n}(p; t) &\approx \frac{\xi^2}{4} \int_{-\pi\xi}^{\pi\xi} \frac{dp'}{2\pi} \tilde{f}(p') \left[ \epsilon\left(k_F - \frac{p}{\xi}\right) + \epsilon\left(k_F - \frac{p+p'}{\xi}\right) \right]^2 [\mathbf{n}(p+p'; t) - \mathbf{n}(p; t)] \\ &\xrightarrow{\xi \to \infty} v^2 \int_{-\infty}^{\infty} \frac{dp'}{2\pi} \tilde{f}(p') \left(p + \frac{p'}{2}\right)^2 [\mathbf{n}(p+p'; t) - \mathbf{n}(p; t)] \\ &= v^2 \int_{-\infty}^{\infty} \frac{dq}{2\pi} \tilde{f}(p-q) \left(\frac{q+p}{2}\right)^2 [\mathbf{n}(q; t) - \mathbf{n}(p; t)]. \end{aligned} \qquad (G17) $$

We can now compute the mean energy density with respect to the ground state $\mathbf{e}_F(\tau)$, which is expanded as

$$ \begin{aligned} \mathbf{e}_F(\tau = t\xi) &= \int_{-\pi}^{\pi} \frac{dk}{2\pi} \epsilon(k) [n_{\tau}(k) - n_0(k)] \\ &\approx \frac{1}{\xi^2} [\tilde{\mathbf{e}}^+(t) + \tilde{\mathbf{e}}^-(t)], \end{aligned} \qquad (G18) $$

$$ \begin{aligned} \tilde{\mathbf{e}}^+(t) &= \lim_{\xi \to \infty} \xi \int_{-\infty}^{\infty} \frac{dp}{2\pi} \epsilon\left(k_F - \frac{p}{\xi}\right) [\mathbf{n}(p, t) - \mathbf{n}(p, 0)] \\ &= -v \int_{-\infty}^{\infty} \frac{dp}{2\pi} p[\mathbf{n}(p, t) - \mathbf{n}(p, 0)] \end{aligned} \qquad (G19) $$

and similarly for $\tilde{\mathbf{e}}_F(t)$ with $k \sim -k_F$.

We now prove the validity of Eq. (49) in the main text. We show that $\tilde{\mathbf{e}}^+(t)$ satisfies a first-order differential equation. Indeed, differentiating Eq. (G18) with respect to the variable $t$ and using Eq. (G17), we obtain

$$ \begin{aligned} \partial_{\tau} \tilde{\mathbf{e}}^{+}(t) &= -v^3 \int_{-\infty}^{\infty} \frac{dp}{2\pi} \frac{dq}{2\pi} p \tilde{f}(p-q) \left(\frac{q+p}{2}\right)^2 \\ &\quad \times [\mathbf{n}(q; t) - \mathbf{n}(p; t)] \\ &= -\frac{v^3}{2} \int_{-\infty}^{\infty} \frac{dq}{2\pi} \frac{dp}{2\pi} \tilde{f}(p-q)(p-q) \left(\frac{q+p}{2}\right)^2 \\ &\quad \times [\mathbf{n}(q; t) - \mathbf{n}(p; t)], \end{aligned} \qquad (G20) $$

where in the last equality we symmetrize the integrand with respect to the exchange $p \leftrightarrow q$. The integral in Eq. (G20) is
finite, but, to proceed further, we need to split the integral in two terms involving $n(q; t)$ and $n(p; t)$, which are both divergent. To avoid this issue, we change variables $p = q + u$ and integrate by parts with respect to $q$ using that $n(q; t) - n(q + u; t)|_{q=-\infty} = 0$ to arrive at

$$ \partial_t \tilde{\mathbf{e}}^+(t) = v^3 \int \frac{dq \, du}{2\pi} \frac{\tilde{f}(u)u}{2\pi} \frac{(q+u/2)^3}{6} \times [\mathbf{n}'(q; t) - \mathbf{n}'(q+u; t)], \quad (G21) $$

where $\mathbf{n}'(q; t) = \partial_q \mathbf{n}(q, t)$. We can now split the integral into two finite terms and replace $q \to q - u$ in the second integral. This leads to

$$ \begin{aligned} \partial_t \tilde{\mathbf{e}}^+(t) &= v^3 \int \frac{dq \, du}{2\pi} \frac{\tilde{f}(u)u}{2\pi} \frac{1}{6} [(q+u/2)^3 - (q-u/2)^3] \mathbf{n}'(q; t) \\ &= v^3 \int \frac{du}{2\pi} \tilde{f}(u) \frac{u^4}{48\pi} \int dq \mathbf{n}'(q; t) \\ &\quad + \frac{v^3}{2} \int \frac{dq}{2\pi} q^2 \mathbf{n}'(q; t) \int \frac{du}{2\pi} \tilde{f}(u) u^2 \\ &= \frac{v^3 f^{(4)}(0)}{48\pi} - v^2 f''(0) \tilde{\mathbf{e}}^+(\tau), \end{aligned} \quad (G22) $$

where we use that

$$ \begin{aligned} \int dq \mathbf{n}'(q; t) &= 1, \\ \int \frac{dq}{2\pi} q^2 \mathbf{n}'(q; t) &= \int \frac{dq}{2\pi} q^2 [\mathbf{n}'(q; t) - \mathbf{n}'(q; 0)] \\ &= \frac{2}{v} \tilde{\mathbf{e}}^+(t), \end{aligned} \quad (G23) \\ \int \frac{du}{2\pi} \tilde{f}(u) u^2 &= -f''(0), \\ \int \frac{du}{2\pi} \tilde{f}(u) u^4 &= f^{(4)}(0). \quad (G24) $$

Note that in Eq. (G23) we use the initial condition (G16) which shows that $\mathbf{n}'(q, 0) = \delta(q)$. We can thus solve

Eq. (G22) with the initial condition $\tilde{\mathbf{e}}^+(t=0) = 0$ and obtain

$$ \tilde{\mathbf{e}}^+(t) = \frac{f^{(4)}(0)v(1 - e^{-v^2 f''(0)t})}{48\pi f''(0)}. \quad (G25) $$

Summing also the equivalent contribution from $\tilde{\mathbf{e}}^-(t)$ and using the definition of $\theta = -v^2 f''(0)/2$, we recover Eq. (40) by setting the central charge $c = 1$, which is expected from universality. Note that the present calculation is from first principles.

## APPENDIX H: ENERGY DISTRIBUTION AT FINITE TEMPERATURE

Let us assume now that the initial state is not at zero temperature. In this case, the operator identity Eq. (35) still holds. As explained in the main text, the random variable $T_\beta = \langle \hat{T}(y, t) \rangle_\beta$ can be expressed in terms of the variables $r$ and $\kappa$ in the small $\ell$ limit via Eq. (42). Since we are interested in the limit $\ell \to 0$, the stochastic equation describing $r$ and $\kappa$ is simply given by Eq. (20). Thus, we have

$$ dT_{\beta} = \frac{\pi c}{12\beta^2 v^2} [2r^2 v dB_1 - f''(0) r^2 v^2 dt] \\ - \frac{cvr^2}{24\pi} \left( -\frac{v}{2} f^{(4)}(0) dt + dB_2 \right). \quad (H1) $$

Using the solution for the variable $r(t) = e^{-\theta t + vB_1(t)}$, we can integrate and obtain

$$ \begin{aligned} \tilde{T}_{\beta}(t) &:= T_{\beta}(t) - T_{\beta}(0) \\ &= \int_{0}^{t} e^{2vB_{1}(s)+v^{2}f''(0)s} \left[ \frac{\pi c dB_{1}}{6\beta^{2}v} - \frac{c dB_{2}v}{24\pi} \right. \\ &\quad \left. + \left( \frac{cf^{(4)}(0)v^{2}}{48\pi} - \frac{\pi c f''(0)}{12\beta^{2}} \right) ds \right], \end{aligned} \quad (H2) $$

where we subtract the initial value $T_\beta(0) = \pi c / (12\beta^2 v^2)$. We once again reparametrize time as in Appendix D, with $B_i(t-s') = \tilde{B}_i(t) - \tilde{B}_i(s')$. In addition to Eq. (D7), we also need

$$ \int_0^t e^{2vB_1(s)+v^2f''(0)s} dB_1(s) = e^{2v\tilde{B}_1(t)+v^2f''(0)t} \left[ \int_0^t e^{-2v\tilde{B}_1(s)-v^2f''(0)s} d\tilde{B}_1(s) + 2vf^{(2)}(0) \int_0^t e^{-2v\tilde{B}_1(s)-v^2f''(0)s} ds \right], \quad (H3) $$

which leads to the alternative expression

$$ \tilde{T}_{\beta}(t) = e^{2v\tilde{B}_1(t)+v^2f''(0)t} \int_0^t e^{-2v\tilde{B}_1(s)-v^2f''(0)s} \left[ \frac{\pi c d\tilde{B}_1}{6\beta^2 v} - \frac{c d\tilde{B}_2 v}{24\pi} + \frac{\pi c f^{(2)}(0) ds}{4\beta^2} - \frac{c f^{(4)}(0) v^2 ds}{16\pi} \right]. \quad (H4) $$

We can now derive a closed stochastic equation satisfied by $\tilde{T}_\beta(t)$. After some manipulations, we arrive at
$$d\tilde{T}_{\beta} = dB_{1}\left(\frac{\pi c}{6\beta^{2}v} + 2\tilde{T}_{\beta}v\right) - \frac{cdB_{2}v}{24\pi} + \frac{1}{48}dt\left(-\frac{4\pi c f^{(2)}(0)}{\beta^{2}} + \frac{c f^{(4)}(0)v^{2}}{\pi} - 48f^{(2)}(0)\tilde{T}_{\beta}v^{2}\right). \quad (H5)$$

We can further simplify this expression and bring it back to the Bougerol form [73]. First, we put together the two noise terms, introducing a new standard Brownian noise dB:

$$d\tilde{T}_{\beta} = -\frac{\pi c dt f^{(2)}(0)}{12\beta^{2}} + \frac{c dt f^{(4)}(0)v^{2}}{48\pi} - dt f^{(2)}(0)\tilde{T}_{\beta}v^{2} + A(\tilde{T}_{\beta})dB, \quad (H6)$$

$$A(\tilde{T}_{\beta}) = \sqrt{-\frac{c^{2}[16\pi^{4}f^{(2)}(0) - 8\pi^{2}\beta^{2}f^{(4)}(0)v^{2} + \beta^{4}f^{(6)}(0)v^{4}]}{576\pi^{2}\beta^{4}v^{2}} + c\tilde{T}_{\beta}\left(\frac{f^{(4)}(0)v^{2}}{6\pi} - \frac{2\pi f^{(2)}(0)}{3\beta^{2}}\right) - 4f^{(2)}(0)\tilde{T}_{\beta}^{2}v^{2}. \quad (H7)$$

Then, we observe that setting

$$\tilde{T}_{\beta}(t) = \frac{c\tilde{\kappa}_{0}}{4\pi}\left(\frac{\omega_{\beta}(t)}{\omega_{0}} - 1\right) - \frac{\pi c}{12\beta^{2}v^{2}} \Rightarrow T_{\beta}(t) = \frac{c\tilde{\kappa}_{0}}{4\pi}\left(\frac{\omega_{\beta}(t)}{\omega_{0}} - 1\right), \quad (H8)$$

where $\tilde{\kappa}_0$ and $\omega_0$ are defined in Eq. (24), the variable $\omega_\beta(t)$ satisfies exactly the same Langevin equation (25). The only difference lies in the initial condition as $\tilde{T}_\beta(0) = 0$ or, equivalently, Eq. (43).

[1] C. Gross and I. Bloch, *Quantum Simulations with Ultracold Atoms in Optical Lattices*, Science **357**, 995 (2017).

[2] I. Bloch, *Ultracold Quantum Gases in Optical Lattices*, Nat. Phys. **1**, 23 (2005).

[3] U. Schneider, L. Hackermüller, J. P. Ronzheimer, S. Will, S. Braun, T. Best, I. Bloch, E. Demler, S. Mandt, D. Rasch, and A. Rosch, *Fermionic Transport and Out-of-Equilibrium Dynamics in a Homogeneous Hubbard Model with Ultracold Atoms*, Nat. Phys. **8**, 213 (2012).

[4] P. N. Jepsen, J. Amato-Grill, I. Dimitrova, W. W. Ho, E. Demler, and W. Ketterle, *Spin Transport in a Tunable Heisenberg Model Realized with Ultracold Atoms*, Nature (London) **588**, 403 (2020).

[5] F. Alet and N. Laflorencie, *Many-Body Localization: An Introduction and Selected Topics*, C. R. Phys. **19**, 498 (2018).

[6] P. Calabrese, F. H. L. Essler, and G. Mussardo, *Introduction to 'Quantum Integrability in out of Equilibrium Systems', J. Stat. Mech. (2016) 064001*.

[7] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, *Weak Ergodicity Breaking from Quantum Many-Body Scars*, Nat. Phys. **14**, 745 (2018).

[8] M. Serbyn, D. A. Abanin, and Z. Papić, *Quantum Many-Body Scars and Weak Breaking of Ergodicity*, Nat. Phys. **17**, 675 (2021).

[9] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho et al., *Controlling Quantum Many-Body Dynamics in Driven Rydberg Atom Arrays*, Science **371**, 1355 (2021).

[10] P. Sala, T. Rakovszky, R. Verresen, M. Knap, and F. Pollmann, *Ergodicity Breaking Arising from Hilbert Space Fragmentation in Dipole-Conserving Hamiltonians*, Phys. Rev. X **10**, 011047 (2020).

[11] V. Khemani, M. Hermele, and R. Nandkishore, *Localization from Hilbert Space Shattering: From Theory to Physical Realizations*, Phys. Rev. B **101**, 174204 (2020).

[12] J. De Nardis, M. Medenjak, C. Karrasch, and E. Ilievski, *Anomalous Spin Diffusion in One-Dimensional Antiferromagnets*, Phys. Rev. Lett. **123**, 186601 (2019).

[13] S. Gopalakrishnan, R. Vasseur, and B. Ware, *Anomalous Relaxation and the High-Temperature Structure Factor of XXZ Spin Chains*, Proc. Natl. Acad. Sci. U.S.A. **116**, 16250 (2019).

[14] E. Ilievski, J. De Nardis, S. Gopalakrishnan, R. Vasseur, and B. Ware, *Superuniversality of Superdiffusion*, Phys. Rev. X **11**, 031023 (2021).

[15] V. B. Bulchandani, S. Gopalakrishnan, and E. Ilievski, *Superdiffusion in Spin Chains*, J. Stat. Mech. (2021) 084001.

[16] A. Bastianello, A. De Luca, B. Doyon, and J. De Nardis, *Thermalization of a Trapped One-Dimensional Bose Gas via Diffusion*, Phys. Rev. Lett. **125**, 240604 (2020).

[17] A. J. Friedman, S. Gopalakrishnan, and R. Vasseur, *Diffusive Hydrodynamics from Integrability Breaking*, Phys. Rev. B **101**, 180302(R) (2020).

[18] A. Bastianello, A. De Luca, and R. Vasseur, *Hydrodynamics of Weak Integrability Breaking*, J. Stat. Mech. (2021) 114003.

[19] V. Khemani, R. Moessner, and S. Sondhi, *A Brief History of Time Crystals*, arXiv:1910.10745.

[20] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, *Discrete Time Crystals*, Annu. Rev. Condens. Matter Phys. **11**, 467 (2020).

[21] F. Machado, D. V. Else, G. D. Kahanamoku-Meyer, C. Nayak, and N. Y. Yao, *Long-Range Prethermal Phases of Nonequilibrium Matter*, Phys. Rev. X **10**, 011043 (2020).
[22] R. P. Bhatt, J. Kilinc, L. Höcker, and F. Jendrzejewski, *Stochastic Dynamics of a Few Sodium Atoms in Presence of a Cold Potassium Cloud*, Sci. Rep. **12**, 2422 (2022).

[23] C. Gross and I. Bloch, *Quantum Simulations with Ultracold Atoms in Optical Lattices*, Science **357**, 995 (2017).

[24] A. Bermudez, M. A. Martin-Delgado, and D. Porras, *The Localization of Phonons in Ion Traps with Controlled Quantum Disorder*, New J. Phys. **12**, 123016 (2010).

[25] R. Blatt and C. F. Roos, *Quantum Simulations with Trapped Ions*, Nat. Phys. **8**, 277 (2012).

[26] C. Noel, P. Niroula, D. Zhu, A. Risinger, L. Egan, D. Biswas, M. Cetina, A. V. Gorshkov, M. J. Gullans, D. A. Huse, and C. Monroe, *Measurement-Induced Quantum Phases Realized in a Trapped-Ion Quantum Computer*, Nat. Phys. **18**, 760 (2022).

[27] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, *Quantum Entanglement Growth under Random Unitary Dynamics*, Phys. Rev. X **7**, 031016 (2017).

[28] T. Zhou and A. Nahum, *Emergent Statistical Mechanics of Entanglement in Random Unitary Circuits*, Phys. Rev. B **99**, 174205 (2019).

[29] A. Chan, A. De Luca, and J. T. Chalker, *Solution of a Minimal Model for Many-Body Quantum Chaos*, Phys. Rev. X **8**, 041019 (2018).

[30] M. Knap, *Entanglement Production and Information Scrambling in a Noisy Spin System*, Phys. Rev. B **98**, 184416 (2018).

[31] M. Bauer, D. Bernard, and T. Jin, *Stochastic Dissipative Quantum Spin Chains (I): Quantum Fluctuating Discrete Hydrodynamics*, SciPost Phys. **3**, 033 (2017).

[32] D. Bernard and B. Doyon, *Diffusion and Signatures of Localization in Stochastic Conformal Field Theory*, Phys. Rev. Lett. **119**, 110201 (2017).

[33] D. A. Rowlands and A. Lamacraft, *Noisy Coupled Qubits: Operator Spreading and the Fredrickson-Andersen Model*, Phys. Rev. B **98**, 195125 (2018).

[34] D. Bernard, *Can the Macroscopic Fluctuation Theory be Quantized?*, J. Phys. A **54**, 433001 (2021).

[35] D. Bernard and T. Jin, *Open Quantum Symmetric Simple Exclusion Process*, Phys. Rev. Lett. **123**, 080601 (2019).

[36] H.-P. Breuer, F. Petruccione et al., *The Theory of Open Quantum Systems* (Oxford University, New York, 2002).

[37] I. Lesanovsky, M. van Horssen, M. Guťá, and J. P. Garrahan, *Characterization of Dynamical Phase Transitions in Quantum Jump Trajectories beyond the Properties of the Stationary State*, Phys. Rev. Lett. **110**, 150401 (2013).

[38] M. Bauer, D. Bernard, and T. Jin, *Equilibrium Fluctuations in Maximally Noisy Extended Quantum Systems*, SciPost Phys. **6**, 45 (2019).

[39] F. Carollo, R. L. Jack, and J. P. Garrahan, *Unraveling the Large Deviation Statistics of Markovian Open Quantum Systems*, Phys. Rev. Lett. **122**, 130605 (2019).

[40] F. Carollo, J. P. Garrahan, and R. L. Jack, *Large Deviations at Level 2.5 for Markovian Open Quantum Systems: Quantum Jumps and Quantum State Diffusion*, J. Stat. Phys. **184**, 13 (2021).

[41] For clarity, we stress that we are using stationary state in a broad sense as a stationary distribution over pure states. This is a particular unraveling of the density matrix which describes the average dynamics.

[42] D. Bernard and T. Jin, *Solution to the Quantum Symmetric Simple Exclusion Process: The Continuous Case*, Commun. Math. Phys. **384**, 1141 (2021).

[43] B. Skinner, J. Ruhman, and A. Nahum, *Measurement-Induced Phase Transitions in the Dynamics of Entanglement*, Phys. Rev. X **9**, 031009 (2019).

[44] Z.-C. Yang, Y. Li, M. P. A. Fisher, and X. Chen, *Entanglement Phase Transitions in Random Stabilizer Tensor Networks*, Phys. Rev. B **105**, 104306 (2022).

[45] S. Choi, Y. Bao, X.-L. Qi, and E. Altman, *Quantum Error Correction in Scrambling Dynamics and Measurement-Induced Phase Transition*, Phys. Rev. Lett. **125**, 030505 (2020).

[46] D. Bernard, F. H. L. Essler, L. Hruza, and M. Medenjak, *Dynamics of Fluctuations in Quantum Simple Exclusion Processes*, SciPost Phys. **12**, 042 (2022).

[47] B. Shapiro, *Cold Atoms in the Presence of Disorder*, J. Phys. A **45**, 143001 (2012).

[48] M. Ippoliti and V. Khemani, *Postselection-Free Entanglement Dynamics via Spacetime Duality*, Phys. Rev. Lett. **126**, 060501 (2021).

[49] J. Marino and A. Silva, *Relaxation, Prethermalization, and Diffusion in a Noisy Quantum Ising Chain*, Phys. Rev. B **86**, 060408(R) (2012).

[50] J. Marino and A. Silva, *Nonequilibrium Dynamics of a Noisy Quantum Ising Chain: Statistics of Work and Prethermalization after a Sudden Quench of the Transverse Field*, Phys. Rev. B **89**, 024303 (2014).

[51] G. Zhang, E. Novais, and H. U. Baranger, *Rescuing a Quantum Phase Transition with Quantum Noise*, Phys. Rev. Lett. **118**, 050402 (2017).

[52] E. G. Dalla Torre, E. Demler, T. Giamarchi, and E. Altman, *Quantum Critical States and Phase Transitions in the Presence of Non-equilibrium Noise*, Nat. Phys. **6**, 806 (2010).

[53] E. G. Dalla Torre, E. Demler, T. Giamarchi, and E. Altman, *Dynamics and Universality in Noise-Driven Dissipative Systems*, Phys. Rev. B **85**, 184302 (2012).

[54] A. Scheie, N. E. Sherman, M. Dupont, S. E. Nagler, M. B. Stone, G. E. Granroth, J. E. Moore, and D. A. Tennant, *Detection of Kardar-Parisi-Zhang Hydrodynamics in a Quantum Heisenberg Spin-1/2 Chain*, Nat. Phys. **17**, 726 (2021).

[55] X. Wen, Y. Gu, A. Vishwanath, and R. Fan, *Periodically, Quasi-periodically, and Randomly Driven Conformal Field Theories (II): Furstenberg's Theorem and Exceptions to Heating Phases*, SciPost Phys. **13**, 082 (2022).

[56] D. Bernard and P. Le Doussal, *Entanglement Entropy Growth in Stochastic Conformal Field Theory and the KPZ Class*, Europhys. Lett. **131**, 10007 (2020).

[57] P. Calabrese and J. Cardy, *Entanglement Entropy and Conformal Field Theory*, J. Phys. A **42**, 504005 (2009).

[58] A. Belavin, A. Polyakov, and A. Zamolodchikov, *Infinite Conformal Symmetry in Two-Dimensional Quantum Field Theory*, Nucl. Phys. **B241**, 333 (1984).

[59] P. Calabrese and J. Cardy, *Quantum Quenches in Extended Systems*, J. Stat. Mech. (2007) P06008.

[60] P. Calabrese and J. Cardy, *Quantum Quenches in 1 + 1 Dimensional Conformal Field Theories*, J. Stat. Mech. (2016) 064003.
[61] P. Francesco, P. Mathieu, and D. Sénéchal, *Conformal Field Theory* (Springer Science & Business, New York, 2012).

[62] D. Bernard, K. Gawedzki, and A. Kupiainen, *Slow Modes in Passive Advection*, J. Stat. Phys. **90**, 519 (1998).

[63] D. Bernard, *Turbulence for (and by) Amateurs*, arXiv:cond-mat/0007106.

[64] G. Falkovich and K. Gawędzki, and M. Vergassola, *Particles and Fields in Fluid Turbulence*, Rev. Mod. Phys. **73**, 913 (2001).

[65] More precisely, the expectation value employed here is conditioned to the value of the variable $r$ at time $t$.

[66] At a fixed time, the random variable $\kappa(t)$ is uniquely defined by its distribution. However, a stochastic process requires a formal relation between the random variables $\kappa(t)$ and $\kappa(t + dt)$, which corresponds to realizations of the noise defined on two different time domains. In mathematical terms, backward and forward evolution correspond to two different filtrations [67].

[67] I. Karatzas and S. Shreve, *Brownian Motion and Stochastic Calculus* (Springer Science & Business, New York, 2012), Vol. 113.

[68] K. Gawedzki and P. Horvai, *Sticky Behavior of Fluid Particles in the Compressible Kraichnan Model*, J. Stat. Phys. **116**, 1247 (2004).

[69] J. Warren, *Sticky Particles and Stochastic Flows*, arXiv:1409.6946 (2014).

[70] Y. Le Jan and O. Raimond, *Integration of Brownian Vector Fields*, Ann. Probab. **30**, 826 (2002).

[71] G. Barraquand and M. Rychnovsky, *Large Deviations for Sticky Brownian Motions*, J. Probab. **25**, 1 (2020).

[72] P. Bougerol, *Exemples de Théorèmes Locaux sur les Groupes Résolubles*, in *Annales de l'IHP Probabilités et Statistiques* (1983), Vol. 19, pp. 369-391.

[73] S. Vakeroudis, *Bougerol's Identity in Law and Extensions*, arXiv:1201.5823.

[74] G. Gallavotti and E. G. D. Cohen, *Dynamical Ensembles in Nonequilibrium Statistical Mechanics*, Phys. Rev. Lett. **74**, 2694 (1995).

[75] P. Calabrese and A. Lefevre, *Entanglement Spectrum in One-Dimensional Systems*, Phys. Rev. A **78**, 032329 (2008).

[76] H. W. J. Blöte, J. L. Cardy, and M. P. Nightingale, *Conformal Invariance, the Central Charge, and Universal Finite-Size Amplitudes at Criticality*, Phys. Rev. Lett. **56**, 742 (1986).

[77] I. Affleck, *Universal Term in the Free Energy at a Critical Point and the Conformal Anomaly*, in *Finite-Size Scaling, Current Physics—Sources and Comments* Vol. 2, edited by J. L. Cardy (Elsevier, New York, 1988), pp. 347-349.

[78] J. D. Nardis, D. Bernard, and B. Doyon, *Diffusion in Generalized Hydrodynamics and Quasiparticle Scattering*, SciPost Phys. **6**, 49 (2019).

[79] D. Bernard and B. Doyon, *Conformal Field Theory out of Equilibrium: A Review*, J. Stat. Mech. (2016) 064005.

[80] A. De Luca, J. Viti, D. Bernard, and B. Doyon, *Non-equilibrium Thermal Transport in the Quantum Ising Chain*, Phys. Rev. B **88**, 134301 (2013).

[81] T. Giamarchi, *Quantum Physics in One Dimension* (Clarendon, Oxford, 2003), Vol. 121.

[82] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, *Transport in Out-of-Equilibrium XXZ Chains: Exact Profiles of Charges and Currents*, Phys. Rev. Lett. **117**, 207201 (2016).

[83] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, *Emergent Hydrodynamics in Integrable Quantum Systems out of Equilibrium*, Phys. Rev. X **6**, 041065 (2016).

[84] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, *Quantum Generalized Hydrodynamics*, Phys. Rev. Lett. **124**, 140603 (2020).

[85] F. M. Surace, P. P. Mazza, G. Giudici, A. Lerose, A. Gambassi, and M. Dalmonte, *Lattice Gauge Theories and String Dynamics in Rydberg Atom Quantum Simulators*, Phys. Rev. X **10**, 021041 (2020).