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 $\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$. 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 arXiv:2110.15303v1 [cond-mat.stat-mech] 28 Oct 2021 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 pertur-bationĤ 1 is turned on, by coupling a space-dependent noise term with the system energy density. The total Hamiltonian then takes the form whereĥ(x) is the hamiltonian density and η(x, t) is the noise characterised by the space-time correlation . 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 [47]. All operators within each chiral family descend from a primary fieldφ ± (x) with given conformal dimension ∆ ± . In particular, the hamiltonian density can be represented where v is the light velocity andT + andT − are the two components of the stressenergy tensor, withT + (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 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 (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 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 [49][50][51] ∂ t P ± t (x|y) = (±v 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 where we assume y 1 > y 2 and we set 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 [52] 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 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.
-In this case, Taylor expanding the function f in (8) one finds the leading behavior at small of each function as 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 (11) where now dB 1 (t), dB 2 (t) are r-independent Gaussian white noises with fixed correlation matrix dB 1 (t) 2 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 [54][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 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,59]. This equation is best expressed with a change of variable from κ to Y whereB(t) is a standard Wiener process with dB(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, 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 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 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 κ, , 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 where G k (r) satisfies the Schrodinger-like equation for with boundary conditions are G k (0) = 1 and lim r→+∞ G k (r) = 0 and the potential 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 [60]. Beyond the asymptotic solution at small given in (16), one can also derive at large → +∞, 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 . We can thus express the entropy production S 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 tensorT + (x) andT − (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) where the second term is proportional to the central change c and the Schwarzian derivative Since in the initial state T (y, 0) = 0, one has simply 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 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 [47]. 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 quan- . Additionally, denoting the noise average (which is space independent) as e(t) ≡ h(x, t), one has from (13)  (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 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 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 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, 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 [66]. Using the Wigner function, we can express the noise averaged energy density with respect to the groundstate One can derive [48] an exact evolution equation for n τ (k) which reads 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) 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 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 (without any prefactor) with = F /ξ. Once again, a good agreement is found and corrections emerge at τ τ * (ξ).
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 [67]. 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.

Supplementary Material
Universal out-of-equilibrium dynamics of 1D critical quantum systems perturbed by noise coupled to energy 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)(59) in Supp. Mat. of [43]) 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.
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 Averaging over the noise we obtain From the chain rule for the derivation with respect to the variables {y i } it is easy to check that which finally leads to: 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 so that we can rewrite (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 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) 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 operatorsq j . Then, we can represent the probability distribution as: where the last equality follows from the hermitian conjugation. Therefore, the following equation must also holds: · 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 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
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 We then perform the change of variables 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 ): 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 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 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 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.

Direct derivation of the stochastic equations for r and κ
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 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) 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 Now we have for κ = log(r 2 /x 1 x 2 ) using Ito's rule where g(r)dt is the Ito drift which can be computed using (S2.8) and depends only on r. We also need the correlations of the noise part (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. 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 To solve the equation, we then discretize time dt → ∆t and define 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 so thatÃ(r),B(r) are both finite in the limit r → 0. So we finally have the discrete evolution equations 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 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 Injecting this solution in the equation for κ and integrating in time we arrive at 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 dB i (s ) = dB i (t − s ). Therefore for both i = 1, 2, we have whereB i (s) are equivalent Wiener processes which move in the opposite direction. We want now to rewrite Eq. (S4.4) in terms of the reversed processesB i . A little bit of care is needed for the stochastic integral. Indeed, writing explicitly the Ito's integral, we havê Applying these transformations to (S4.4) we arrive at Using Ito's lemma, dκ = κ 2 dκ 1 + κ 1 dκ 2 + dκ 1 dκ 2 , which, after collecting different contributions, leads to The correlations of the dB j (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 betweenB i 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 With the change of variable κ = κ 0 (ω/ω 0 − 1) where we defined as in the main text (14) Note that ω 0 is real and positive as it is guaranteed by the positivity of the Fourier transformf (k) > 0 of f (x). Indeed, which is a consequence of the Cauchy-Schwartz inequality. The SDE for ω becomes Accounting for the initial condition ω(t = 0) = ω 0 , the solution of this SDE takes the form which is then a Bougerol variable with drifted Brownian motion in the exponent [59]. It is useful to do the change of variable (S4.13) One has 16) and from Ito's rule, it follows Eq. (15) in the main text.
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 Q k (r 0 = 0) = 1 , Q k (r 0 → +∞) = 0 (S5.7) 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 If one chooses φ k (r) so that then one finds that G p (r) satisfies the Schrodinger equation with the potential on the positive half-space r ≥ 0 with the boundary condition 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 where κ 0 is explicitly Hence, from Eq. (S5.8) 16) where in the last equality we changed variable k = u − i/2. It implies in particular that 17) or equivalently the symmetry (reminiscent of a Nishimori condition, or a Galavotti-Cohen theorem)

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 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 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 where we set V ∞ = lim r→∞ V (r) and enforced the boundary conditions (S5.13). This implies This is still not enough because we require the expansion (S5.19) for r = 1. However, fixing δ > 0 and r, we can write 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 γ, Finally, integrating over r

Small limit
At small one finds that the potential is a harmonic oscillator and 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 where we have defined as in the text 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 and one can check that the Fourier transform ofP dk 2π e ikκP (κ) = Q k (r = 1) (S5.32) as given in (S5.30). This can be seen using the identitŷ

Large limit
At large , under the hypothesis that f (x) and its derivatives decay at infinity, the potential term reaches a constant value 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 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)) √  in agreement with (S5.35) (θ = 1/2 for (S5.38)). The limit in (S5.41) can be easily obtained using (S5.39) using that lim 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 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 H = vˆdx(1 + η(x, t))(T + (x) +T − (x)) . (S6.1) To derive the time evolution of an operator, we introduce the infinitesimal generator of time evolution dĤ = vˆdx(dt + dW t (x))(T + (x) +T − (x)) . Note that because of the Ito's convention we need to keep terms up to the double commutators.