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

We consider critical one dimensional quantum systems initially prepared in their groundstate 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 2-pt chiral correlation function can be obtained from solving two coupled ordinary stochastic differential equations. In contrast with the general expectation of heating, we demonstrate that the system reaches a non-trivial and universal stationary state characterized by broad distributions. As an example, we analyse 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 obtain a similar result for the entanglement entropy production associated to a given interval of size . The corresponding stationary distribution has a 3/2 right tail for all , and converges to a one-sided Levy stable for large . Our results are benchmarked via analytical and numerical calculations for a chain of non-interacting spinless fermions with excellent agreement.

Introduction. -The coherent dynamics of macroscopic quantum systems has attracted a lot of interests in the last years, both for its fundamental importance and for its relevance to experimental setups [1][2][3][4]. From the theoretical point of view, recent advances have led to a much deeper level of understanding about the mechanism of equilibration and thermalisation of isolated manybody quantum systems brought out-of-equilibrium. In most cases, the eigenstate thermalisation 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 thermalisation has pinpointed a few phenomena, such as many-body localization (MBL) [5], equilibration towards generalized Gibbs ensembles due to integrability [6], quantum scars [7][8][9] and Hilbert space fragmentation [10,11]. In particular, for integrable systems, anomalous transport has been observed beyond the expected ballistic, with superdiffusive behavior [12][13][14][15]. Nevertheless, weak integrability breakings have been shown to eventually lead to thermalisation [16][17][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 in discrete time has recently also attracted a lot of interest, with several exact results involving random unitary circuits [22][23][24]. Although the finite time dynamics exhibits interesting connections to growth processes, at large time the system relaxes to a trivial infinite temperature ensemble. Other results were obtained the context of stochastic dynamics in continuous time [25][26][27], on the Fredrickson-Andersen model [28], and the quantum simple symmetric exclusion process (QSSEP) [29,30]. In these case, by looking at the average dynamics over noise realisations, one obtains an effective Lindbladian description, which once again admits only the infinite temperature state as stationary point. It is not obvious however that the average dynamics is representative of the typical one. Recent studies have thus explored the statistical behavior beyond the average dynamics [26], including large deviations [31][32][33][34]. In particular, for the finite QSSEP with periodic boundary conditions, it was shown that the stationary state is uniformly distributed among Gaussian states with the same occupancy of the initial state [32]. In the presence of appropriate open boundary conditions, a non-trivial steady state can be attained in the QSSEP [30,35]. Another remarkable mechanism involves the use of quantum measurements which compete with the inner unitary dynamics of the system to produce non-trivial stationary states visible in the statistics of entanglement [36]. In general, an important question for quantum stochastic dynamics is whether non-trivial stationary distribution can emerge when the thermodynamic limit is taken before the large time one. In this direction, the finite time fluctuations were studied for QSSEP in [37].

In this paper, we consider critical one dimensional quantum systems initially prepared in their groundstate. 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 [38][39][40], in particular investigating their stability properties under 1/f noise [41,42]. Here, we introduce a spatially-smooth white noise at t > 0 coupled to the energy density and bring the system out-ofequilibrium by evolving it under the corresponding unitary dynamics. Here, we show that the full distribution of correlation functions reaches a non-trivial stationary limit, not visible at the level of noise averages which instead exhibit apparent heating. By using conformal field theory, we deduce a universal description of the outof-equilibrium dynamics. Following Ref. [43], the noise coupled to the energy density is interpreted as a random metric in the CFT formulation. As a result, the dynamics of chiral operators can be solved in terms of stochastic trajectories. In [43] the focus was 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. By studying in detail the stochastic trajectories we show that any 2-pt chiral correlation function can be obtained from solving the two coupled ordinary stochastic differential equations in (7). Remarkably, this leads to a stationary state characterized by broad distributions. We analyse in particular the local energy density. We are able to obtain the first moment at all time and show that, as a consequence of conformal invariance, it diverges exponentially fast in time; nevertheless, the average is not indicative of the typical behavior: the stationary distribution reaches a simple and universal form with a fat tail with 3/2 decay exponent (see (16)) and thus no finite integer moments. We obtain a similar result for the entanglement entropy production associated to a given interval of size . The corresponding stationary distribution still has a 3/2 right tail for all , and converges to a one-sided Levy stable for large , see (23), whose physical origin can be related to the return probability of a Brownian process.

To test the universality of this theory, we study analytically and numerically a chain of non-interacting spinless fermions at half-filling. 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 diverges and the CFT predictions are recovered, as confirmed numerically by computing the local energy and entanglement entropy on the lattice.

Model and CFT. -We consider a 1d model at a second order quantum phase transition, initially prepared in the groundstate |0 of its gapless Hamiltonian Ĥ0 . To simplify the notation we will indicate the groundstate averages as . . . = Ψ 0 | . . . |Ψ 0 . To simplify, we will assume continuous space, but the treatment can be readily extended to lattice systems. At time t = 0, a perturbation Ĥ1 is turned on, by coupling a space-dependent noise term with the system energy density. The total Hamiltonian then takes the form

Ĥ = Ĥ0 + Ĥ1 = ˆdx(1 + η(x, t)) ĥ(x) (1)
where ĥ(x) is the hamiltonian density and η(x, t) is the noise characterised by the space-time correlation η(x, t)η(x , t ) = δ(t -t )f (x -x ). The function f (x) parameterises the noise correlation, and has the dimension of a (turnover) time; it is even and has positive Fourier transform and we choose it smooth, monotonously decreasing for x > 0, with a fast decay at infinity when x 1. In the following, we will indicate as Ō the average of any quantity O over the noise realisations. An effective low-energy description of Ĥ0 can be obtained using the scale invariance which holds at the secondorder critical point. This in turns implies an emergent conformal symmetry of the unperturbed theory [44]. A powerful description can then be obtained by means of conformal field theory (CFT), which has been successfully employed even in out-of-equilibrium dynamics and quantum quenches [45,46]. In particular, all local operators splits into chiral (∼ right moving) and antichiral (∼ left moving) components organised into families [START_REF] Francesco | Conformal field theory[END_REF]. All operators within each chiral family descend from a primary field φ± (x) with given conformal dimension ∆ ± . In particular, the hamiltonian density can be represented as ĥ(x) ∼ v( T + (x) + T -(x)) where v is the light velocity and T + and Tare the two components of the stressenergy tensor, with T + (x) -T -(x) the momentum density. This implies that under Ĥ0 , chiral primary fields simply translate in time φ± (x, t) = φ(x ∓ vt). In [43], it has been recently shown that the time evolution of primary fields under Ĥ in (1) can be interpreted as a conformal transformation. In practice, one introduces the stochastic trajectories q ± (s) as solution of the Langevin equation

dq ± (s) ds = ±v(1 + η(q ± (s), s)) (2)
where the Ito convention is assumed. We define the functions X ± t (y) as the initial condition for (2) at t = 0 (i.e. q ± (0) = X ± t (y)) such that q ± (t) = y. Then, the evolution of a primary field is simply given by φ± (y, t) = (X ± t (y)) ∆ ± φ(X ± t (y), 0) and arbitrary correlation functions at time t can be reduced to those in the initial state via

n i=1 φ+ i (y i , t) = J (y 1 , . . . , y n ) n i=1 φ+ i (X + t (y i ))
(3) where the factor J (y 1 , . . . , y n ) = i (X + t (y i )) ∆ + i accounts for the jacobian of the conformal transformation. An analogous equation can be written for the antichiral component with X + t → X - t and ∆ + i → ∆ - i . To study the correlation functions in (3) and their sample to sample fluctuations, we thus need the joint probability distribution function (jpdf) of the set of 2n random variables X ± t (y j ), j = 1, . . . , n. Let us first focus on the jpdf of X ± t (y j ) = x j for a fixed chirality, choosing either ±, which we denote P ± t (x|y). Given n trajectories q ± j (s) satisfying Eq. ( 2) with endpoints q ± j (0) = x j , q ± j (t) = y j , P ± t (x|y) is thus the jpdf of the initial positions x = (x 1 , . . . , x n ) of these n trajectories conditioned to the positions of their final point y = (y 1 , . . . , y n ). As shown in [48], it satisfies the Fokker-Planck (FP) equation also studied in the context of turbulence and passive scalar [START_REF] Bernard | [END_REF][50][51]

∂ t P ± t (x|y) = (±v n i=1 ∂ i + v 2 2 n i,j=1 ∂ i ∂ j f (x i -x j ))P ± t (x|y) (4
) where ∂ i = ∂/∂x i . Since all trajectories are evolving according to Eq. ( 2) within the same realisation of the noise, they are correlated, which appears as an interaction in (4). The martingale property from the initial time implies x i = y i ∓ vt; additionally the trajectories cannot cross one another, so that the coordinates y and x are always ordered in the same way.

Two-point correlations. -Consider a primary field Φ(x, t) = φ+ (x, t) × φ-(x, t) with a scaling dimension ∆ = ∆ + + ∆ -. According to (3), the time evolution of the 2-point correlation function satisfies

C(y 1 , y 2 , t) ≡ a 2∆ 0 Φ(y 1 , t) Φ(y 2 , t) = = C(y 1 , y 2 , t = 0)e -∆+κ + -∆-κ - (5)
where we assume y 1 > y 2 and we set

κ ± (y 1 , y 2 , t) = ln (X ± t (y 1 ) -X ± t (y 2 )) 2 (y 1 -y 2 ) 2 X ± t (y 1 )X ± t (y 2 ) (6)
This expression gives access to the full statistics of the correlation function. We first focus on either κ + or κ -. Indeed the trajectories corresponding to the chiral and antichiral components are typically separated by a distance ∼ 2vt. Therefore, at large time for 2vt 1, the noises they feel become uncorrelated, and we expect the two components to decorrelate (this is discussed in more details below). Let us define := y 1 -y 2 > 0 and the ratio r = (X ± t (y 1 )-X ± t (y 2 ))/ . Using spatial homogeneity the one point pdf of κ = κ ± (y 1 , y 2 , t) is only a function of and t and 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, κ) of κ and r.

By considering four trajectories, i.e. n = 4 in (4), one first obtains a FP equation for the jpdf of r, X ± t (y 1 ) = x 1 and X ± t (y 2 ) = x 2 . This is achieved through the linear change of variable r = (x 1 -x 2 )/ , x 1 = (x 3 -x 1 )/ , x 2 = (x 4 -x 1 )/ , and taking the limit → 0 (the center of mass variable x 1 + x 2 decouples). Remarkably, performing another change of variable, one finds that the random variables r and κ = log(r 2 /(x 1 x 2 )) are solution [START_REF]Solution of (7) provides the distribution for κ in (5) at a fixed time t, but the two quantities can be different as stochastic processes in t[END_REF] of the Stochastic differential equation (SDE) in Ito's form [48]

dr = v dW 1 (t) , dκ = v 2 g(r)dt + v dW 2 (t) (7)
where v 2 g(r) is a drift term and dW 1 (t), dW 2 (t) are two centered Gaussian white noises in time of r-dependent variances dW 1 (t) 2 = 2A(r)dt, dW 2 (t) 2 = 2C(r)dt and cross correlation dW 1 (t)dW 2 (t) = B(r)dt [53]. Here, we introduced

A (r) = f (0) -f ( r) 2 , B (r) = 2f ( r) + 4(f (0) -f ( r)) 2 r , C (r) = 4(f (0) -f ( r) + rf ( r)) 2 r 2 -f (0) -f ( r) , g (r) = -f (0) - 2(f (0) -f ( r)) 2 r 2 . ( 8
)
Eqs. ( 7) must be solved with the initial condition r = 1 and κ = 0 at t = 0. Since the equation for r does not involve κ, one may first solve for r(t) and then insert the solution for r(t) in the equation for κ(t). For finite , Eq. ( 7) cannot be solved explicitly for an arbitrary f (x). Nevertheless, one can understand its behavior at finite time in two regimes = y 1 -y 2 1 and = y 1 -y 2 1, as well as in the large time limit.

Small separation 1.

-In this case, Taylor expanding the function f in (8) one finds the leading behavior at small of each function as

A (r) - 1 2 f (0)r 2 B (r) 2 6 f (4) (0)r 3 (9)
C (r) -

4 72 f (6) (0)r 4 g (r) 2 12
r 2 f (4) (0) (10) We see that (7) can be rewritten, after a redefinition of the noises dW 1 (t) = rdB 1 (t) and dW 2 (t) = -2 6 r 2 dB 2 (t), in the form

dr = rvdB 1 (t) , dκ = - 2 6 r 2 v(- v 2 f (4) (0)dt + dB 2 (t)) (11)
where now dB 1 (t), dB 2 (t) are r-independent Gaussian white noises with fixed correlation matrix dB

1 (t) 2 = -f (0)dt, dB 2 (t) 2 = -f (6) (0)dt, 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 1. One can first solve the stochastic equation for r and obtain after an application of Ito's lemma: r(t) = e -θt+vB1(t) . We have defined θ = -v 2 f (0)/2 > 0, which is the inverse of a characteristic time. Alternatively, one can change variable to ρ = log r, which obeys the stochastic equation dρ = -θdt + vdB 1 (t), implying that ρ(t) is a Brownian motion with a drift. Hence P t (r) is a log-normal distribution for the variable r, with ln r = -θt , Var[ln r] = 2θt .

Since θ > 0, this shows, interestingly, that the trajectories X t (y 1 ), X t (y 2 ) tend to get closer as time grows, a manifestation of the phenomenon of sticky particles observed in turbulent fluids [START_REF] Gawedzki | [END_REF][55][56][57]. On the other hand, r = 1 holds independently of time, which shows that although the typical value of r decreases to zero, the distribution of r is broadening with time. Hence higher moments, such as r 2 , grow with time.

Using this result it is easy to calculate the noise average of κ = κ ± by simply averaging (11) using that dB 2 = 0. Evaluating r 2 e 2θt and integrating over time one finds

κ = v 2 f (4) (0) 2 24θ (e 2θt -1) + O( 4 ) (13)
This result allows to obtain the noise average of ln C(y 1 , y 2 , t), irrespective of the possible correlations between κ + and κ -, by averaging the logarithm of (5). It is possible to calculate the higher integer moments κ n and one finds [48] that they all grow exponentially with time within the validity of the small regime. However, upon averaging (7) over the noise, we observe that κ ≤ 2θt is an exact bound since g(r) ≤ g(∞) = -f (0). Thus, while the moments are still diverging at large time, the exponential growth is only valid when 2 e 2θt 1.

Nevertheless, as we now show, the pdf of κ remarkably converges to a stationary distribution at large time. This distribution is very broad and consistently does not possess any integer moments for n ≥ 1. To obtain the pdf of κ, we proceed in two steps: first we plug the solution r(t) = e -θt+vB1(t) into the equation (11) for κ; secondly, we use time-reversal [48] to recast the resulting stochastic equation in the standard form studied by Bougerol [58,[START_REF] Vakeroudis | Bougerol's identity in law and extensions[END_REF]. This equation is best expressed with a change of variable from

κ to Y κ = 2 κ0 sinh Y ω 0 -1 , κ0 = - 1 12
f (4) (0) f (0) ( 14
)
where we defined ω 0 = f (6) (0)f (0)

f (4) (0) 2 -1 -1/2 > 0. This leads to dY = -2θ tanh Y dt + √ 8θd B(t) (15)
where B(t) is a standard Wiener process with d B(t) 2 = dt. Eq. ( 15) describes the Langevin motion of a particle in a confining potential U (y) = 2θ log cosh y 2θ|y| at temperature 4θ. Hence it reaches an (equilibrium) stationary measure at large time,

P stat (Y ) = C/ cosh(Y ) with C = √ 2π/Γ(1/4) 2 . Thus, in the limit → 0, the scaled variable ω ≡ sinh Y = ω 0 (1 + κ/( 2 κ0 )) is de- scribed by the stationary distribution B(ω) ≡ C (1 + ω 2 ) 3/4 (16)
Consistently, the power-law tail ∝ |ω| -3/2 implies that ω (and κ as well) does not have finite integer moments. We remark that, because of the time-reversal transformation, the stochastic process (15) for Y and the original one for κ in (11) are not equivalent: the former is ergodic in time, while the latter has a finite (random) limit κ(t → ∞). Nevertheless, they are equivalent in law at fixed t.

Large interval

1.

-The leading behavior at large of each function up to O(1/ 2 ) reads

A (r) f (0) 2 B (r) 4 f (0) 2 r (17) C (r) -f (0) + 4 f (0) 2 r 2 g (r) -f (0) - 2f (0) 2 r 2
where we assume that f (x) decays faster than a power law. At leading order one can set B (r) 0 which implies that the equation ( 7) for κ becomes independent of r.

Using g (r) -f (0) and C (r) -f (0) to leading order one obtains

κ = 2θt + √ 4θW (t) , r = 1 + 2f (0)v W (t) (18)
where W (t), W (t) are two uncorrelated standard Brownian motion with dW (t) 2 = d W (t) 2 = dt. Note that the growth of κ saturates the exact bound κ < 2θt. To obtain this result we have assumed not only that 1 but also that r 1. Although this condition holds for finite time, since r is undergoing diffusion, we see from ( 18) that for t ∼ 2 /f (0), r(t) may become close to zero and the condition is violated.

Large-time limit. -To search for a stationary distribution for any , we derive in [48] an evolution equation for the characteristic function of

κ, Q k (r 0 , t) = e -ikκ r0
, where the superscript r 0 = r(t = 0) indicates the initial condition for the variable r in (7), ultimately setting r 0 = 1. Looking for a time dependent solution in the large time limit Q k (r 0 , t) → Q k (r 0 ), we obtain that

Q k (r) = 2 r 2 f (0) 2(f ( r) -f (0)) ik G k (r) . ( 19
)
where G k (r) satisfies the Schrodinger-like equation for r ≥ 0

-G k (r) -k(k + i)V (r)G k (r) = 0 (20)
with boundary conditions are G k (0) = 1 and lim r→+∞ G k (r) = 0 and the potential

V (r) = - d 2 dr 2 log[f (0) -f ( r)] + 2 f (0) f (0) -f ( r) . (21)
Studying this equation (see [48]), we find that the stationary distribution P stat (κ) (obtained by Fourier inversion of Q k (r = 1)) depends non-trivially on and f (x), with however a symmetry valid in all cases

P stat (-κ 0 + β) P stat (-κ 0 -β) = e β , κ 0 = -log 2(f ( ) -f (0)) 2 f (0) (22)
reminiscent of a fluctuation theorem [START_REF] Gallavotti | [END_REF]. Beyond the asymptotic solution at small given in (16), one can also

derive at large → +∞, κ = θ 2 χ/f (0) + O( ) with χ distributed according to L(χ) ≡ 1 √ 2π e -1 2χ χ 3/2 Θ(χ) (23)
i.e. to the stable one sided Levy distribution of index 1/2. Eqs. ( 16) and ( 23) characterise asymptotic behavior of the stationary distribution for κ at small and large respectively. For intermediate values of , an explicit expression is not available for generic f (x). However we provide some cases which are explicitly solvable for any , e.g. for f (x) = 1/ cosh n x with n = 1, 2.

One can show that for sufficiently smooth f (x),

Q k (r) = 1 + O( √ k) at small k, irrespectively of .
This translates onto a -3/2 power-law tail on the positive κ side. For the left tail κ → -∞, it follows from Eq. ( 22) that P stat (κ) decays exponentially for any for > 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 < , log r is approximately Brownian with a negative drift (see (12)). In this regime, κ remains almost constant as dκ ∝ r 2 dt. However, r starts from 1 and has a finite probability to move right towards r > ; in this case, κ increases by ∼ ∆t (see Eq. ( 18)) where ∆t is the time interval before r(t) hits again . Since in this regime, r(t) is approximately an unbiased Brownian, ∆t is the corresponding first-passage time, distributed as 1/(∆t) 3/2 .

Finally, let us recall for time t 1/(2v) the two chiral components κ + and κ -decorrelate hence their joint distribution reaches in the large time limit the factorized form P stat (κ + )P stat (κ -). We expect that κ ± (t) defined in Eq. ( 6) as a stochastic process in time has P stat (κ ± ) as the one-time ergodic measure. This is not in contradiction with the fact that κ defined in the auxiliary stochastic system (7) has a (random) finite limit κ ± (t → ∞).

Entanglement production. -An interesting application of the above results is the calculation of the entanglement entropies. Let us define as ρ A,t the reduced density matrix for the interval A = [y 1 , y 2 ] at time t. Then, the Renyi entropies are defined as S

(n) t = 1/(1 -n) ln Tr ρ n A,t .
Introducing the twist fields Φn (y, t) [61], one can identify Tr ρ n A,t ∝ Φn (y 1 , t) Φn (y 2 , t) . We can thus express the entropy production

S (n) t ≡ S (n) t -S (n) t=0 using (5) as S (n) t = (n + 1)c 24n (κ + + κ -) (24)
where we used that for twist fields

∆ ± = c(n -1/n)/24.
The Von Neumann entropy corresponds to n = 1 and we denote it simply as S t . Eq. (24) shows that all the Renyi 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 [62], it implies that the entanglement spectrum retain the same form in each noise realization, i.e. the density of log λ/ log λ max is independent of the noise, while log λ max = -c/24(κ + +κ -).

Using that the noise average S (n) t = (n+1)c 12n κ, we see that at early times two different growth regimes exists: for large intervals ( 1) the average entropy production grows linearly with time as in (18), while for small intervals it grows as in (13). Finally, at large time, we find that the entropy production (24) reaches a stationary distribution given up to a scale by the convolution P stat * P stat determined above, still with a -3/2 power law tail and no integer moments.

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 T + (x) and T -(x). Their time evolution can be obtained either by direct calculation of the commutators [ Ĥ, T ± (x)] in the Heisenberg equation [48]. 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 (for simplicity, we omit the ± superscript)

T (y, t) = |X t (y)| 2 T (X t (y), t = 0) - c 24π (S • X t )(y) (25)
where the second term is proportional to the central change c and the Schwarzian derivative

(S • X t )(y) = 2X (y)X (y) -3X (y) 2 2X (y) 2 (26)
Since in the initial state T (y, 0) = 0, one has simply

T (y, t) = - c 24π (S.X t )(y) (27)
The time evolution of this quantity can be simply obtained from the above analysis of κ ± in the regime → 0. Indeed, expanding (6) in small y 1 -y 2 we see that

lim →0 1 2 κ(y 2 + , y 2 , t) = - 1 6 (S.X t )(y 2 ) . (28)
Eq. ( 28) is well-known in the CFT context: it reflects the occurence of the conformal anomaly in the transport of the stress-tensor in a gaussian free field theory, see [START_REF] Francesco | Conformal field theory[END_REF]. Here, it implies that the distribution T ± can be obtained from the one of κ ± in the limit of small . In particular, 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) ≡ ĥ(x, t) = v( T + + T -), one has h(x, t) = lim →0 vS t /(π 2 ). Additionally, denoting the noise average (which is space independent) as e(t) ≡ h(x, t), one has from ( 13)

e(t) = cv 3 f (4) (0) 48πθ e 2θt -1 . (29)
FIG. 1. Top: The average energy eF (τ ) vs t = τ /ξ, the scaled time, for different values of the noise correlation length ξ. Continuous line are obtained from the numerical solution of the Wigner function equation (34), while the markers corresponds to the exact dynamics of (31) for L = 2048. The dotted dashed line is the CFT result (29), which from ( 35) is predicted to hold for large ξ. Bottom: The median of the distribution of ξ 2 ( ĥi τ -ĥi 0) vs the scaled time τ /ξ. In the limit of ξ → ∞, the median is expected to decrease towards the negative asymptotic value predicted by CFT (dot dashed horizontal line). For finite ξ, the median starts to grow at large times, suggesting that heating may eventually dominate on the lattice.

For t 1/(2v) we expect T ± to decorrelate; thus from (14), we extract the exact stationary distribution of the one-point energy density

lim t→∞ h(x, t) in law = vcκ 0 4π (Ω/ω 0 -2) ( 30
)
where Ω = ω + + ω -and ω ± are independent random variables both distributed according to B(ω) in ( 16).

Comparison with free fermions. -We consider a model of non-interacting spinless fermions, where we denote τ the time in the lattice model with Hamiltonian

ĤF = i (1 + η i (τ )) ĥi , ( 31
) ĥi = -J(â † i âi+1 + â † i+1 âi ) .
Note that the local energy density is only defined up to a total derivative. A good way to fix such an ambiguity and preserve the continuous limit is to impose paritywhich is verified in (31) [63][64][65]. On the contrary, in the presence of a finite chemical potential, a different definition is required (see [48]). In the absence of noise, this corresponds to a dispersion relation (k) = -2J cos(k).

For the noise, we choose the correlation

η i (τ )η j (τ ) = τ 0 δ(τ -τ )F (i -j) , F (j) = f (j/ξ) (32)
where ξ 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. The dynamics induced by ( 31) is better studied in terms of the noise averaged Wigner function n τ (k) = j â † j+j âj τ e ikj where . . . τ denotes the quantum average at time τ under the evolution ĤF in (31). We choose the initial state as the groundstate so that n τ (k) does not depend on the lattice site j. Also,

n τ =0 (k) = Θ(k+k f )-Θ(k-k f )
, where Θ(z) is the Heaviside function and k f is such that (k f ) = 0 and k f = π/2, which corresponds to half filling. The system is critical and can be described with a conformal field theory with central charge c = 1 [START_REF] Giamarchi | Quantum physics in one dimension[END_REF]. Using the Wigner function, we can express the noise averaged energy density with respect to the groundstate

e F (τ ) ≡ ĥi τ -ĥi 0 = ˆdk 2π (k)(n τ (k) -n 0 (k)) (33)
One can derive [48] an exact evolution equation for n τ (k) which reads

∂ τ n τ (k) = τ 0 ˆdk 2π F (k ) (k + k /2) 2 (n τ (k + k ) -n τ (k)) (34)
The study of this equation [48] shows that around the Fermi points n τ (k) takes the scaling form n τ (k) = n((k f ∓ k)ξ, t = τ /ξ). In turns this leads to the result for fixed t as ξ → +∞ lim ξ→∞ ξ 2 e F (tξ) = e(t) (35) where e(t) is given in Eq. ( 29). Hence the CFT predicts, upon rescaling, the mean energy for the fermion system. This explicit calculation suggests that in the scaling limit of large ξ, the noisy dynamics in ( 31) is fully captured by the universal description provided by the CFT, upon rescaling space and time as j = xξ, τ = tξ and set τ 0 = ξ.

In order to validate this hypothesis, we have computed numerically the two-point correlation matrix C ij (τ ) ≡ a † i a j τ . Since the model ( 31) is non-interacting and the initial groundstate is Gaussian, for each realisation of the noise, all quantities can be expressed via the Wick theorem in terms of the coefficients C ij . In Fig. 1 top, we show the convergence for ξ → ∞ of the noise-averaged energy density e F (τ ) to its CFT prediction, consistently with (35). We observe the emergence of a characteristic time scale τ * (ξ), diverging with ξ, after which the CFT (1) (τ ) vs t = τ /ξ evolving in time under (31), for increasing values of ξ and fixed ratio F /ξ = = 1/2. The dotted line is obtained from the numerical solution of Eqs. (7) and using (24) for n → 1, i.e. SF (τ ) = cκ/6 with c = 1. Middle: Distribution of the scaled energy density ξ 2 ĥi(τ ) at ξ = 64.. For the analytical prediction, we used that in the scaling limit ξ 2 ĥj → h(x, t) and at large time h(x, t) is distributed as (30), which in the present case reduces to h in law = (3Ω -5)/(24π). In the inset, the right tail of the distribution is shown in log-log scale, showing the predicted ∝ h -3/2 tail. Right: Distribution of the entanglement entropy SF (τ ) at ξ = 64. for an interval of size F = 32. For the analytical prediction, we used that in the stationary limit limτ→∞ SF (τ ) in law = c(κ + + κ -)/12, with κ ± independently distributed according to Pstat(κ). The stationary distribution Pstat(κ) is obtained by numerically inverting its Fourier transform Q k (1) for = 1/2, defined in (19) (without any prefactor) with = F /ξ. Once again, a good agreement is found and corrections emerge at τ τ * (ξ).

ξ =16, ℓ F =8 ξ =32, ℓ F =16 ξ =64, ℓ F =32 CFT 0 1 2 3
Our prediction is that the CFT describes also the distribution of these quantities. We show in Fig. 2 middle, the one-point PDF for the local energy density ξ 2 ĥi τ at the largest time τ available. As one sees it compares reasonably well, with no free parameters, with the prediction from the CFT, i.e. the convolution P stat * P stat where P stat was obtained in (16). This confirms that with the chosen f (x), at this observation time, 2v F τ /ξ is large enough, so that the two chiral components have decoupled. As we see on Fig. 1 top, the average energy grows with time, consistent with P stat having infinite first moment. The median of the energy distribution, shown in Fig. 1 bottom, is thus a better probe of the typical behavior. Remarkably, it is found to decrease with time, approaching at large ξ a stationary value compatible with the CFT prediction e median stat = -c 2π κ0 < 0, see Eq. ( 14). However, for any finite ξ, we expect lattice ef-fects to eventually break the CFT description. Although we have not studied it in detail, it is expected that the ultraviolet cutoff induces heating and stationarity does not hold anymore, as hinted by the rebounce in the median observed at larger times. Finally, in Fig. 2 right, as a representative of the finite-behavior, we compare the distribution of the entanglement entropy at different times for intervals of size F = ξ, and = 1/2, with the analytic prediction obtained from CFT.

Conclusion.

-We have identified an out-ofequilibrium protocol which leads to a non-trivial stationary state for generic gapless one-dimensional system.

Several questions and directions remain open. First, 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. 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 non-interacting systems, in particular to observe the role played by interactions.

Finally, the surprising existence of a stationary state 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 correspond-ing lattice system, similarly to what has been observed in the context of many-body quantum scars [START_REF] Surace | [END_REF]. More generally, it remains an open question whether lattice effects or the presence of finite correlation length are compatible with the emergence of non-trivial steady states in the thermodynamic limit. In this section we derive the Fokker-Planck equation ( 4) of the text, for the joint PDF of the backward stochastic trajectories x 1 = X + t (y 1 ), . . . , x n = X + t (y n ) associated to the Langevin equation ( 2). We thus consider only a given chirality, here we choose +, but the same Fokker-Planck equation holds for the chirality -, with v → -v. We do not consider here the joint PDF of both chiralities. So here we denote simply X + t → X t . One can show that Eq. ( 2) translates into a stochastic equation for the variable X t (y i ) as a function of t which takes the form (see Eq. (58)[START_REF] Vakeroudis | Bougerol's identity in law and extensions[END_REF] in Supp. Mat. of [43])

dX t (y i ) = v 2 2 X t (y i )f (0)dt -X t (y i )v dt + dW t (y i ) (S1.1)
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 Eq. ( 1) and (2) via W t (y) = ˆt 0 ds η(y, s) , dW t (y)dW t (y ) = dtf (y -y ) (S1.2)

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 (2), we define:

g t (y 1 . . . y n ) = G X t (y 1 ) . . . X t (y n ) (S1.3)
By Ito calculus the time variation dg t = g t+dt -g t of this observable is obtained by expanding up to second order.

dg t = n j=1 ∂ j G X t (y 1 ), . . . , X t (y n ) dX t (y j ) + 1 2 n j,m=1 ∂ j ∂ m G X t (y 1 ), . . . , X t (y n ) dX t (y j )dX t (y m ) (S1.4)
Here we shortened the notation setting

∂ xj . . . ∂ xm G(x 1 . . . x n ) = ∂ j . . . ∂ m G(x 1 . . . x n )
. Using (S1.1) and (S1.2), 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 ) . (S1.5)
Averaging over the noise we obtain

dg t dt = -v n j=1 ∂ j G X t (y 1 ), . . . , X t (y n ) X t (y j ) + v 2 f (0) 2 n j=1 ∂ j G X t (y 1 ), . . . , X t (y n ) X t (y j )+ + v 2 2 n j,m=1 f (y j -y m )∂ j ∂ m G X t (y 1 ), . . . , X t (y n ) X t (y j )X t (y m ) (S1.6)
From the chain rule for the derivation with respect to the variables {y i } it is easy to check that

∂ yj ∂ ym g t = X t (y j )X t (y m )∂ j ∂ m G + δ j,m X t (y j )∂ j G (S1.7)
which finally leads to:

dg t dt = -v n j=1 ∂ yj + v 2 2 n j,m=1 f (y j -y m )∂ yj ∂ ym g t . (S1.8)
It is useful to re-express (S1.8) 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 ω(y) as q j • ω(y) = y j ω(y) and p j • ω(y) = -ı∂ j ω(y). For these conjugate variables the canonical quantization holds [q i , p j ] = ıδ i,j . We can then define

H FP ≡ -iv n i=1 p i + v 2 2 ij
p i p j f (q i -q j ) , (S1.9) so that we can rewrite

dg t dt = -H † FP • g t , H † FP ≡ iv n i=1 p i + v 2 2 ij f (q i -q j )p i p j . (S1.10)
To deduce the Fokker-Planck equation for P (x|y), defined in the text, we need one more step. From the definition of P (x|y), we have

g t (y) = ˆdx G(x )P t (x |y) (S1.11)
Since G is an arbitrary function we choose G(x ) = δ(x -x) to recover g t (y) → P (x|y) the jpdf for the initial points x = (x 1 , . . . , x n ) of the stochastic trajectories and finally deduce from (S1.10)

∂ t P t (x|y) = -H T FP [y] • P t (x|y) (S1.12)
This equation can be formally solved starting from the initial condition at t = 0 is P t=0 (x|y) = δ(x -y). It is useful to employ the bra y| and ket |x notation for the eigenstates of the position operators qj . Then, we can represent the probability distribution as:

P t (x|y) ≡ y|e -tH † FP |x = x|e -tHFP |y (S1.13)
where the last equality follows from the hermitian conjugation. Therefore, the following equation must also holds:

∂ t P t (x|y) = -H FP [x] • P t (x|y) (S1.14)
where the action of the differential operator H FP [x] is now over the variables x. More explicitly using (S1.9), we arrive at

∂ t P t (x|y) = v   n i=1 ∂ i + v 2 n i,j=1 ∂ i ∂ j f (x i -x j )   P t (x|y) . (S1.15)
which coincides with Eq. ( 4) given in the main text.

Appendix 2: Joint evolution of r and κ

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

Derivation of the evolution equation for

Pt(r, x 1 , x 2 )
We will 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 , x 2 defined in the text. The starting point is the application of the Fokker-Planck equation (S1.15) in the case of 4 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

x 1 = X ± t (y 1 ), x 2 = X ± t (y 2 ), x 3 = X ± t (y 1 + ), x 4 = X ± t (y 2 + ) (S2.1)
We then perform the change of variables x 1 , x 2 , x 1 = (x 3 -x 1 )/ , x 2 = (x 4 -x 2 )/ together with the limit= → 0. By applying this change of variables in the Fokker-Planck equation above and taking the limit → 0, one obtains the following equation for P t = P t (x 1 , x 2 , x 1 , x 2 ):

∂ t P t = v (∂ x1 + ∂ x2 ) + v 2 f (0)(∂ 2 x1 + ∂ 2 x2 ) + v∂ x1 ∂ x2 f (x 1 -x 2 ) - v 2 F (0)(∂ 2 x 1 (x 1 ) 2 + ∂ 2 x 2 (x 2 ) 2 ) +v∂ x2 ∂ x 1 x 1 f (x 1 -x 2 ) -v∂ x1 ∂ x 2 x 2 f (x 1 -x 2 ) -v∂ x 1 ∂ x 2 x 1 x 2 f (x 1 -x 2 ) P t (S2.2)
Then, we exploit the invariance under translation by making the change of variables from (

x 1 , x 2 , x 1 , x 2 ) to (R = x1+x2 2 , r = (x 1 -x 2 )/ , x 1 , x 2 )
. Finally, integrating out the center of mass variable R, we obtain

∂ t P = v 2 1 2 ∂ 2 r (f (0)-f (r ))- f (0) 2 (∂ 2 x 1 (x 2 ) 2 +∂ 2 x 1 (x 2 ) 2 )-f (r )∂ x 1 ∂ x 2 x 1 x 2 - 1 ∂ r f (r )(∂ x 1 x 1 +∂ x 2 x 2 ) P (S2.3)
From this result we now obtain the Fokker-Planck equation for the joint distribution P t (κ, r) of the variable κ and the variable r defined as

P t (κ, r) ≡ ˆdx 1 dx 2 P (r, x 1 , x 2 )δ κ -ln r 2 x 1 x 2 . (S2.4)
Interestingly, P t (κ, r) satisfies a closed evolution equation. To obtain it, we compute the time derivative of (S2.4) using (S2.3), then we integrate by part and obtain

2 v 2 ∂ t P = ∂ 2 r (f (0) -f ( r)) + 4 f (0) -f ( r) + rf ( r) r 2 ∂ 2 κ + 2 f (0)(∂ κ -∂ 2 κ ) -2 f ( r)∂ 2 κ +2 ∂ r ∂ κ f ( r) + 4∂ r (f (0) -f ( r)) 1 r ∂ κ + 2 f (0) -f ( r) r 2 ∂ κ P . (S2.5)
Then, it is easy to see that the FP equation in Eq. (S2.5) is equivalent to the stochastic equations for r and κ (in Ito convention) given in (7) of the main text. 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 (4). It provides an alternative way to derive Eq. ( 7) 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

dx i (t) = ∓vdt + dW i (t) , dW i (t)dW j (t) t = f (x i (t) -x j (t))dt (S2.6)
where the W i (t)'s are mutually correlated Wiener processes in time and shall not be confused with W t (y) introduced in (S1.2). Formally, the expectation value O t used here is conditioned to the position of the trajectories x i (t) at time t. The ∓ refers to ± chiralities, but the drift term plays no role in the following, so will be 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. ( 2).

Let us set = 1 and restore later. One performs the linear change of variable x 3 = x 1 + x 1 , x 4 = x 2 + x 2 and x 1 -x 2 = r, which leads to (we keep the time dependence implicit for convenience)

dr = dW 1 -dW 2 , dx 1 = 1 (dW 3 -dW 1 ) , dx 2 = 1 (dW 4 -dW 2 ) (S2.7)
The correlations can be expressed in terms of r, x 1 , x 2 (the center of mass decouples) and performing the limit → 0 one finds with j = 1, 2

dr 2 = 2(f (0) -f (r))dt , drdx j = -f (r)x j dt , (dx j ) 2 = -f (0)(x j ) 2 dt , dx 1 dx 2 = -f (r)x 1 x 2 dt (S2.8)
Now we have for κ = log(r 2 /x 1 x 2 ) using Ito's rule

dκ = 2 dr r - dx 1 x 1 - dx 2 x 2 + g(r)dt (S2.9)
where g(r)dt is the Ito drift which can be computed using (S2.8)

- 1 r 2 dr 2 + 1 2 (dx 1 ) 2 (x 1 ) 2 + 1 2 (dx 2 ) 2 (x 2 ) 2 = g(r)dt , g(r) = 2 r 2 (f (0) -f (r)) -f (0) (S2.10)
and depends only on r. We also need the correlations of the noise part

dκdr = 2 dr 2 r - drdx 1 x 1 - drdx 2 x 2 = B(r)dt , B(r) = 4 r (f (0) -f (r)) + 2f (r) (S2.11) dκ 2 = 4 dr 2 r 2 + (dx 1 ) 2 (x 1 ) 2 + (dx 2 ) 2 (x 2 ) 2 -4 drdx 1 rx 1 -4 drdx 2 rx 2 +2 dx 1 dx 2 x 1 x 2 = 2C(r)dt , C(r) = 4 f (0) -f (r) + rf (r) r 2 -f (0)-f (r) (S2.
12) to which we must add dr 2 = 2A(r) with A(r) = f (0) -f (r). Restoring one finds the equations in the text.

To solve the equation, we then discretize time dt → ∆t and define

∆W 1 = 2A(r)β 1 , ∆W 2 = B(r) 2 2A(r) β 1 + 2C(r) - B(r) 2 2A(r) β 2 (S3.3)
where β 1 , β 2 are independently gaussian random variables with zero average and ∆t variance. To rewrite the equation in a more numerically stable way, we extract from the quantities A(r), B(r), their leading behavior at small r Ã(r) = A(r)/r 2 , B(r) = B(r)/r 3 . (S3.4) so that Ã(r), B(r) are both finite in the limit r → 0. So we finally have the discrete evolution equations

ρ(t + ∆t) = ρ(t) + v 2 Ã(r)∆tβ 1 -v 2 Ã(r)∆t (S3.5) κ(t + ∆t) = κ(t) + v 2 g(r)∆t + vr 2 B(r) 2 2 Ã(r) ∆tβ 1 + 2C(r) - r 4 B(r) 2 2 Ã(r) ∆tβ 2 (S3.6)
Appendix 4: Analysis of the short distance regime for κ

In the limit 1, one has Eq. ( 11) in the main text that we report for convenience

dr = rvdB 1 , dκ = - 2 6 r 2 v(- v 2 f (4) (0)dt + dB 2 ) (S4.1)
dB 1 dB 1 = -f (0)dt , dB 2 dB 2 = -f (6) (0)dt , dB 1 dB 2 = -f (4) (0)dt (S4.2)
We will now show from these equations, that κ 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 vB 1 (t) + 1 2 v 2 f (0)t (S4 .3)
Injecting this solution in the equation for κ and integrating in time we arrive at

κ(t) = - 2 6 ˆt 0 e 2vB1(s)+v 2 f (0)s v - v 2 f (4) (0)ds + dB 2 (s) (S4.4)
This equation gives already a closed representation for κ(t).

We can further simplify Eq. (S4.4) by making use of the time-reversal symmetry.

We denote s = t -s and introduce new processes as d Bi (s ) = dB i (t -s ). Therefore for both i = 1, 2, we have

B i (t -s ) = ˆt-s 0 dB i (s ) = ˆt 0 dB i (s ) + ˆt-s t dB i (s ) = B i (t) - ˆs 0 d Bi (s ) = Bi (t) -Bi (s ) (S4.5)
where Bi (s) are equivalent Wiener processes which move in the opposite direction. We want now to rewrite Eq. (S4.4) in terms of the reversed processes Bi . A little bit of care is needed for the stochastic integral. Indeed, writing explicitly the Ito's integral, we have

ˆt 0 e 2vB1(s)+v 2 f (0)s dB 2 (s) = lim n→∞ i e 2vB1(si)+v 2 f (0)si (B 2 (s i+1 ) -B 2 (s i )) = = e 2v B1(t)+v 2 f (0)t lim n→∞ j e -2v B1(sj )-v 2 f (0)sj ( B2 (s j ) -B2 (s j-1 )) (S4.6)
where the s i 's are a partition of [0, t] and sj = t -s i defines an equivalent partition, with j = N -i. Clearly, the last expression does not converge to a stochastic integral in the Ito form. We thus expand the exponent arriving at ˆt 0 e 2vB1(s)+v 2 f (0)s dB 2 (s) = e 2v B1(t)+v 2 f (0)t ˆt 0 e -2vB1(s)-v 2 f (0)s d B2 (s) + 2vf (4) (0) ˆt 0 e -2vB1(s)-v 2 f (0)s ds (S4.7)

Applying these transformations to (S4.4) we arrive at

κ(t) = - 2 6 e 2v B1(t)+v 2 f (0)t κ1 ˆt 0 e -2v B1(s)-v 2 f (0)s v 3v 2 f (4) (0)ds + d B2 (s) κ2 (S4.8)
Using Ito's lemma, dκ = κ 2 dκ 1 + κ 1 dκ 2 + dκ 1 dκ 2 , which, after collecting different contributions, leads to

dκ = v 2 2 12 f (4) (0) -f (0)κ dt - 2 v 6 d B2 (t) + 2vd B1 (t)κ (S4.9)
The correlations of the d Bj (t) are the same as the ones of the dB j (t) in Eq. (S4.2). However, κ(t) defined by Eq. (S4.8) is a different process in t than κ(t) defined by Eq. (S4.4). At fixed t the two random variables have the same law, but as t is varied the trajectories are different (since the relation between Bi and B i involves t explicitly). As a consequence, the two stochastic equations (S4.9) and (S4.4) are inequivalent although they lead to the same singletime distribution for κ(t). One illustration of that is that while the second process converges, i.e. κ(t) → κ ∞ , where the distribution of κ ∞ is given below, the first process is ergodic (with the same law) as we show below. We can recast (S4.9) as an equation with a single Brownian process dB (we are using that

a 1 dB 1 + a 2 dB 2 = -a 2 1 f (0) -a 2 2 f (4) (0) -2a 1 a 2 f (6) (0)d B, where d B is a new Wiener process with standard normalization, d B2 = dt) dκ = v 2 2 12 f (4) (0) -f (0)κ dt + vd B 6 -4 f (6) (0) -24 2 f (4) (0)κ + 144f (0)κ 2 (S4.10)
With the change of variable κ = κ 0 (ω/ω 0 -1) where we defined as in the main text ( 14)

κ 0 = - 2 12
f (4) (0) f (0) , ω 0 = 1 f (6) (0)f (0) f (4) (0) 2 -1 (S4.11)
Note that ω 0 is real and positive as it is guaranteed by the positivity of the Fourier transform f (k) > 0 of f (x). Indeed,

f (6) (0)f (0) -f (4) (0) 2 = ˆk2 f (k) ˆk6 f (k) -ˆk4 f (k) 2 > 0 (S4.12)
which is a consequence of the Cauchy-Schwartz inequality.

The SDE for ω becomes

dω = 2θωdt + √ 8θ 1 + ω 2 dW (S4.13)
Accounting for the initial condition ω(t = 0) = ω 0 , the solution of this SDE takes the form

ω(t) = e √ 8θB(t)-2θt ω 0 + √ 8θ ˆt 0 e - √ 8θB(s)+2θs dγ s (S4.14)
which is then a Bougerol variable with drifted Brownian motion in the exponent [START_REF] Vakeroudis | Bougerol's identity in law and extensions[END_REF]. It is useful to do the change of variable (S4. 13)

ω = sinh Y , Y = Y (ω) = argsinhω (S4.15)
One has

Y (ω) = 1 √ 1 + ω 2 = 1 cosh Y , Y (ω) = - ω (1 + ω 2 ) 3/2 (S4.16)
and from Ito's rule, it follows Eq. ( 15) in the main text.

with the potential

V (r) = 2 f ( r) 2 + (f (0) -f ( r)) (f ( r) + f (0)) (f (0) -f ( r)) 2 = - d 2 dr 2 log[f (0) -f ( r)] + 2 f (0) f (0) -f ( r) (S5.12)
on the positive half-space r ≥ 0 with the boundary condition

G k (0) = 1 , lim r→∞ G k (r) = 0 . (S5.13)
Before solving various cases, let us indicate a nice symmetry property. One notes that the dependence in k in (S5.11) is only via the prefactor γ ≡ k(k + i) of the potential: G k (r) ≡ G(γ; r). This implies that ˆdκP stat (κ)e -ikκ = φ k (1)G k (1) = e ikκ0 G(γ; 1) (S5. 14) where κ 0 is explicitly

κ 0 = 1 ik log φ k (1) = -log 2(f ( ) -f (0)) 2 f (0) (S5.15)
Hence, from Eq. (S5.8)

P stat (κ) = ˆdk 2π e ik(κ+κ0) G(k(k + i); 1) = e κ+κ 0 2 ˆdu 2π e iu(κ+κ0) G(u 2 + 1 4 ; 1) (S5.16)
where in the last equality we changed variable k = u -i/2. It implies in particular that Here, we show that for a smooth noise correlation f (x), the right tail of the stationary distribution is always a power-law P stat (κ) ∝ κ -3/2 for κ → +∞, independently of . Thanks to (S5.1), it is enough to prove the following expansion at small k for its Fourier transform

P stat (κ) = e
Q k (r = 1) = 1 + C √ k + O(k). (S5.19)
with C is a constant (see below). To prove Eq. (S5.19), we proceed as follow. First of all, since φ k (r) is analytic in k, using (S5.9), we can focus on the small k expansion of G k (r). The small k behavior of the solution of Eq. (S5.11) can be obtained setting x = r √ γ, with γ = k(k + i). Then, setting G k (r) = g(r √ γ), we can rewrite (S5.11) in the limit k → 0 as

-g (x) -V ∞ g(x) = 0 ⇒ g(x) = e - √ -V∞x (S5.20)
where we set V ∞ = lim r→∞ V (r) and enforced the boundary conditions (S5.13). This implies

G k (r) ∼ g(r √ γ) = e - √ -V∞γr , ∀r = O(γ -1/2 ) (S5.21)
This is still not enough because we require the expansion (S5.19) for r = 1. However, fixing δ > 0 and r, we can write

|G k (r) -G k (δ/ √ γ)| = ˆr δ/ √ γ dr G k (r ) ≤ γ ˆr δ/ √ γ dr |V (r )G k (r )| ≤ Kγ(r -δ/ √ γ) (S5.22)
where we set K = sup r |V (r)G k (r)|, which is finite for a sufficiently smooth f (x) ∈ C 6 . As a consequence, at small γ,

G k (r) = G k (δ/ √ γ) + O(δ √ γ) = --V ∞ γ (g(δ) + O(δ)) δ→0 -→ --V ∞ γ + O(γ) (S5.23)
Finally, integrating over r

G k (1) = 1 + ˆ1 0 dr G k (r) = 1 --V ∞ γ + O(γ) (S5.24)
which proves Eq. (S5. 19) with C = √ -V ∞ .

Small limit

At small one finds that the potential is a harmonic oscillator

V (r) = 4 r 2 f (4) (0) 2 -f (6) (0)f (0) 36f (0) 2 + O( 6 r 4 ) (S5.25)
and

κ 0 = κ0 2 + O( 4 ) , κ0 = - f (4) (0) 12f (0) (S5.26)
We see that κ 0 is O( 2 ) as → 0, so the random variable κ = O( 2 ) in this limit. Therefore, we can obtain a independent limit by scaling k = k/ 2 . In terms of this variable the potential term in the Schrodinger equation has a finite limit

t -k(k + i)V (r) 4 k2 r 2 κ2 0 ω 2 0 (S5.27)
where we have defined as in the text

ω 0 = 1 f (0)f (6) (0) f (4) (0) 2 -1 (S5.28)
A solution of Eq. ( 20) with the potential (S5.25) and the boundary conditions (S5.13) can be expressed in terms of Bessel function as

G k (r) = 2 3/4 Γ(1/4) κ0 | k| ω 0 1/4 √ rK 1 4 | k|r 2 κ0 ω 0 (S5.29)
and the prefactor has been fixed imposing that G k (r = 0) = 1. This leads to

Q k (r = 1) = 2 3/4 Γ(1/4) κ0 | k| ω 0 1/4 e i kκ0 K 1 4 | k|κ 0 ω 0 (S5.30)
which allows to determine P stat (κ) by Fourier inversion from (S5.8).

One can check that this coincides with the result in the text, identifying κ = κ0 ω ω0 -1 . Equivalently, we obtain the scaling form

P stat (κ) 1 1 2 P ( κ 2 ) , P (κ) ≡ Cω 0 κ0 1 + ω 2 0 κ + κ0 κ0 2 -3/4 (S5.31)
and one can check that the Fourier transform of P ˆdk

2π e i kκ P (κ) = Q k (r = 1) (S5.32)
as given in (S5.30). This can be seen using the identity

ˆdxe ikx 1 (1 + x 2 ) 3/4 = √ 2π (2|k|) 1/4 K 1 4 (|k|) Γ 3 4 (S5.33)

Large limit

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

-k(k + i)V (r) -k(k + i) 2 f (0) f (0) ∼ 2i kθ f (0) . (S5.34) f (x) V (r)/ 2 f (k) Smoothness e -|x| 1 (1-e -r ) 2 2 1+k 2 C 0 2e -|x| -e -2|x| -2 1-e -r 12 k 4 +5k 2 +4 C 2 1 cosh x -tanh( r) 2 π cosh(kπ/2) C ∞ 1 (cosh x) 2 -2 tanh( r) 2 πk sinh(kπ/2) C ∞ TABLE S1.
A few examples of noise correlation functions f (x) leading to Schrödinger equations with a solvable potential V (r).

where we used once again the scaling k = k/ 2 , which implies κ = 2 κ, but with → ∞ in this case. From this potential, we immediately derive the solution respecting the boundary conditions (S5.13) in the form

G k (r) = e - √ 2i kθ/f (0)r (S5.35)
Note that at large

κ 0 →∞ = log( 2 (-f (0))/2f (0)) (S5.36)
so that kκ 0 →∞ -→ 0 and therefore Q k (r) = G k (r). Inverting the Fourier transform (S5.8), we obtain once again the stationary distribution

P stat (κ) 1 1 2 p κ 2 , p(κ) ≡ θ 2πf (0) e -θ 2f (0)κ κ3/2 Θ(κ) (S5.37)
Equivalently, denoting κ = θ 2 ξ/f (0) one finds that ξ is distributed according to L(ξ) in Eq. (23) in the text, i.e. the stable one sided Levy distribution of index 1/2.

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 S1, we list a few interesting cases. Here, we focus on the case

f (x) = 1/ cosh(x) (S5.38)
which is analytic and fastly decaying. Setting γ = k(k + i), the solution G k (r) respecting the boundary conditions (S5.13) can be expressed in terms of hypergeometric function

G k (r) = e -√ γ r (1 + tanh( r)) √ γ 2 F 1 √ γ 2 -1 2 γ + 1 4 + 1 4 , √ γ 2 + 1 2 γ + 1 4 + 1 4 ; √ γ + 1; 1 cosh( r) 2 2 F 1 √ γ 2 -1 2 γ + 1 4 + 1 4 , √ γ 2 + 1 2 γ + 1 4 + 1 4 ; √ γ + 1; 1 (S5.39)
Equivalently, this expression can be represented in terms of generalised Legendre functions

G k (r) = P - √ γ 1 2 ( √ 1+4γ-1) (tanh( r)) P - √ γ 1 2 ( √ 1+4γ-1) (0) (S5.40)
As a first check, we verify that the solution Eq. (S5.39) reproduces the known solutions in the small/large limits.

Asymptotic limits → ∞

At large , we simply have lim

→∞ Q k / 2 (r) = lim →∞ G k / 2 (r) = e -r √ ik (S5.41)
Consider the particular case of Ô = T ± (y), the stress energy tensor at the point y. We have for the equal-time commutator

[ T + (x), T + (y)] = -ı(2δ (x -y) T + (y) -δ(x -y) T + (y) - c 24π δ (x -y)) (S6.5) [ T -(x), T -(y)] = ı(2δ (x -y) T -(y) -δ(x -y) T -(y) - c 24π δ (x -y)) (S6.6)
From this we deduce

[d Ĥ, T ± (y)] = ±2ı T ± (y)vdW t (y) ± ı(dt + dW t (y))v T ± (y) ∓ ı c 24π vdW t (y) (S6.7)
where dW t (x) = ∂ x dW t (x) (the noise is smooth in space) and similarly for the higher derivatives. This implies the evolution equation for the operators T ± (y, t), which reads

d T ± (y, t) = ∓(2 T ± (y, t)vdW t (y) + (dt + dW t (y))v T ± (y, t) - c 24π vdW t (y))+ + cv 2 f (4) (0) 48π dt -f (0)v 2 T ± (y, t)dt + 1 2 f (0)v 2 ∂ 2 y T ± (y, t)dt . (S6.8)
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 T ± (y, t), which is still a quantum operator. Its evolution is obtained by taking the noise average of (S6.8) and reads

∂ t T ± (y, t) = cv 2 f (4) (0) 48π -f (0)v 2 T ± (y, t) + 1 2 f (0)v 2 ∂ 2 y T ± (y, t) ∓ v∂ y T (y, t) (S6.9)
Another observable is the quantum expectation T ± (y, t) ≡ Ψ 0 | T ± (y, t)|Ψ on any translational invariant state |Ψ in a given noise realisation. It satisfies a stochastic differential equations obtained by taking the quantum expectation of (S6.8) (which leads to a similar equation as (S6.8) but now for a scalar T ± (y, t) ). Note that it describes the coupled stochastic evolution of the two fields T ± (y, t) , a complicated problem. Here we will only focus on (i) noise moments, which have a solvable dynamics (ii) the one-point PDF of T ± (y, t) . 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 ± superscript and denote T (t) ≡ T ± (y, t) the corresponding random variable. Then, one can show that its PDF can be obtained from the stochastic equation (S6.10)

dT = 2T vdB 1 (t) - c 24π vdB 2 (t) + cv 2 f (4) (0) 48π dt -f (0)v 2 T dt (S6.10)
with dB 1 (t) = ∓dW (y, t) and dB 2 (t) = ∓dW (y, t). This can be seen intuitively: indeed, we expect that the spatial derivative terms in (S6.8) are irrelevant since the one point PDF of T ≡ T ± (y, t) does not depend on y. More formally, one can prove that (S6.10) and (S6.8) lead to the same evolution equation for the noise average Z(T ) for any smooth function Z. Although Eq. (S6.10) was obtained here by a completely different method, one can check that Eq. (S6.10) is equivalent to (S4.9) with the correspondence .11) noting also that the noise satisfies (S4.2).

T = lim →0 c 4π 2 κ , ( S6
For the translationally invariant initial state chosen here, it is easy to obtain from (S6.10) the recursive equation for the moments (over the noise) of the quantum expectation of T as

∂ t T n = - v 2 c 2 f (6) (0)n(n -1) 1152π 2 T n-2 + cv 2 f (4) (0)(4n -3)n 48π T n-1 -n(2n -1)v 2 f (0)T n ; (S6.12)
These can be solved with T (t = 0) = 0, which corresponds to |Ψ = |Ψ 0 being the ground state. This leads to the first two moments

T = c f 4 (0) 48πf (0) 1 -e -v 2 f (0)t (S6.13) T 2 = c 2 3456π 2 1 2 f (4) (0) f (0) 2 (e -6v 2 f (0)t -6e -v 2 f (0)t + 5) + f (6) (0) f (0) (e -6v 2 f (0)t -1) (S6.14)
As a result, a single correlator survives. This is expected since the full Hamiltonian (S7.13) is quadratic, so the correlation matrix â † j âj 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

∂ τ n τ (k) = - τ 0 2 m,j,j e ikm 2(ζ j ζ j + ζ * -j ζ * -j ) (S7.9) ×(F (j) -F (m + j)) + (ζ j ζ * -j + ζ * -j ζ j )(F (0) + F (j + j ) -F (m) -F (m + j + j )) â † m+j+j â0 τ
Using that ζ(k) * = j e ikj ζ * -j we finally obtain

∂ τ n τ (k) = τ 0 ˆq F (q)|ζ(-k) + ζ(-q -k) * | 2 (n(k + q) -n(k)) (S7.10)
where

F (q) = j e iqj F (j) (S7.11)
Consider first half-filling, i.e. µ = 0, as in the main text. To recover the coupling with the noise defined there, we must choose ζ(k) = -Je ik . With this choice, |ζ(-k) + ζ(-q -k) * | 2 = (k + q/2) 2 (S7.12)

and one recovers (34) given in the main text. For finite -2 < µ < 2, the Hamiltonian in (S7.1) remains critical and described by a c = 1 CFT. In order to preserve parity and time-reversal symmetry, a simple choice is ζ(k) = (k)/2 with (k) = -J(2 cos k -µ), which corresponds to the coupling to the noise of the form Ĥ = Ĥ0 + 1 2 j η j (τ )( ĥj + ĥj-1 ) (S7. 13) In this case, we obtain after replacing ζ(k) → (k)/2 in (S7.10) ∂ τ n τ (k) = τ 0 4 ˆq F (q)( (k) + (q + k)) 2 (n(k + q) -n(k)) . (S7.14)

Scaling limit

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

n τ (k) n(ξ(k f -k), τ /ξ) + n(ξ(k + k f ), τ /ξ) (S7.15)
The initial condition which corresponds to the ground state reads

n τ =0 (k) = Θ(k + k f ) -Θ(k -k f ) = -1 + Θ(k + k f ) + Θ(k f -k) (S7.16)
which gives for the initial scaling function n(p, t = 0) = -1/2+Θ(p). We can now derive directly an evolution equation for the scaling function n(p, t). In order to do this we replace in Eq. (S7. and similarly for e - F (t) with k ∼ -k f . We now prove the validity of Eq. (35) in the main text. We show that ẽ+ (t) satisfies a a first order differential equation. Indeed, differentiating (S7.18) w.r.t. the variable t and using (S7.17), we obtain

∂ t ẽ+ (t) = -v 3 ˆ∞ -∞ dp 2π
dq 2π p f (p-q) q + p 2 2 (n(q; t)-n(p; t)) = -v 3 2 ˆ∞ -∞ dq 2π dp 2π f (p-q)(p-q) q + p 2

2

[n(q; t)-n(p; t)] (S7.20) where in the last equality we symmetrized the integrand with respect to the exchange p ↔ q. The integral in (S7.20) 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=-∞ = 0 to arrive at ∂ t ẽ+ (t) = v 3 ˆdq 2π du 2π f (u)u (q + u/2) 3 6 (n (q; t) -n (q + u; t)) (S7.21)

where n (q; t) = ∂ q n(q, t). We can now split the integral into two finite terms and replace q → q -u in the second integral. This leads to ∂ t ẽ+ (t) = v 3 ˆdq 2π du 2π f (u)u 1 6 (q + u/2) 3 -(q -u/2) 3 n (q; t) = = v 3 ˆdu 2π f (u) u 4 48π ˆdq n (q; t) + v 3 2 ˆdq 2π q 2 n (q; t)

ˆdu 2π f (u)u 2 = v 3 f (4) (0) 48π -v 2 f (0)ẽ + (τ ) (S7.22)
where we used that ˆdq n (q; t) = 1 , ˆdq 2π q 2 n (q; t) = ˆdq 2π q 2 [n (q; t) -n (q; 0)] = 2 v ẽ+ (t) (S7.23) ˆdu Note that in (S7.23) we used the initial condition (S7. 16) which shows that n (q, 0) = δ(q). We can thus solve (S7.22) with the initial condition ẽ+ (t = 0) = 0 and obtain ẽ+ (t) = f (4) (0)v 1 -e -v 2 f (0)t 48πf (0) (S7.25)

Summing also the equivalent contribution from ẽ-(t) and using the definition of θ = -v 2 f (0)/2, we recover (29) setting the central charge c = 1, which is expected from universality. Note that the present calculation is from first principles. In principle the full scaling function n(p; t) should be universal, i.e. independent of the microscopic details of the discrete model, but depending on f (x), but we have not attempted to obtain it analytically.

To better characterize the time evolution of the scaling function n(p; t), it is also possible to look at the time evolution of the moments defined by and clearly one has from (S7.23) M 0 (t) = 1 and M 2 (t) = 2ẽ + (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. 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) and this provides an estimation for the breakdown of the CFT description at the level of average quantities. Indeed, the Fermi points of k ∼ ±k F independently broaden with time, up to a time τ * (ξ), where their width ∼ ξ √ M 2 is comparable with their initial distance ∼ 2k F . Using M 2 (t) = 2ẽ + (t)/v and (S7.

FIG. 2 .

 2 FIG.2. Left: Noise-averaged Von Neumann entanglement entropy production S(1) (τ ) vs t = τ /ξ evolving in time under(31), for increasing values of ξ and fixed ratio F /ξ = = 1/2. The dotted line is obtained from the numerical solution of Eqs.(7) and using(24) for n → 1, i.e. SF (τ ) = cκ/6 with c = 1. Middle: Distribution of the scaled energy density ξ 2 ĥi(τ ) at ξ = 64.. For the analytical prediction, we used that in the scaling limit ξ 2 ĥj → h(x, t) and at large time h(x, t) is distributed as(30), which in the present case reduces to h

  and (S5.39) of [48]. All simulations are performed on systems of total length L = 2048 and are repeated for N sample = 800 samples. description breaks down. A simple argument based on the broadening of the Wigner function around the Fermi points suggests τ * (ξ) ∝ ξ ln ξ [48], which is in agreement with our numerics. The correlation matrix C can also be used to compute explicitly the Renyi entropies for any interval I of size F . Indeed, setting C I as the F × F matrix obtained restricting the indexes of C(τ ) to I, we have for the Von Neumann entropy of the interval I: S F (τ ) ≡ -Tr[C I ln C I + (1 -C I ) ln(1 -C I )]. In Fig. 2 left, we show the noise averaged entanglement entropy production S F (τ ) ≡ S F (τ ) -S F (0) for the fermion system for intervals of various sizes F on the lattice. Our prediction is that it should equal at large ξ the CFT value S (1) t

2. 2 .

 2 Direct derivation of the stochastic equations for r and κ

κ+κ 0 2 P 5 . 3 .

 253 (|κ + κ 0 |) , P (y) = ˆdu 2π e iuy G(u 2 + 1 4 ; 1) (S5.17)or equivalently the symmetry (reminiscent of a Nishimori condition, or a Galavotti-Cohen theorem)P stat (-κ 0 + β) P stat (-κ 0 -β) = e β (S5.18) Proof of the right 3/2-tail

2π f (u)u 2 =

 2 -f (0) , ˆdu 2π f (u)u 4 = f (4) (0) . (S7.24)

M

  n (t) ≡ ˆ∞ -∞ dp 2π p n n (p, t) = -n ˆ∞ -∞ dp 2π p n-1 (n(p, t) -n(p, 0)) . (S7.26)

∂

  t M n (t) = k α k M n-2k (t) .(S7.27)

25 ) 2 F ξ 2 ⇒

 2522 This leads to ẽ+ (τ * /ξ) = vk

  14) ζ(k) with the energy dispersion relation (k), which gives Eq.(34) in the text. Next, we inject the scaling form Eq. (S7.15) around each Fermi point (choosing +k f ) and replace k = k f -p/ξ, k = -p /ξ. From (32) and (S7.11), we obtain that at large ξ, F (k) ξ f (ξk), which leads to (using τ 0 = ξ)We can now compute the mean energy density with respect to the ground state e F (τ ), which is expanded as e F (τ = tξ) =

		ẽ+ (t) ≡ lim ξ→∞	ˆπ -π ξ ˆ∞ -∞ dk 2π dp (k)[n τ (k) -n 0 (k)] 2π (k f -p ξ )[n(p, t) -n(p, 0)] = -v 1 ξ 2 (ẽ + (t) + ẽ-(t)) ˆ∞ -∞ dp 2π	p [n(p, t) -n(p, 0)]	(S7.18) (S7.19)
	∂ t n(p; t) v 2	ξ 2 4 ˆ∞ -∞ ˆπξ -πξ dp 2π	dp 2π f (p ) p + f (p )	p 2	k f -	p ξ	+ k f -	p + p ξ	ξ→∞ -→

2

(n(p + p ; t) -n(p; t))

2 (n(p + p ; t) -n(p; t)) = v 2 ˆ∞ -∞ dq 2π f (p -q) q + p 2 2

(n(q; t) -n(p; t)) (S7.17)

Acknowledgments. -PLD acknowledges support from the ANR grant ANR-17-CE30-0027-01 RaMaTraF. DB acknowledges support from the project 'ESQuisses', ANR-20-CE47-0014-01.

Supplementary Material

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

In this supplementary material we provide additional details about the calculations in the Letter.

CONTENTS

References 8

Appendix 3: Numerical solution of Eq. (7) In this section, we provide some details about how to numerically solve the stochastic equations for κ and r. The main difficulty arises because r can become typically very small (see Eq. ( 12)). It is thus useful to rewrite such a system of SDE in terms of another variable: ρ = ln r , r = e ρ (S3.1)

Then, using Ito's lemma, the stochastic equation for ρ takes the form

Appendix 5: Stationary measure for κ for any finite

Backward method

We now apply the backward method to the full stochastic equation for r(t), κ(t) in (7). We define

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

expanding with Ito's lemma and averaging we arrive at

with the boundary conditions

and we recall that

The second condition in (S5.4) comes from the fact that dκ = 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 κ is pure diffusion at large r 0 . Let us denote Q k (r) the stationary solution of (S5.3). It thus satisfies

with the boundary conditions

From this stationary solution one obtains the stationary measure P stat (κ) for κ by Fourier inversion

Schrodinger equation for the stationary measure

We can further simplify Eqs. (S5.4) and (S5.6) by removing the first derivative term. This can be achieved setting

then one finds that G p (r) satisfies the Schrodinger equation

in agreement with (S5.35) (θ = 1/2 for (S5.38)). The limit in (S5.41) can be easily obtained using (S5.39) using that

irrespectively of x.

Small check

In this case, the limit is less trivial as k = k / 2 becomes large in the limit of small . 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

Then, we use that lim

to recover, after some manipulations, Eq. (S5.30).

Appendix 6: Distribution of the stress-energy tensor

As we saw in the main text, the distribution of the stress energy tensor can be deduced from the one of κ in the limit of small , simply using (28). An alternative approach is to compute explicitly the time evolution of the stress energy tensor directly from the expression of the Hamiltonian.

Direct derivation

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

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

and define the time evolution U t operator up to time t with the equation

The evolution of an operator Ô(t) ≡ U † t ÔU t takes the form

Note that because of the Ito's convention we need to keep terms up to the double commutators. Here, we will derive the evolution equation Eq. ( 34) 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 the noise

where we have introduced a chemical potential µ (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

where ζ(k) = j ζ j e ikj is the Fourier transform of ζ j and consider the coupling

where the correlation of the noise η j (τ ) is defined in the main text. We will keep ζ(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 will focus on the noise average Wigner function which is defined as Although we are interested in the groundstate of Ĥ0 , these considerations apply to any translational invariant gaussian initial state |Ψ 0 . In this case, because of the noise average, Eq. (S7.4) is independent of the position j.

We now consider the quantum evolution of τ . It is easy to verify that, after the noise average, one obtains the Lindblad form

We can use (S7.4) and Eq. (S7.5) to obtain an evolution equation for n τ (k). Note that the first term in (S7.5) does not contribute because of translational invariance and (S7.2). Using (S7.2) we obtain explicitly