Exact large-scale correlations in integrable systems out of equilibrium

Using the theory of generalized hydrodynamics (GHD), we derive exact Euler-scale dynamical two-point correlation functions of conserved densities and currents in inhomogeneous, non-stationary states of many-body integrable systems with weak space-time variations. This extends previous works to inhomogeneous and non-stationary situations. Using GHD projection operators, we further derive formulae for Euler-scale two-point functions of arbitrary local fields, purely from the data of their homogeneous one-point functions. These are new also in homogeneous generalized Gibbs ensembles. The technique is based on combining a fluctuation-dissipation principle along with the exact solution by characteristics of GHD, and gives a recursive procedure able to generate $n$-point correlation functions. Owing to the universality of GHD, the results are expected to apply to quantum and classical integrable field theory such as the sinh-Gordon model and the Lieb-Liniger model, spin chains such as the XXZ and Hubbard models, and solvable classical gases such as the hard rod gas and soliton gases. In particular, we find Leclair-Mussardo-type infinite form-factor series in integrable quantum field theory, and exact Euler-scale two-point functions of exponential fields in the sinh-Gordon model and of powers of the density field in the Lieb-Liniger model. We also analyze correlations in the partitioning protocol, extract large-time asymptotics, and, in free models, derive all Euler-scale $n$-point functions.


Introduction
The nonequilibrium dynamics of integrable many-body systems has received a large amount of attention recently, especially in view of experimental realizations in cold atomic gases [1][2][3]. It is known that in situations with slow, large-scale variations in space and time, the principles of hydrodynamics hold [4][5][6]. The recently developed generalized hydrodynamics (GHD) [7,8] applies these principles to the presence of infinitely-many conservation laws afforded by integrability. The original works [7,8] strongly suggest that GHD, in the quasi-particle formulation, has wide applicability within quantum systems, including quantum chains and quantum field theory (QFT), requiring only a restricted set of dynamical and kinematical data. These data arise from the thermodynamic Bethe ansatz (TBA) [9][10][11]. In the quantum context (and omitting the simple cases of free particles), GHD has been explicitly worked out in general integrable QFT with diagonal scattering (such as the sinh-Gordon and Lieb-Liniger models) [7,12,13], in the XXZ quantum chain [8,14], and in the Hubbard model, which displays "nested Bethe ansatz" [15], and is expected to apply to all known integrable QFT and Bethe-ansatz integrable models. The structure of GHD, however, transcends its origin from the Bethe ansatz, and GHD can be shown to apply to an even larger variety of models, including classical integrable field theory [16] and classical gases such as the hard rod model [13,[17][18][19] and soliton gases [20][21][22][23][24][25]. The theory has been quite successful, see for instance [26][27][28][29][30][31][32][33]. GHD, as developed until now, is valid at the Euler scale, but viscous and other corrections have been considered, see [19,29,[34][35][36][37][38]. In the present paper, we restrict to the Euler scale. An important problem is that of evaluating dynamical correlations. For definiteness, let an initial state 〈· · ·〉 ini be of the form (1.1) (for any observable O). Here q i (x), i ∈ form a basis of local and quasi-local densities [39] of homogeneous, extensive conserved quantities Q i = dx q i (x) in involution, and β i (x) are parameters, which can be interpreted as generalized local temperatures or local chemical potentials of the integrable hierarchy. We use a continuous space notation x, and the trace notation Tr. This is for convenience, and the problem is posed in its most general setting, for classical (where the trace means a summation over classical configurations) or quantum models, on a one-dimensional infinite space that can be continuous or discrete.
The state (1.1) is an inhomogeneous version of a generalized Gibbs ensemble [40][41][42]. Let the evolution of a local observable O(x) be generated by some homogeneous dynamics that is integrable, for instance with Hamiltonian H, (1.2) Then one would like to evaluate the set of dynamical connected correlation functions 1 〈O 1 (x 1 , t 1 ) · · · O n (x n , t n )〉 c ini (1.3) for local observables O k (x k , t k ). The problem can be divided into two classes. First, if all β i (x) = β i are independent of position, then the initial state is (homogeneous) GGE. The evaluation of exact correlations functions within GGEs is a difficult problem, and in classical models has been little studied. One-point functions of conserved densities q i are directly accessible from the TBA, and those of conserved currents j i (with ∂ t q i + ∂ x j i = 0) were obtained as part of the development of GHD [7,8]. There is also the Leclair-Mussardo formula for GGE one-point functions of generic local fields in integrable QFT [43,44], based on form factors [45][46][47], and formulae for certain one-point functions in the Lieb-Liniger model [48,49] and the sinh-Gordon model [50][51][52]. For GGE two-point functions, various types of spectral expansions exist [53][54][55][56][57], including new results of the Leclair-Mussardo type [58], as well as exact results in free-particle models based on integrable partial differential equations [59][60][61][62][63] (mostly Gibbs states are considered, but the techniques are extendable to GGEs). In integrable quantum spin chains, expressions for correlation functions in Gibbs states [64,65] and in GGEs [66][67][68] have been obtained, but large space-time asymptotics are still to be fully addressed. Stronger results exist in the hydrodynamic regime: Lieb-Liniger particle density correlations from form factors [69], and more generally a set of efficient formulae for two-point functions of all local densities and 1 Here and below, the superscript c means "connected".
Second, more interestingly, let β i (x) depend on the position x in a weak enough fashion. This may arise, in good approximation, as initial ground states or finite-temperature states of quantum or classical systems in weakly varying potentials, or after a (short) local-relaxation time in the partitioning protocol of non-equilibrium steady states [72]. In this case, much less is known. GHD gives direct access to local GGEs describing the mesoscopic fluid cells, hence to all space-time dependent one-point functions of observables whose GGE averages can already be evaluated. However, for two-and higher-point functions, results only exist in the context of free field theory. Importantly, this includes Luttinger Liquids, and gives access, using the local density approximation and related hydrodynamic ideas, to the low-temperature limit of inhomogeneous integrable models, such as the Lieb-Liniger model in inhomogeneous potentials or the Heisenberg chain with in homogeneous interaction coupling, see for instance [73][74][75][76][77][78][79]. Inhomogeneous two-and higher-point functions have never been studied in more general interacting integrable systems.
In this paper, we provide both a first step in the study of correlation functions in inhomogeneous situations, and further develop the theory of correlation functions in (homogeneous) GGEs. We evaluate Euler-scaled dynamical connected correlation functions in inhomogeneous, non-stationary states, in the generality of GHD (without inhomogeneous force fields). The results apply not only to conserved densities and currents, but also to more general local fields where correlation function formulae are obtained purely from the knowledge of GGE one-point functions. The latter are new also when specialized to GGEs.
More precisely, the objects we study are as follows. Consider the scaled initial state , (1.4) for smooth functions β i (x). The scaling in λ guarantees that the Lagrange parameters of the initial state depend weakly on the position. Let us denote by N λ (x, t) a mesoscopic fluid cell: this can be taken as a space-time region whose extent scales as λ ν for some ν 0 < ν < 1, around the scaled point λx, say N λ (x, t) = {( y, s) : ( y − λx) 2 + (s − λt) 2 < λ ν }. The value of ν 0 depends on the subleading corrections to Euler hydrodynamics; if they are diffusive, then we would expect ν 0 = 1/2. Let us also denote by |N λ | = N λ (x,t) d yds its volume. The "Eulerian scaling limit" for correlation functions is defined as the limit for fixed x k 's and t k 's. Here the superscript c means that we take connected correlation functions, and n 0 represents the initial GHD occupation function, which characterizes the initial state at the Euler scale (the GGEs of the initial fluid cells). Fluid-cell averaging, d y k ds k |N λ | · · · , is necessary in order to avoid non-Eulerian oscillations, and averaging can be performed in various ways (see [16] for a discussion of fluid-cell averaging and oscillations). For one-point functions, numerical observations and exact calculations in free models suggest that fluid-cell averaging is not necessary, and one has 〈O(x, t)〉 Eul [n 0 ] = lim λ→∞ 〈O(λx, λt)〉 ini,λ . We propose a generating function method in order to evaluate (1.5), based on combining an Euler-scale fluctuation-dissipation principle with the "nonlinear method of characteristics" introduced in [33]. We expect the generating function method to be valid whenever equal-time correlations vanish fast enough in space. It is expected to work in all quantum and classical systems that have been shown to be accessible by GHD, and applies to conserved densities q i and currents j i . In the cases of two-point functions, we show that the method provides explicit nonlinear integral equations which can in principle be solved numerically, and from which various special cases can be extracted. The results on two-point functions agree with the GHD projection operators derived in [13], and in homogeneous states, reproduce the formulae found in [13,15].
Further, using hydrodynamic projections, we find formulae for Euler-scale two-point functions of arbitrary local fields, expressed purely in terms of their homogeneous GGE averages. To every local field we associate a hydrodynamic spectral function obtained from its GGE averages, which enters the two-point function formula. Combining with the Leclair-Mussardo expansion in integrable QFT (or its counterpart in classical field theory [80]), we obtain form factor series for Euler-scale dynamical two-point functions for any local field. Using the Bertini-Piroli-Calabrese simplification of the Negro-Smirnov formula [50][51][52] we also obtain explicit results for two-point functions of exponential fields in the sinh-Gordon model, and using Pozsgay's formula [81], of powers of the density operator in the Lieb-Liniger model. These constitute the first such exact results not only in inhomogeneous, non-stationary states, but also in homogeneous GGEs.
Finally, we obtain all Euler-scale n-point functions in free models, study two-point functions of conserved densities in the partitioning protocol, obtaining a number of new results for its solution by characteristics, and study the large-time asymptotics of two-point functions from arbitrary inhomogeneous initial conditions. The paper is organized as follows. In Section 2, we review the basics of GHD, with emphasis both on the general framework accounting for all known examples, and on aspects which are important for the study of dynamical correlation functions. In Section 3, we present the main results about correlation functions, including the generating function method, the two-point functions of conserved densities and currents, the hydrodynamic projection interpretation, and the extension to generic local observables. In Section 4, we give examples of the main formulae, in the sinh-Gordon and Lieb-Liniger models, and in free-particle models. In Section 5 we provide some discussion and analysis of the results, including a study of two-point functions in the partitioning protocol, and a precise analysis of the large-time asymptotics of two-point functions for a large class of initial states. Finally, we conclude in Section 6. The details of the computations are reported in appendices.

Review of GHD
Making full sense of the state (1.1) is not a trivial matter. If the infinite sum in the exponential truncates, then -at least in classical and quantum chains -there is a well developed mathematical theory [82][83][84]. In the case of homogeneous states, β i (x) = β i , there are many studies that discuss the precise terms that must be included within the infinite series i β i Q i in various situations, and its convergence in terms of averages of local observables, see the review [41]. A mathematically rigorous framework has been given [85] showing that the infinite sum can be interpreted as a decomposition in a basis of the Hilbert space of pseudolocal charges; in particular, the infinite series itself is a pseudolocal conserved charge. Later, it was understood how GGEs connect to the quasi-particle description of TBA [86], and an in-depth analysis of finite-series truncations and convergence of local averages was given [67].
Here we concentrate on the quasi-particle description of GHD as originally developed [7,8]. The generality of GHD has been claimed in various works and the same basic ingredients extracted, see e.g. [13,15,33]. In order to establish the notation, which follows [7], we recall these ingredients. We further provide general notions concerning correlation functions, and we make a full account of situations with non-symmetric differential scattering phase (or TBA kernel), making apparent the invariance under quasi-particle reparametrization. It has been noted that this general framework needs small adjustments in order to deal with spin-carrying quantities in the massive regime of the XXZ Heisenberg chain, see [14]; we will not consider this subtlety here.

GGEs in the quasi-particle formulation
We denote by S the spectral space of the model. The space S can be seen roughly as the space of all quasi-particle characteristics admitted in the thermodynamics of the model; it is the space of excitations emerging after diagonalizing the scattering in the thermodynamic limit. In general, S is decomposed into disconnected components: each component represents a quasi-particle type, and is a continuum representing the allowed momenta for this quasi-particle type. The spectral space, therefore, has the form of a disjoint union S = ∪ a∈A I a , where A is the set of quasi-particle types, and I a are continuous subsets of copies representing the continua of momenta for each particle type. We will parametrise each continuum by a variable θ ∈ I a , which we will refer to as the rapidity 2 . One may write a spectral parameter as θ = (θ , a) with θ ∈ I a and a ∈ A. We will use the notation (2.1) Besides the set S, the model is specified by giving the momentum and energy functions p(θ ) and E(θ ) respectively, and the differential scattering phase (or more generally the TBA kernel occurring after diagonalization of the scattering) ϕ(θ , α), a function of two spectral parameters. The momentum function p(θ ) defines physical space and specifies the parametrisation used. Without loss of generality, by faithfulness of the parametrisation we assume that it satisfies where p (θ ) = dp(θ )/dθ (here and below the prime denotes a rapidity derivative). The energy function, on the other hand, defines physical time, and equals the "one-particle eigenvalue" (or the equivalent in classical systems) of the conserved charge that generates time translations (the Hamiltonian), see for instance [12]. The differential scattering phase, of course, specifies the interaction. All equations below are independent of the momentum parametrisation θ used. This invariance involves certain transformation properties of the objects introduced, which are either scalar fields or vector fields. Under rapidity reparametrisations, the differential scattering phase ϕ(θ , α) transforms as a vector field (i.e. as ∂ /∂ θ ) in θ , and a scalar field in α, that is For instance, the differential scattering phase is defined, in diagonal scattering models, as ϕ(θ , α) = −i dS(θ , α)/dθ where S(θ , α) is the two-body scattering matrix. The momentum and energy functions are scalar fields, while their derivatives, p (θ ) and E (θ ), are vector fields. Also given is a set of one-particle eigenvalues, scalar fields h i (θ ) for i ∈ associated to the conserved charges Q i . The space spanned by these functions is assumed to be in bijection with a dense subspace of the Hilbert space of pseudolocal conserved charges (this Hilbert space is induced by the inner product defined via integrated correlations, see [85] and the Remark in Subsection 2.2).
The important dynamical quantities, which specify the GGE in the TBA quasi-particle formulation, are an occupation function n(θ ), a pseudo-energy ε(θ ), a particle density ρ p (θ ) and a state density ρ s (θ ), which are all related to each other [10,11]. The former two are scalar fields, the latter vector fields. Associated to these is the dressing map h → h dr [n] , which is a functional of n(θ ) and a linear operator on (an appropriate space of) spectral functions h. We define it, in general, differently for its action on vector fields and on scalar fields: it is defined by solving the linear integral equations The dressing operation preserves the transformation property under rapidity reparametrization. For lightness of notation in this paper, omitting the index [n] means dressing with respect to the occupation function denoted n(θ ), that is h dr = h dr [n] . It will be convenient to employ an integral-operator notation. We introduce the scattering operator T , with kernel T (θ , α) = ϕ(θ , α)/(2π), acting on spectral functions h as (2.4) as well as its transposed T T with kernel T T (θ , α) = ϕ(α, θ )/(2π). By a slight abuse of notation, we also sometimes use n for the diagonal operator acting as multiplication by n(θ ). In these terms, Both the occupation function and the particle density may be taken as characterising a thermodynamic state (a GGE). Other state quantities are related to them: where p (θ ) is a vector field 3 . The relation between the pseudo-energy ε(θ ) and the occupation function n(θ ) depends on the type of excitation mode considered: it is different for quantum fermionic or bosonic degrees of freedom (as discussed in [10]), for classical particlelike modes such as solitons (as discussed in [87,88]) or hard rods (as discussed in [13,29]), and for classical radiative modes occurring for instance in classical field theory (the GHD of classical field theory is developed in [16] based on [87,88]). We have n(θ , a) = ∂ F a (ε)/∂ ε | ε=ε(θ ,a) where the free energy function F a is given by (a is a classical particle) 1/ε(θ ) (a is a radiative mode) (2.7) (recall that the mode type is encoded within the particle type a of the spectral parameter θ = (θ , a)). Note that the free energy function determines the "generalized free energy" of the GGE, given by dθ p (θ ) F a (ε(θ )).
Averages in GGEs will be denoted by 〈O〉 [n] , functionals of the state variable n(θ ). Averages of conserved densities and currents are found to be [7,8] The effective velocity is [7,8,90] v eff (θ ) = (E ) dr (θ ) (p ) dr (θ ) . (2.10) Here we recall that h i (θ ) are scalar fields and E (θ ) and p (θ ) are vector fields. The Lagrange parameters {β i } of a GGE fix the state, formally, via the trace expression One can recover the occupation function n(θ ) from the set {β i : i ∈ }, and vice versa, via a set of nonlinear integral equations: one defines the GGE driving term w(θ ) = i β i h i (θ ), which involves the one-particle eigenvalues h i (θ ) associated to the conserved charges Q i , and one solves ε(θ ) = w(θ ) + (dγ/2π) ϕ(γ, θ )F b (ε(γ)) (where γ = (γ, b)). For our purposes, we mainly need the derivative of n(θ ) with respect to β i . Again the result depends on the type of excitation mode considered, and may be written as where the statistical factor of the mode is (radiative modes). (2.13) The quantities ε(θ ), ρ s (θ ), ρ p (θ ) v eff (θ ) and f (θ ) are all functionals of an occupation function; below we use these symbols for the quantities associated to the occupation function denoted n(θ ).

Generalized fluids in space-time
Recall that the Eulerian scaling limit (1.5) is defined as a large-scale limit, with fluid cell averaging, of connected correlation functions. This exactly extracts the information about the correlations that is present in the physics of Euler fluids. In order to describe it, we need to construct fluid configurations where at every Euler-scale space-time position (x, t) ∈ × lies a GGE. We thus need a family of state functions, which we denote equivalently as n x,t (θ ) ≡ n t (x; θ ), with θ ∈ S the spectral parameter. The function n x,t (θ ), as a function of θ for x, t fixed, is the occupation function describing the GGE in the fluid cell at (x, t). Below we will use the index [n x,t ] for averages in the GGE at the space-time point (x, t), which are functionals of this function of θ . On the other hand, n t (x; θ ) seen as a function of the doublet (x, θ ) for t fixed, is the fluid state on the time slice t. We will use the index [n t ] for functionals that depend on this function of (x; θ ). For instance, the Eulerian scaling limit (1.5) is a functional of the initial state n 0 , while by definition, evolving for a (Euler-scale) time t gives (2.14) Recall that the dressing operation (2.3) as well as the various TBA quantities are all functionals of an occupation function.
For readability, we will use the notation h dr (x, t; θ ) = h dr [n x,t ] (θ ), as well as ρ s (x, t; θ ), ρ p (x, t; θ ), v eff (x, t; θ ) and f (x, t; θ ) for the quantities associated to the occupation function n x,t (α) (as a function of α for (x, t) fixed).
The fluid state on any time slice t takes a factorized form, where on each fluid cell lies a GGE. That is, at large scales correlation functions factorize as is the average of the local (Schrödinger-picture) operator O k (x k ), in the GGE n x k ,t which lies at Euler-scale space-time position (x k , t). In order for the results below to be valid, we in fact require that equal-time, space-separated connected correlation functions vanish fast enough 4 , In particular, according to (2.12), it satisfies the functional derivative equation In accordance with the factorized form (2.15) and especially (2.16), equal-time scaled connected correlation functions have support only at coinciding points. In fact, taking the Eulerian scaling limit, they can be written in the form By integration, one can identify the pre-factor as the full integral of the connected correlation function in the homogeneous local state at x 1 , . (2.20) Note that the scaling factor λ N −1 exactly cancels that coming from the re-scaling of the integration variables, and that thanks to the space integration, it is not necessary anymore to average over fluid cells.
Remark. For every GGE n(θ ), there is a Hilbert space formed by the completion, under the natural topology, of the space of local observables with the (x, t)-dependent inner "hydrodynamic inner product" There is a sub-Hilbert space formed by the set of conserved densities O 1 , O 2 ∈ {q i : i ∈ } within this Hilbert space. The space spanned by h i (θ ), i ∈ is required to be dense within this sub-Hilbert space, and this, for all n = n t (x). This, generically, imposes the inclusion of quasi-local conserved densities. See the review [39] for quasi-local densities, and [85] for a rigorous description of these Hilbert spaces and the way they are involved in generalized thermalization.

Time evolution
Consider a generalized fluid in space-time that is obtained, after the Eulerian scaling limit, by evolving an initial state (1.1) using a homogeneous dynamics as in (1.2), (1.3). This satisfies an Eulerian fluid equation [7,8]. This is the main equation of GHD, which can be written as the convective evolution equation Its "solution by characteristics" was discovered in [33]. Given the initial condition n 0 (x; θ ), one introduces the characteristics, a function u(x, t; θ ), which one evaluates along with the evolved state n t (x; θ ) by solving the following set of equations: In these equations, x 0 is an "asymptotically stationary point": it must be chosen far enough on the left in such a way that n s (x; θ ) = n 0 (x; θ ) for all x < x 0 and s ∈ [0, t] (typically, one should think of it as x 0 = −∞). This provides the evolution from the initial condition n 0 for a time t.
It is worth noting that the function u(x, t; θ ) has the simple interpretation as the position, at time 0, from where a quasi-particle trajectory of spectral parameter θ would reach the position x at time t. Indeed, it solves Thus, defining the trajectory x(t), starting at x(0) = y, via we find Below we assume the following: (i) the state density ρ s (θ ) is positive for all θ , and (ii) the equations (2.23) have a unique solution. Thanks to these assumptions, differentiating with respect to x the second equation in (2.23), we have the inequality (2.27) which imply that the function u(x, t; θ ) is invertible with respect to the position.
Remark. Note that if we assume that the effective velocity v eff (θ ) is a monotonically increasing function of the rapidity θ , then (Appendix A) u (x, t; θ ) < 0 (2. 28) so that that the function u(x, t; θ ) is invertible with respect to the rapidity. The latter condition is satisfied for instance in Galilean of relativistic field theories. This condition slightly simplifies some of the considerations, and in particular it guarantees that (v eff ) (θ ) = 0. In fact, if the latter inequality is not satisfied, then some of the asymptotic results below do not apply. Yet, we will not make use of the monotonicity assumption, but we will implicitly assume that (v eff ) (θ ) = 0 when it appears in denominators, keeping the discussion of how a vanishing derivative of the effective velocity may change some results for the conclusion.

Correlation functions
Despite the factorization properties (2.15) and (2.16) on equal-time slices, scaled connected correlation functions (1.5) are nontrivial when fields do not all lie on the same time slice. That is, a connected dynamical N -point function vanishes, at the Euler scale, as λ 1−N with a generically nonzero coefficient, which is extracted (after fluid-cell average) by taking the Eulerian scaling limit (1.5).
In this section, we develop a recursive procedure that generates all scaled dynamical correlation functions (1.5). The procedure is based on linear responses and an extension of the fluctuation-dissipation theorem to Euler scale correlations. We identify the propagator, propagating from time 0 to time t, as (simply related to) the linear response of n t to variations of the initial condition n 0 . We explain how, in the cases of two-point functions involving conserved densities q i (x, t) and currents j i (x, t), one can obtain from this procedure explicit integral equations. We also explain how one can extend these formulae, combining hydrodynamic projection principles with the Leclair-Mussardo formula, to two-point functions involving other local fields. We finally state the general results for scaled n-point functions in free models.
It is worth noting that in general, correlation functions depend on much more than the information present in the Euler hydrodynamics. For instance, although the knowledge of the GGE equations of states is sufficient to determine the full thermodynamics and Euler hydrodynamics, it cannot be sufficient to determine correlation functions of the type 〈q i (x, t)q j (0)〉 c ini . Indeed, GGE equations of state give information about conserved charges Q i = dx q i (x), but conserved densities q i (x) are defined from these only up to total spatial derivatives of local fields. Thus any result from GHD for two-point correlation function 〈q i (x)q j (0)〉 c ini cannot depend on the precise definition of q i (x). The Eulerian scaling limit (1.5) only probes large wavelengths, and derivative corrections to q i (x) are expected to give vanishing contributions. This is why it is possible to obtain exact results purely from GHD for this scaling limit. Any correction to the Eulerian scaling limit necessitates additional information, hence cannot lie entirely within the present GHD framework.
Euler-scaled dynamical correlations can be seen as being produced by "waves" of conserved quantities ballistically propagating in the fluid between the fields involved in the correlation function. The problem can thus be seen as that of propagating Euler-scale waves from the initial delta-function correlation (2.19), essentially using the evolution equation (2.22). This form of the problem is made more explicit in the case of two-point functions in Subsection 3.3 using hydrodynamic projection theory.

Generating higher-point correlation functions
The main idea of the method is to use responses to local (in the Euler sense) disturbance in order to generate dynamical correlations. Indeed, consider the state (1.1). The response to a small change of the local potential β i (x) at the point x should provide information about the correlation between the observable O (which can be a product of local observables) and the local conserved density q i (x). At the Euler scale (1.5), the functional differentiation with respect to β i (x) brings down the density q i (x), and does nothing else. This is clear in classical models as it follows from differentiation of the exponential function. In quantum models, terms coming from nontrivial commutators between local conserved densities are negligible at the Euler scale: they only give rise to derivatives of local operators, see [12, eqs. 91-93], which can be neglected in Eulerian correlation functions 5 . Therefore, We see that Eulerian dynamical correlation functions are related to response functions. This constitutes a generalisation, both out of equilibrium and to the presence of the higher conserved charges of integrable models, of the fluctuation-dissipation theorem. Let us consider Euler scale correlation functions (1.5) involving charge densities and currents. The one-point functions are given by (2.8), (2.9). Evolving in time and taking the Eulerian scaling limit is simple, Higher-point functions with many insertions of conserved densities are obtained recursively as follows. Let be a product of local observables at various space-time positions. It is convenient to assume that t N = 0, without loss of generality as we can always evolve in time using (2.23) is known as a functional of n 0 (x; θ ). This is the case for N = 1 with O 1 being a conserved density or current (see below for other one-point functions). From this, we may obtain correlation functions 〈 This is of the same form as the correlation at order N : it contains N + 1 observables, where all N previous local observables have been evolved for a time t, and a new conserved density has been inserted at time t N +1 = 0. We obtain: We have used (3.1), (2.14) and (2.18). In this expression, δn t (z; θ )/δn 0 ( y; α) is the functional derivative of the time-evolved occupation function n t (z; θ ) with respect to variations of the initial condition n 0 ( y; α) from which it is evolved. Density-density two-point functions take a particularly simple form thanks to the general formula obtained in [13], where µ is any parameter on which a GGE state n(θ ) may depend, and g(θ ), h(θ ) are any spectral functions (either g is a vector field and h is a scalar field, or vice versa). The functional derivative on the right-hand side in (3.4) may be evaluated using this along with (3.2) (specialized to t = 0), giving Recall that ρ s (x, t; θ ) is the state density (2.6) evaluated with respect to the occupation function at space-time position (x, t). The density-current two-point function can be obtained similarly. Higher-point functions are obtained using (3.4) by further functional differentiation, using similar techniques. The crucial objects in these formulae are the functional derivatives of the time-evolved occupation function n t (x; θ ) with respect to its initial condition n 0 ( y; α). These describe the dynamical responses of the fluid at time t to a change of initial condition. The two-point function only involves the first derivative, while higher-point functions will involve higher derivatives.
Below it will be convenient to define the propagator as a simple conjugation of the first derivative of the evolution operator: In terms of the propagator, the density-density two-point function takes the form Note that the propagator is a vector field as a function of its first argument, and a scalar field as a function of its second. In the following, we concentrate on two-point functions: we explain how to evaluate the propagator via integral equations, and how to go beyond correlation functions involving conserved densities. It turns out that the propagator, as defined in (3.7), satisfies a linear integral equation whose source term and kernel stay well defined even at points where the occupation function vanish. We leave for future studies the developments of expressions for higher-point functions and the evaluation of higher-derivatives of the time evolved occupation function.

Exact two-point functions of densities and currents
The derivation of the following formulae, based on the techniques introduced above, is presented in Appendix B.1. Here we describe the main results.
In order to express the results, it is convenient to introduce the "star-dressing" operation, which for a vector field g(θ ) and a GGE occupation function n(θ ) is defined by Note that without interaction, we have g * dr = 0. We will also need the effective acceleration a eff . This is a functional of n 0 (x; θ ) (seen as a function of (x, θ )). It is defined as a eff is a scalar field, the TBA driving term of the GGE n 0 (x; θ ) (see (2.17)). For our purpose, we may write it in the equivalent forms The effective acceleration encodes the inhomogeneity of the fluid state n 0 (x; θ ).
It will be convenient to see the propagator Γ ( y,0)→(x,t) (θ , α) as the kernel of a linear integral operator acting on scalar fields via contraction on the spectral parameter α: . (3.11) This can be interpreted as bringing the spectral function g from the point ( y, 0) to the point (x, t) starting in the initial state n 0 . We show in Appendix B.1 that the propagator satisfies, for x, y > x 0 (recall (2.23) for the quantity x 0 ), the following integral equation: α). In this and other equations below, functions such as ρ s (x, t; θ ) and f (x, t; θ ) with omitted spectral argument θ , are to be seen as diagonal integral operators, acting simply by multiplication by the associated quantity.
Remark that if the initial state is homogeneous, in which case the evolution is trivial n t (x; θ ) = n(θ ), then we have a eff [n 0 ] (u; θ ) = 0 and u = x − v eff (θ )t, and we find (homogeneous states). (3.13) Here and below, δ S (θ − α) = δ(θ − α)δ b,a for α = (α, a) and θ = (θ , b). In the absence of interaction, is the group velocity, and Further, at t = 0, one obtains The propagator (3.12) allows one to evaluate two-point functions of conserved densities in inhomogeneous states as per (3.8). For two-point functions involving currents, the results are simple modifications of the above, where the effective velocity multiplies the dressed oneparticle eigenvalues. The results are These expressions are similar to those obtained in [13], except for the nontrivial propagator In the homogeneous case, using (3.13), we indeed recover the result of [13]. It is natural to separate the propagator into two terms, (3.20) We will refer to the first term as the direct propagator, and the last as the indirect propagator.
As explained in Appendix B.2, the indirect propagator satisfies the following linear integral equation: Here the dressed scattering operator is and T dr (z, t; θ , γ) is, as a function of θ , γ, the kernel of T dr (z, t) (with dressing with respect to the state [n z,t ]). T dr (z, t; θ , γ) is a vector field as function of the θ , and a scalar field as function of γ. The root set is If the effective velocity is monotonic with respect to the rapidity, by virtue of (2.28), the function u(x, t; θ ) is locally invertible on the rapidity θ , wherefore the set θ (x, t; y) contains at most one element θ = (θ a , a) per particle type a ∈ A. In general, however, the set may contain more solution per particle type. In terms of (3.20), we have for instance in the first term on the right-hand side).

Connection with hydrodynamic projections
The main ideas of hydrodynamic projections, in the cases of two-point functions at the Euler scale, can be gathered within two statements. First, correlations are transported solely by (ballistically propagating) conserved densities. Second, the overlap between a local observable and such a propagating conserved density, in the fluid cell containing the local observable, is obtained by the hydrodynamic inner product (2.21) within this cell. Using this, the expressions (3.17)- (3.19), involving currents, are in fact consequences of (3.16) using hydrodynamic projection theory. Here we first re-write the expressions obtained above in the hydrodynamic-projection form. We then show how taking this form implies the Euler-scale fluctuation-dissipation principle we have used to derive the expressions (3.17)-(3.19).

Re-writing in hydrodynamic-projection form
First, the integral operator S ( y,0)→(x,t) , acting on scalar fields and giving vector fields, that generates the charge two-point function as is given by Relation (3.27) is simply a re-writing of (3.16). Note that by symmetry of the correlation func- where T denotes transpose. From hydrodynamic projection, it is known that correlation functions involving currents are obtained by using the linearized Euler operator, which, in a GGE state n(θ ), is given by [13] it acts on vector fields and gives vector fields. Re-writing (3.17)- (3.19), currents correlations are obtained as: (for both observables being currents). (3.29) These are indeed expressions that are expected form hydrodynamic projection principles [13].
In the homogeneous case we have that n t = n 0 and, using (3.13), that 0) . In this case one usually denotes the operator as S(x − y, t), and we recover S(x − y, t) T = S(x − y, t).
Further, according to hydrodynamic projection principles, one would expect S ( y,0)→(x,t) to solve the evolution equation , for a GGE state n(θ ), is the correlation operator. Its matrix elements, in the space of conserved densities, are the connected integrated two-point functions C (2.20) and (2.21)), and as an operator it is [13] Eq. (3.30) is the generalisation to space-time dependent states of the equation that was solved in [13] in order to obtain Euler-scale correlations in homogeneous states. It is an explicit expression of the problem of propagating Euler-scale waves from the initial delta-function correlation (2.19), in the case of two-point correlations. It is simple to verify that indeed (3.30) follows from the results obtained: the initial condition holds by using (3.15), and (3.30) follows from (3.16) and (3.17) and the conservation law for local densities.

From hydrodynamic projections to Euler-scale fluctuation-dissipation principle
Above, we saw the hydrodynamic projection evolution problem emerging as a consequence of defining space-time correlation functions using an Euler-scale fluctuation-dissipation principle. Let us now reverse the logic: let us take (3.26) with (3.30), along with the appropriate correct initial condition as stated after (3.30), as a definition of the scaled dynamical twopoint functions of conserved densities. From this, let us show that the Euler-scale fluctuationdissipation principle (3.1) holds for two-point functions of conserved densities. The relation = 0 follows from the basic GHD results. Taking functional derivatives with respect to β( y), it implies The result for these functional derivatives are the expressions on the right-hand sides of (3. 16) and (3.17). Using these in the above equation, we indeed find that (3.27) satisfies (3.30) with the correct initial condition. Hence, we conclude that

Two-point correlations of generic local observables
Let n(θ ) be a GGE. We consider the inner product 〈O|q i 〉 [n] , and it is assumed that spatial correlations within the state n(θ ) decay faster than the inverse distance. Without loss of generality we assume O to be hermitian. If the average of O is known in a generic GGE n(θ ), then we can use In integral operator form, by linearity of the result in h i , there exists a scalar field V O (θ ) (a hydrodynamic spectral function associated to the local field O), which is also a functional of n(θ ), such that Here for later convenience we introduced the factors ρ p (θ ) f (θ ) and used h dr i (θ ) instead of h i (θ ). For instance, according to results of [13], we have the overlap, within a fluid cell of GGE n(θ ), between conserved densities, as per (2.21) -this is the correlation matrix of conserved densities, the case N = 2 of (2.20). In integral operator form, this is the correlation operator (3.31) introduced above. Let us now consider a generalized fluid, with space-time state described by n x,t (θ ) ≡ n t (x; θ ). According to hydrodynamic projection principles, Euler-scale correlation functions can be written as The first equality is explained as follows. Reading from the right to the left, we first overlap the observable O with a complete set of conserved quantities q l , with respect to the inner product (2.21) for the state n y,0 (at the space-time point ( y, 0), where the observable O lies). Because the conserved quantities q l 's don't necessarily form an orthonormal set, we introduced the inverse correlation matrix C [n y,0 ] at the space-time point ( y, 0). These two factors represent the amplitude for O to produce Euler propagating waves of conserved quantities. We then "transport" these waves from ( y, 0) to (x, t) by using the dynamical two-point function 〈q j (x, t)q k ( y, 0)〉 [n 0 ] between conserved densities. Finally, we represent the amplitude for the transported wave to correlate with O by overlapping with O with respect to the inner product at x, t, introducing the inverse correlation matrix C [n x,t ] for orthonormality. The second equality is a re-writing in terms of integral operators, using (3.33) and (3.26). Finally, the last equality is obtained by replacing with the expressions (3.27) and (3.31). As a check, note that using (3.34), the above indeed reproduces the formulae (3.16)- (3.19). In particular, in homogeneous states, we use (3.13) and find where ξ = x/t and θ (ξ) is the set of solutions to v eff (θ ) = ξ. It turns out that, in integrable QFT, there exists a formula for the averages 〈O〉 [n] in GGEs, for any local field O [43,44]. This formula, called the Leclair-Mussardo formula, involves an infinite summation over multiple integrals of form factors of the field O. Nevertheless, its truncations can be used to numerically approximate expectation values. The Leclair-Mussardo formula was proven in [81], and it was used in [7] in order to provide further evidence for the proposed one-point averages of currents.
The formula has the structure of a sum over all numbers of particles k of "connected" diagonal matrix elements F O k (θ 1 , . . . , θ k ) = 〈θ 1 , . . . θ k |O|θ 1 , . . . , θ k 〉 conn. of the field O (this is defined in [43]). Consider a GGE n(θ ). Then the formula is Recall that in our simplified notation, θ represents the combination of a rapidity and any particle type the model may admit. Importantly, in this formula, the information about the state is fully contained within the integration measure dθ n(θ ). Using (3.32) and (2.12), we therefore find The function F O k is symmetric in all its arguments, and we may identify where ρ s (θ ) is the density of state, given in (2.6). The state dependence is within the integration measure and the density of state; the regularized diagonal matrix element is purely a property of the field O. It is interesting to re-specialize to O being a conserved density or current in order to verify that one indeed recovers (3.34) from (3.39). This is done in Appendix C.
Using (3.27) and reverting to the explicit notation Γ ( y,0)→(x,t) (θ , α) for the propagator via It is remarkable that such a complete formula exists in integrable field theory, for very general dynamical Euler-scale two-point correlation functions of local fields in inhomogeneous, non-stationary states. This formula is new both in the inhomogeneous case, and in the case of a homogeneous GGE; in the latter case, recall that the propagator simplifies to (3.13).
Remark. It is very likely that the form factor series (3.40), in homogeneous GGEs, can be obtained directly using an appropriate spectral expansion of the two-point function. Indeed, the structure of this series is extremely suggestive of the techniques introduced in [53,56], based on the ideas of the Gelfand-Naimark-Segal construction. From these techniques, the trace expression representing the GGE average of a product of local fields is expressed as an expansion in "GGE form factors" very similar to the form factor expansion of vacuum twopoint functions, in which each GGE form factor is itself a GGE trace of a single local field with additional particle creation / annihilation operator inserted. The leading term at the Euler scale is that with one particle and its hole at the same rapidity, so that, pictorially, Each single-field trace can be evaluated using the Leclair-Mussardo formula, giving a righthand side similar to that of (3.40). We hope to come back to this problem in a future work.

Sinh-Grodon model
The sinh-Gordon model is an integrable relativistic QFT with Lagrangian density for a real scalar field Φ, where g is a coupling parameter and M is a mass scale. Its TBA description contains a single particle of Fermionic type, so that S = , θ = θ and f (θ ) = 1 − n(θ ).
We may choose θ as the rapidity, with p(θ ) = m sinh θ and E(θ ) = m cosh θ , and the physical mass and differential scattering phase are given by This includes the density of momentum (h 1 − h −1 ) and the density of energy (h 1 + h −1 ), as well as higher-spin local conserved densities. The formulae derived in subsections 3.1 and 3.2 immediately give correlation functions for these densities in inhomogeneous, non-stationary situations (generalizing the homogeneous, stationary two-point function formulae found in [13]). One may also obtain two-point correlation function formulae for other local fields that are not local conserved densities and currents, using the results of subsection 3.4. There exist explicit results for one-point function of certain exponential fields in GGEs, avoiding the complicated LM series, which thus can be used to extract V O (θ ) as defined in (3.32), (3.33). It was found in [50][51][52] that 〈e k gΦ 〉 [n] = 1 + 4 sin(πa(2k + 1)) dθ 2π where is (in our interpretation) the k-dressing of the function e −1 (θ ) = e −θ seen as a vector field: the dressing with respect to a different, k-dependent scattering kernel given by . (4.6) Defining is the k-dressing of the function e(θ ) = e θ seen as a scalar field. Using (3.33), we then find, for all k ∈ , By the 2 symmetry 6 Φ → −Φ we have V −k (θ ) = V k (θ ). Further [52], there is a symmetry 〈e k gΦ 〉 [n] = 〈e (k+a −1 )gΦ 〉 [n] , and thus V k (θ ) = V k+a −1 (θ ), which, for irrational couplings a, allows us to reach arbitrary values of k ∈ . The resulting V k (θ ) can be inserted into (3.35) and (3.36) in order to get Euler-scale two-point correlation functions of fields e kgΦ and e k gΦ for any k, k . The generalized hydrodynamics of classical limit of the sinh-Gordon model was investigated in [16], where the classical limit of V k (θ ) was derived. Euler-scale two-point functions obtained by (3.36) were verified to agree with direct numerical simulations.

Lieb-Liniger model
The repulsive Lieb-Liniger model is defined (for mass equal to 1/2) by the second-quantized Hamiltonian for a single complex bosonic field Ψ, where c > 0 is a coupling parameter. It is Galilean invariant, and its TBA description contains a single quasi-particle type, so we take S = and write θ = θ . One may choose the parametrization given by the momentum, θ = p ∈ (so that p (θ ) = 1 and ρ s = 1 dr /(2π)). There are various TBA descriptions possible, but in one convenient description, the quasi-particle is of Fermionic type, hence f (p) = 1 − n(p). In this description, the differential scattering phase is given by (4.12) Again, as a set of natural local fields, one may consider the local conserved densities and currents of the model; they correspond to the spectral functions h r (p) = p r−1 , r = 1, 2, 3, . . . (4.13) This includes the density of particles (r = 1), the density of momentum (r = 2) and the density of energy (r = 3). Again, the formulae derived in subsections 3.1 and 3.2 give correlation functions for these densities in inhomogeneous, non-stationary situations.
One may also obtain two-point correlation formulae for other local fields that are not local conserved densities and currents, using the results of subsection 3.4. Consider the K th power of the particle density, (4.14) It was shown in [49] that in a homogeneous state characterized by the occupation function n(p), its average takes the form Taking the β i -derivative is simple, by using (2.12) and the general formula (3.5): where is defined as a function of p s with p r =s fixed parameters. From this we identify where 1 dr (p) is the dressed constant function 1. This gives two-point functions by insertion in (3.35) and, in the homogeneous case, in (3.36). We note finally that in the very recent paper [91] new expressions for expectation values of the fields O K are obtained using the non-relativistic limit of the sinh-Gordon model and the results of [50-52] recalled above. These appear to be more efficient. By using the methods shown here, this can in turn be used to obtain different expressions for V O K (θ ). This is worked out in [92].

Free particle models
In free particle models, formulae (3.16) -(3.19) simplify. Using (3.14), the fact that the dressing operator is trivial, and n t (x; θ ) = n 0 (x − v gr (θ )t, θ ), one obtains the simple expression Here (v gr ) −1 (ξ) = {θ : v gr (θ ) = ξ}. The integral form has the clear physical interpretation of a correlation coming from the ballistically propagating particles on the ray connecting the two fields. It is evaluated by summing over all solutions to v gr (θ ) = (x − y)/t, of which there is at most one for every particle type. One similarly obtains current correlations by multiplying by factors of v gr (θ ). For instance, in the quantum Ising model (a free Majorana fermion), where there is a single particle type, one has p(θ ) = m sinh θ , E(θ ) = m cosh θ , v(θ ) = tanh θ and f ( y, 0; θ ) = 1 − 2πρ p ( y, 0; θ )/(m cosh θ ). If the initial state is locally thermal with local inverse temperature β( y), then . (4.20) In this case, the energy density dynamical two-point function (writing q 1 = T 00 , the time-time component of the stress-energy tensor) is zero outside the lightcone, and otherwise is where s = t 2 − (x − y) 2 is the relativistic time-like distance between the fields. Similarly, consider the correlation function of particle densities (writing the density as q 0 = n) in the free nonrelativistic, spinless fermion, evolved from a state with space-dependent temperature β( y) and chemical potential µ( y). For instance, this describes the Tonks-Girardeau limit of the Lieb-Liniger model. We find 〈n(x, t)n( y, 0)〉 Eul In free particle models, it is possible to develop the full program outlined in Subsection 3.1, and to obtain explicit expressions for every N -point correlation functions, at least for conserved densities. The procedure is quite straightforward, and the results can be expressed as follows. Let w(x; θ ) = ∞ j=0 β j (x)h j (θ ) be the TBA driving term of the GGEs (2.17). Then, for N = 2, 3, 4, . . . we have Here the functions g N ( y; θ ) are defined via generating functions as where a is the quasi-particle type associated to θ = (θ , a), and where the free energy function F a (w) is given in (2.7). In (4.23), the delta-functions on the right-hand side constrain the , which equal v gr (θ ). Therefore, we may replace the argument ( or by x k − v gr (θ )t k for any j = k.
In fact, all Euler-scale correlation functions, for N = 1, 2, 3, 4, . . ., can be obtained formally by using generating functionals over the generating parameters k (x) via where on the right-hand side θ = (θ , a).

Conjecturally, correlation functions involving currents are obtained by replacing factors
We note that (4.19) and (4.23) give general, explicit expressions for Euler-scale N -point correlation functions of conserved densities in free theories. There is no integral over spectral parameters: for every particle type, there is a single velocity that contributes to the connected Euler-scale correlation, which is the velocity of the particle propagating from the initial to the final point. For similar reasons, correlation functions for N ≥ 3 have a delta-function structure which imposes colinearity of all space-time positions. Connected Euler-scale correlations can only arise from single quasi-particles travelling through each of the space-time points. Due to this, all correlation functions depend on the initial state only through the local state at a single position: the position, at time 0, crossed by the single ray passing through all space-time points (this is y in (4.19) and more generally x k − v gr (θ )t k in (4.23)). Therefore, the only effect of the weak inhomogeneity is to give a dependence on the state via this single position. All these properties are expected to be broken in inhomogeneous states of interacting models. The dependence is not solely on the state at a single position, as the knowledge of the state at other positions is necessary in order to evaluate the effect of a disturbance on the quasiparticle trajectories (hence to evaluate the response function). Similarly, we do not expect a delta-function structure for higher-point functions in interacting models.

Formulae (3.16)-(3.19) can be given a relatively clear interpretation.
A correlation function is expressed as an integral over all spectral parameters of the product of the quantity of charge of the first observable, h dr i (x, t, θ ), carried by the spectral parameter θ and dressed with respect to the the local bath n t (x), times the propagation to the point (x, t), of the quantity of charge of the second observable h dr j (0, y; α) dressed by the local bath n 0 ( y). The propagation fac- representing the effect of quasi-particles ballistically propagating from ( y, 0) to (x, t), as well as the density ρ p (x, t; θ ), which weigh this effect with the quantity of quasi-particles actually propagating. It also includes the factor f (x, t; θ ), which modulates the weight according the quasi-particle statistics; for instance, for fermions, this factor forbids the entry of a new quasi-particles if the local occupation function is saturated to n t (x; θ ) = 1, making the correlation effect of this quasi-particle vanish. The same structure occurs for more general local observables in (3.40), where the only complication is in the dressing by the local bath, which involves a sum over all form factors. The nontrivial physics of the propagation of correlations -the response of the operator at (x, t) to a disturbance by the observable at ( y, 0) -is fully encoded within the propagator. There is an apparent asymmetry between the initial position ( y, 0) and the final position (x, t), as the factors ρ p (x, t; θ ) f (x, t; θ ) are only present for the latter position. However we note that the points ( y, 0) and (x, t) are not independent, they are related by the evolution equation (2.22): thus every quasi-particle at (x, t) has an antecedent at ( y, 0). In nontrivial (inhomogeneous, interacting) cases, asymmetry is also explicit in the equation defining the propagator (3.12), where quantities pertaining to the initial state density appear naturally. In these cases, the propagator is not an intrinsic, state-independent property of quasi-particle propagation: it is affected by the initial state in nontrivial ways. It is possible to write the two-point functions in more symmetric ways, such as in (3.26), but the choice (3.7) for the propagator has the advantage that (i) it specializes to delta-functions in simple cases (3.13), (3.14), (3.15), and (ii) its defining integral equation (3.12) only involves quantities that are explicitly non-divergent in any GGE.
The propagator is composed of two elements, as per (3.20). The first, the direct propagator, comes from the direct propagation of the disturbance of the initial state due to the observable at ( y, 0). At the Euler scale, this travels with quasi-particles along their characteristics described by the function u(x, t; θ ). Only particles with just the right spectral parameter will travel from ( y, 0) to (x, t), and thus this element should indeed give a delta-function contribution to the propagator.
The second element, the indirect propagator ∆ ( y,0)→(x,t) (θ , α), is more subtle. It comes from the change of trajectories of quasi-particles due to the disturbance at ( y, 0). In the explicit calculation in Appendix B.1, it is seen as the change of the characteristics u(x, t; θ ) upon differentiation with respect to the Lagrange parameter β j (0). The indirect propagator ∆ ( y,0)→(x,t) (θ , α) is still applied to the local dressed quantity h dr j (0, y; α), as it is this quantity that is to travel on the slightly modified trajectory in order to create correlations. However all spectral parameters α may generically participate, instead of a single one, because all are involved in determining the trajectory.
As has been noted above, the indirect propagator vanishes in homogeneous states and in free-particle models (see (3.13) and (3.14)). The above interpretation makes these fact clear: in homogeneous states, it does not matter if the quasi-particle trajectories are modified, as the state is everywhere the same; and in free models, the trajectories do not depend on the local states, thus are not affected by the disturbance at ( y, 0).
One can see that the indirect propagator is largely controlled by the effective acceleration a eff [n 0 ] (u(x, t; θ ); θ ), and in particular that it vanishes if the latter does. Recall that the effective acceleration was initially introduced in order to describe force terms due to external, space-dependent fields (that is, weakly inhomogeneous evolution Hamiltonians) [12]. Here, it instead encodes the (weak) spatial inhomogeneity of the initial state. The space-dependent GGEs in the fluid cells of the initial fluid state are associated with an inhomogeneous "Hamiltonian" i β i (x)Q i , and it would be an evolution with respect to this that would generate force terms controlled by the acceleration field a eff [n 0 ] (z; θ ). Here we see that the effective acceleration instead determines, in part, the way in which characteristics are modified due to disturbances.
It is in principle possible to numerically evaluate the expressions (3.16)- (3.19). Recall that the exact solution (2.23) can be solved very efficiently by iteration, as explained in [33]. Therefore, we may assume n 0 (z; θ ), n t (z; θ ) and u(z, t; θ ) to be readily numerically available for all z, θ . It is then straightforward to evaluate dressed quantities, which can be done by solving (2.3) either by iteration, or by discretizing the linear integral equation and inverting the resulting matrix 1 − T n. Therefore, the only ingredient in (3.16)-(3.19) that is not readily numerically available from previous works is the propagator Γ ( y,0)→(x,t) (θ , α). For this, we write it in the form (3.20). We can then evaluate the indirect propagator by solving (3.21), where the function g(θ ) is chosen as either h dr j ( y, 0; θ ) or v eff ( y, 0; θ )h dr j ( y, 0; θ ) depending on the correlator sought for. The source term can be evaluated from the quantities already numerically available, and (3.21) is a linear integral equation which can be solved, for instance, by iterations. One difficulty might lie in the evaluation of derivatives, for instance u (x, t; γ). One might find it more efficient to differentiate the integral equation (2.23) and solve for u (x, t; γ) instead of directly taking the derivative numerically.

Partitioning protocol (domain wall initial condition)
Consider the evolution from an initial density operator where the state is spatially separated between two different homogeneous states on the left and right. This is referred to as the partitioning or cut-and-glue protocol, or as the evolution with domain wall initial condition, and has been studied extensively, see the review [72]. Even though the initial condition is not smooth, as the initial generalized temperatures display an abrupt jump at the origin, profiles quickly smooth out and the fluid approximation is very accurate after a small relaxation time. The GHD solution [7,8], obtained with initial fluid state of the form gives extremely accurate predictions, as verified in the XXZ model [8] and in the hard rod gas [29]. The solution is a set of ray dependent states where ξ = x/t. Below we will denote scaled correlation functions in this protocol. Let us analyze certain correlation functions in this setup. Naively, one might think that scaled two-point correlation functions (3.16)- (3.19) for two fields lying on the same ray should equal those in the homogeneous state of this ray, as obtained using (3.13): correlations should be carried by particles traveling along this ray alone. This is however incorrect. In order to see this, we consider two different situations. For simplicity we concentrate on charge-charge correlations (3.16), but a similar analysis holds in other cases. See Appendix D for a study of the characteristics in the partitioning protocol.

Correlations on a ray away from connection time
Consider the initial domain-wall state to be at time −t 0 < 0 (so n 0 is the fluid state after the evolution by t 0 from the domain wall), and let (x, t) and ( y, 0) lie on the same ray emanating from (−t 0 , 0), that is ξ = x/(t + t 0 ) = y/t 0 . Then the fluid state is the same at (x, t) and at ( y, 0). Thus we have First let us look at the contribution from the direct propagator δ( y − u(x, t; θ ))δ S (α − θ ) (see (3.20)). For (x, t) and ( y, 0) on the same ray, the only solutions θ to y = u(x, t; θ ) are the solutions to y = x − v eff (θ )t. However, we cannot replace δ( y − u(x, t; θ )) by δ(x − y − v eff (ξ; θ )t), as would be required to reproduce the homogeneous correlators according to (3.13). Indeed, the variation, with respect to θ , of u(x, t; θ ) is not the same as that of v eff (ξ; θ )t, due to the second equation in (D.14). Instead, the direct-propagator contribution to the two-point function is The contribution from the indirect propagator ∆ ( y,0)→(x,t) (α; θ ) gives an additional correction. This contribution is generically nonzero, in particular the state at t = 0 is not homogeneous and thus the effective acceleration a eff [n 0 ] (x; θ ) is nonzero. Therefore, as compared to the homogeneous correlator obtained using (3.13) in (3.16), there are two corrections: the factor 1/V (θ ) in the direct-propagator contribution (5.4), and the indirect-propagator contribution. We have not shown that these two corrections don't cancel each other, but this seems unlikely. Both corrections are due to the fact that the insertion of an observable in a correlation function perturbs the state as seen by other observables, and that due to the nonlinearity of GHD, this perturbation generically affects the trajectories of quasi-particles. Thus other rays are explored, and the two-point function is not that in the homogeneous state of a single ray.

Correlations with one observable at connection time
Second, consider the initial domain wall to be at t = 0. In this case, the state is locally homogeneous at ( y, 0) for any y ∈ \ {0}, therefore a eff [n 0 ] ( y; θ ) = 0. As a consequence only the direct propagator remains, The expression for the scaled two-point function simplifies to a finite sum, as per (3.25). Then we have, for any x, t, y and with ξ = x/t, ; γ), (5.6) whereũ(ξ; θ ) = u(x, t; θ )/t (see Appendix D). Taking y → 0 ± , we can use again (D.14) and the argument above to obtain (5.7) This again looks very similar to the two-point function in a homogeneous state, except for two differences: the factor V (θ ), and the fact that the state at ( y = 0 ± , 0) is not equal to that on the ray ξ that emanates from the origin: it is instead the initial condition, equal to the state at ξ = ±∞. Thus, again, the two-point function on a ray is not that in the homogeneous state of that ray.
The question of the two-point function with y = 0, that is, with one observable within the original discontinuity, is more subtle and answered below.

Long-time asymptotics
Consider an initial state n 0 (x; θ ). Suppose it has well-defined asymptotic behavior at large distances, where it becomes homogeneous: In particular, we suppose that x 0 can be set to −∞ in (2.23). Suppose also that the asymptotic is uniform enough, so that the following integrals converge absolutely: (here and below we denote by ρ ± s (θ ) = lim x→±∞ ρ s (x, 0; θ ) the asymptotic forms of the initial state density). For instance, the initial state could be a state that varies nontrivially only on some finite region. Consider the long-time limit t → ∞ of scaled two-point functions (3.16)- (3.19), along rays x = ξt with y, ξ fixed. By a simple scaling argument, they should decay like 1/t. We provide a derivation of the coefficient of this decay: Again we concentrate on the charge-charge two-point function as the derivation and result is easy to generalize to the currents. In order to derive this result, we further assume that in the limit t → ∞ along any ray x = ξt, one obtains the state n(ξ; θ ), given by (5.2), of the partitioning protocol with initial condition specified by n ± 0 (θ ): lim t→∞ n t (ξt; θ ) = n(ξ; θ ), n(ξ; θ ) from the partitioning protocol with n R,L (θ ) = n ± 0 (θ ). (5.11) We provide in Appendix E a proof under certain more basic assumptions of uniform convergence. Here and below, for lightness of notation, we take the convention that GHD functions explicitly evaluated on a ray, say ξ, instead of a space-time doublet (x, t), are understood as the functions obtained in this limit, for instance ρ s (ξ; θ ) = lim t→∞ ρ s (ξt, t; θ ). These are set by the solution to the partitioning protocol (5.2), see also Appendix D. From this viewpoint, we note that the exact initial condition n 0 (x; θ ) provides a regularization of the initial discontinuity at x = 0 of the partitioning protocol.
An important observation of the result below is that A i j (ξ; y) is not determined solely by the partitioning protocol; in particular it is not the coefficient obtained in either of the two situations studied in Subsection 5.2, and it depends on the point y and on the details of n 0 (x; θ ). What this means is that, from the viewpoint of the partitioning protocol, correlation functions on a single ray, with one space-time point being at the initial time t = 0 and lying on the initial discontinuity of the protocol, explicitly depend on the regularization n 0 (x; θ ) of this initial discontinuity, and on the exact position y, within the regularized region, of the observable at initial time.
Consider the following quantity, which encodes the difference between the regularized initial condition n 0 (x; θ ) and the discontinuous one determined by n R,L (θ ). Given a point y ∈ , a ray ξ and a time t, we look for the spectral parameters θ of quasi-particles starting at y that reaches the position ξt at time t, under the full initial condition n 0 (x; θ ). If the effective velocity is monotonic with respect to the rapidity, then thanks to (2.28), this is unique once the quasi-particle type is determined. In general, we simply consider the set of such θ . We then look for the position r of a quasi-particle θ that would reach the same point (x t, t), but in the partitioning protocol. See Fig. 1. Finally we take the limit t → ∞ of this position. In Figure 1: A pictorial representation of how to evaluate the quantity r, given ξ, t, y. Start at the point y, and find the quasi-particle's rapidity which is such that its trajectory joins y with (ξt, t), in the full problem with initial condition n 0 (x; θ ). Then, using the same quasi-particle type and rapidity, evaluate the backward trajectory from (ξt, t) in the partitioning protocol. The value of r is the position obtained at time 0. In this picture, the shade indicates the space-time region where the fluid states in the full problem and in the partitioning protocol are substantially different, thus affecting the trajectories. this limit, θ ∈ θ (ξ) (that is v eff (ξ, θ ) = ξ). Given this value of θ , the ray ξ is known uniquely (see Appendix D), and thus it fully encodes the ray ξ. The result is r( y; θ ), which depends on both y and on this limiting value of θ . This is defined for all values of θ (there is always a solution to v eff (ξ, θ ) = ξ). In formulae, this is expressed as follows in terms of the functioñ u(ξ; θ ) of the partitioning protocol, which has the explicit form (D.9). We define θ t (whose depence on ξ, y we keep implicit) as u(ξt, t; θ t ) = y, and then r( y; θ ∞ ) = lim t→∞ tũ(ξ; θ t ) with θ ∞ = lim t→∞ θ t ∈ θ (ξ).
The above defines r( y; θ ) in a very delicate way, that involves the full time evolution: one needs to evaluate the finite difference between the end-points of two trajectories that start far in time and stay near to each other for a long time. In order to go further, we need to make certain assumptions about r( y; θ ), which appear to be natural but which we do not know how to verify explicitly. The main assumption is simply that r( y; θ ) is finite. This seems natural if the space-time region where the effects of the regularized partitioning is felt, is of finite extent, as pictorially suggested in Fig. 1. Other more subtle assumptions relate to the exchange of y-derivative and large-time limit, see Appendix E.2.
Under these assumptions, we show in Appendix E.2 that the following integral equation holds: dz ρ s (z, 0; θ ).

(5.12)
Recall the dressed scattering operator (3.23). Equation (5.12) is a powerful result, as it determines r( y; θ ) entirely in terms of initial data, without the need for time evolution. Even more powerful is the fact that, although the left-hand side depends on θ , it is independent of y. Thus, by uniform convergence (5.9), we must have This is simply saying that for y far from the regularization region, there is no difference with the partitioning protocol. Further, differentiating with respect to y, we obtain ∂ r( y; θ ) ∂ y = ρ s ( y, 0; θ ) ρ sgn(r( y;θ )) s (θ ) > 0, (5.14) thus r( y; θ ) is monotonic in y. This means that r( y; θ ) has a unique zero y (θ ), which determines its sign: It is this zero that plays a fundamental role for the long-time asymptotics of correlation functions. Consider the sign function σ( y; θ ) = sgn y − y (θ ) = sgn(r( y; θ )).
An equation determining this zero is inferred from (5.12): The function r( y; θ ) then takes the simple form dz ρ s (z, 0; θ ) for y ≷ y (θ ).

(5.18)
We show in Appendix E.3 that: This provides the long-time asymptotic coefficient explicitly in terms of initial data. In the special case where both sides have the same asymptotics, the partitioning protocol is homogeneous, and the formula simplifies to: where GHD quantities depending only on the spectral variable are to be evaluated in the asymptotic GGE state n(θ ). Here we have used the fact that V (θ ) = 1 in the homogeneous case, as v eff (θ ) = ξ (θ ) (see (D.15)).
As a check, we can verify that the limit | y| → ∞ of (5.21) gives the two-point correlation function in the homogeneous case (3.16) with (3.13) (whose full dependence on time is in the factor t −1 ): (same left and right asymptotics).
(5.22) Indeed, this follows from (using (5.8)): We see that the homogeneous correlation function is at the point y = 0: this is natural, as on the left-hand side, the limit t → ∞, x = ξt is taken before | y| → ∞.
We can similarly check that the limit y → ∞ of (5.19) gives the the two-point correlation function (different left and right asymptotics) (5.24) in the partitioning protocol, see (5.7). We therefore find the natural result that the limit | y| → ∞ of the regularized partitioning protocol, where y starts within the inhomogeneous region that regularizes the discontinuity and goes away from it, exactly agrees with the limit y → 0 ± of the exact partitioning protocol, where y starts within the homogeneous region and goes towards the discontinuity.
Remark. We note a somewhat surprising result that is derived in Appendix E.3, and that leads to the particular form of the results expressed above. It can be expressed equivalently as a "sum rule" 25) or as an "occupation equipartition" relation, The latter relation is extremely nontrivial, as it relates the zero y (θ ) to state properties at the asymptotics and at the point y (θ ) only, while (5.17) defines the zero in terms of states at other points as well. Also, it is not a priori obvious that (5.26) has a solution at all. The derivation we provide in Appendices E.2 and E.3 imply that there is at least one solution. We believe this is deeply related to the requirements of finiteness of r( y; θ ) and the possibility of exchanging y-derivative and large-t limit.

Conclusion
In this paper we have obtained exact expressions for dynamical connected correlation functions at the Euler scale in non-equilibrium integrable models. These represent correlations obtained under unitary time evolution from inhomogeneous density matrices in quantum models, or deterministic evolution from random initial configurations in classical models. The time evolution is taken to be the homogeneous evolution of the integrable model (and thus this excludes the cases of evolutions in external potentials). The results are expressed solely in terms of quantities that are available within the thermodynamic Bethe ansatz framework. They are valid at the Euler scale, where variations of averages of local fields occur on very large scales. Interestingly, this shows that hydrodynamic ideas provide, in principle, all large-scale correlation functions. The range of applicability of our results is the same as that of GHD, and thus includes a wide variety of integrable models. Our derivation is based on a natural Euler-scale fluctuation-dissipation principle combined with the exact GHD general solution found in [33].
We showed that our results agree with the general principles of the hydrodynamic projection theory, with in particular the hydrodynamic operators found in [13].
We also showed how two-point functions of arbitrary observables can be obtained from the knowledge of their one-point functions using the hydrodynamic projection theory. From the Leclair-Mussardo formula, valid for one-point functions, we therefore obtained Euler-scale two-point functions as infinite form factor series. This formula is new both in the inhomogeneous case, and in homogeneous GGEs. We also remark that recently, an exact recursion relation was obtained for expectation values of vertex operators of the form e aφ in the sinh-Gordon model [50][51][52]. This can be used to extract some of the spectral function V e aφ , and deduce their Euler-scale two-point functions using (3.35).
The general Euler-scale hydrodynamic argument presented in this paper supports the assumption that N -point correlation functions vanish, under scaling by λ, as λ 1−N , as per the formula (1.5). In particular, two-point functions vanish as 1/t at large times. This is clear in various formulae established, for instance in the homogeneous case (3.36), in free models (4.19), in the partitioning protocol (5.6), and in the long-time asymptotics (5.10). However, in these formulae, the quantity |∂ θ v eff (θ )| appears in denominators, evaluated in particular states and at particular values of θ (for instance, |∂ θ v eff (ξ, θ )| in the state at ray ξ of the partitioning protocol, evaluated at θ such that v eff (ξ; θ ) = ξ). This quantity may vanish if the effective velocity is not strictly monotonic with respect to the rapidity, and may thus lead to singularities (except, in some cases, if there is zero density of quasi-particles at this rapidity, or, in fermionic systems, of quasi-holes). In such situations, the asymptotic formulae we show do not apply, and we expect a modification of the large-time limit, naively as 1/ t. The physical intuition is that, if there is a finite quasi-particle density at a rapidity for which ∂ θ v eff (θ ) = 0, then, as the effective velocity is stationary in θ , there is an accumulation of quasi-particles around this effective velocity. If, for instance, the observation ray in the partitioning protocol is along this velocity, then this accumulation may increase the correlation. This might happen at the boundary of the "light cone" emanating from the connection point, if there is a maximal velocity. Similar effects might appear in fully inhomogeneous situations if the rapidity derivative of the characteristic function u(x, t; θ ) vanishes, as again singularities may occur for instance in (3.22) and (3.25). It would be interesting to further investigate this aspect.
Comparing exact hydrodynamic predictions for two-point functions with numerics is a very important problem. Steps forwards are made in this direction in [16], where the classical sinh-Gordon model is studied, both for one-point functions in the partitioning protocol and correlation functions in GGEs. In particular, the classical spectral functions for an infinite family of vertex operators are evaluated, and numerical comparisons are made. Comparison with quantum field theories are however more challenging.
It would be interesting to investigate if hydrodynamic ideas provide more than the Eulerscale part of correlation functions, the least decaying part found along ballistic rays. For instance, the recent works [76][77][78][79] suggest that it is possible to combine hydrodynamics with a more detailed knowledge of local observables in order to go further.
Other ways of deriving Euler-scale correlation functions in homogeneous cases are based on form factors. This was done in [69] based on form factors obtained in [93]. The form factor techniques of [53,56] might also be applicable as explained in the Remark in Subection 3.4. We note that using the general results of [58] for space-like two-point functions in arbitrary homogeneous GGEs, one may combine this with the spectral function method of subsection 3.4 in order to get various configurations of dynamical higher-point functions. It would be interesting to see if form factors can be used to derive results in the inhomogeneous situations considered here. It would also be interesting to obtain correlations in situations with evolution in weakly varying external potentials or temperature fields. The GHD theory for such situations was developed in [12], however the equivalent of the solution by characteristics (2.23) has not yet been written. We leave this for future works.

Acknowledgments
I am grateful to Alvise Bastianello, Olalla Castro Alvaredo, Jacopo De Nardis, Jérôme Dubail, Balázs Pozsgay, Herbert Spohn, Gerard Watts and Takato Yoshimura for discussions, comments and encouragement, and for collaborations on related aspects. I am grateful to the Institut d'Étude Scientifique de Cargèse, France and the Perimeter Institute, Waterloo, Canada for hospitality during completion of this work. I thank the Centre for Non-Equilibrium Science (CNES).

A Invertibility of the function u(x, t; θ )
Here we show that u (x, t; θ ) < 0 if the effective velocity is monotonic in the rapidity. This is natural: recall that the function u(x, t; θ ) represents the position, at time 0, from where a quasi-particle trajectory of spectral parameter θ would reach the position x at time t. Therefore, since the effective velocity is monotonic with the rapidity of quasi-particles, a positive change of the velocity associated to θ will occasion a negative change of the initial position of the trajectory that reaches the space-time point (x, t).
More precisely, for the formal proof, we assume that the effective velocity is monotonic (v eff ) (θ ) > 0 (this is the case in many situations in QFT for instance, see [7]), and that (2.23) has a solution. We also recall that u(x, 0; θ ) = x.

B Propagator and derivation of two-point function formulae B.1 Main formulae
Here we present the derivation of formulae (3.16) and (3.17) using the technique explained in Subsection 3.1. Formulae (3.18) is obtained by symmetry, and we note that (3.17) and (3.18) agree with hydrodynamic projection principles. Formula (3.19) is then obtained by using hydrodynamic projections as explained in the Subsection 3.3.
We start with the one-point function. We use the general formula (3.5). Differentiating with respect to β j ( y), we find In order to evaluate η j (x, t; y; θ ), we use the solution (2.23). Two terms occur: the first is the derivative of n 0 (u; θ ) with respect to β j ( y) at u fixed, the second involves the derivative of u(x, t; θ ) with respect to β j ( y). Using (2.18), the first term gives The second term is evaluated using (3.10), giving In this latter expression, the derivative of u(x, t; θ ) occurs. This cannot be evaluated explicitly, but we can obtain an integral equation involving it by differentiating the second equation of (2.23): (B.5) Let µ be again a generic parameter on which a GGE state n may depend. Then the following general formula holds, for spectral functions g(θ ): (B.6) Therefore, using (2.18), and definition (B.2) implies where u = u(x, t; θ ). Replacing η j (x, t; y; θ ) = n t (x; θ ) f (x, t; θ ) Γ ( y,0)→(x,t) h dr j ( y, 0) (θ ), which follows from the definitions (3.7) and (B.2), along with the chain rule and (2.18), we obtain the defining integral equation (3.12).
Formula (3.17) is obtained in a similar way starting from the one-point function 〈j i 〉 [n] , giving

B.2 Indirect propagator
We show the integral equation (3.21) for the indirect propagator ∆ ( y,0)→(x,t) g (θ ). According to the second term on the left-hand side of (3.12), we need to evaluate the star-dressing of the function ρ s (z, t; θ ) f (z, t; θ )δ( y − u(z, t; θ ))g(θ ) as a function of θ . For this purpose, we note that the application T n t (z) on it, which is required as per definition (3.9), gives where the root set θ (z, t; y) is defined in (3.24). This then needs to be dressed, but the only dependence in θ is via the differential scattering phase ϕ(θ , α). Consider the dressed scattering operator (3.23). In components, it is Using (2.23), this can be simplified slightly to (3.22).

C Verification of the Leclair-Mussardo spectral function
In this appendix we verify that the Leclair-Mussardo spectral function (3.39) indeed reproduces the conserved-density and conserved-current spectral functions (3.34), when the diagonal matrix elements involved in the Leclair-Mussardo formula are specialized to those of conserved densities and currents. Here for simplicity we specialize to the sinh-Gordon model with unit mass. We refer to the explanations in [43,44,94] for the initial studies, and to [7, App D] for the explicit diagonal matrix elements af all conserved densities and currents. The results are F q i k (θ 1 , . . . , θ k ) = ϕ(θ 1,2 ) · · · ϕ(θ k−1,k ) h i (θ 1 ) cosh(θ k ) + permutations, (C.1) and and similarly for j i . Explicitly, we have Identifying the *-dressing operation, this is Therefore, using cosh dr (θ ) = 2πρ s (θ ), we find that (C.3) reproduces the first equation of (3.34). A similar calculation reproduces the second.

D The characteristics in the partitioning protocol
Recall the partitioning protocol, and in particular (5.1) and (5.2). In this case, the solution (5.2) is obtained without the use of characteristics (2.23). Nevertheless, it is useful, when analyzing correlations, to have an understanding of the function u(x, t; θ ). First, we show that we can re-write where the functionũ depends on x, t only through the ratio ξ = x/t. Indeed, thanks to the exact solution (5.2) it is clear that the effective velocity v eff (ξ; θ ) likewise only depends on ξ.
With the change of variable, the differential equation in (2.24) leads to Further, the initial condition becomes the asymptotic condition lim t→0 (t/x)ũ = 1, Since the equation and asymptotic condition only involve the variable ξ, the solution likewise only do. Next, we may solve exactly this equation. Integrating over ξ, we obtaiñ This can be re-written, for any A(θ ) and B(θ ), as The solution of course does not depend on ξ 0 . We may therefore take the limit ξ 0 → −∞ on the right-hand side. Choosing A(θ ) = 1, we may use (D.3) as well as the fact that lim η→−∞ v eff (η; θ ) = v eff L (θ ) (that is, the limit is finite) in order to see that the right-hand side has a limit that gives Since u(x, t; θ ) is strictly increasing with x (see (2.27)), it has at most a single zero as function of x. Comparing the exact solution (5.2), (5.1) with the first equation in (2.23), we find that this zero must be at the solution to the equation x/t = v eff (x/t; θ ). This has the simple physical interpretation that the quasi-particle whose trajectory reaches the point (x, t) at an effective velocity equal to the ray ξ, is the one that goes along the ray ξ and thus originates from x = 0 at t = 0. Let us define the function ξ (θ ) by this solution This gives a convenient explicit form of the functionũ(ξ; θ ). Note that the integrand on the right-hand side of (D.9) is composed of two terms both of which have a unique pole at the same position. Since ∂ ξũ (ξ; θ ) exists and is (finite and) nonzero by (2.27), thenũ(ξ; θ ) must have a simple zero at ξ = ξ (θ ). This implies that the poles of the two terms in the integrand cancel each other. Therefore Taking the θ -derivative of (D.7), we conclude that .

E Long time limit E.1 A proof of the emergence of the partitioning solution
We make the assumptions stated in the first paragraph of Subsection 5.3, and only assume, instead of those made in the second paragraph, that the limit lim t→∞ n t (ξt; θ ) exists and is of the form n(ξ; θ ). We also assume that the state density is uniformly bounded away from zero and infinity in space time, and we assume that the integral converges to that of the pointwise limit of its integrand at large t. We show (5.11) as follows. Consider the integral equation (2.23). Using the stated assumptions, we have that lim t→∞ ρ s (ζt, t; θ ) is of the form ρ s (ζ; θ ), and that lim t→∞ ρ s (ζt, 0; θ ) = ρ sgn(ζ) s (θ ) and we find, to leading order in t, where ξ = x/t. Since the state density is (uniformly) positive, in order for the equality to hold u(ξt, t; θ ) must scale proportionally to t at large t, except for the possible values of (ξ, θ ) where the left-hand side vanishes. Given a ξ, only a finite number of values of θ might make this happen. Since these do not affect spectral integrals, they do not affect the evaluation of the state density from the occupation function. Therefore, using (5.9) we find

E.2 The function r( y; θ ) and its integral equation
Choose ξ = ξ (θ ), for spectral parameter θ = (θ , a). We assume that (v eff ) (ξ; θ ) = 0. Define θ t ( y; θ ) = (θ t ( y; θ ), a) some element in θ (ξt, t; y); if the effective velocity is monotonic with respect to the rapidity, then this is the unique element with particle type a; but otherwise it is an element which continuously depends on t, and which, given θ , can generically be made unique for large enough t by its large-time limit. That is, we have u(ξ (θ )t, t; θ t ( y; θ )) = y. Taking the large-t limit, we getũ(ξ; θ ∞ ) = 0, and thus θ ∞ = lim t→∞ θ t ( y; θ ) = θ , and the choice of θ determines (generically) the element θ t (t; θ ). The function r( y; θ ) is defined as r( y; θ ) = lim t→∞ tũ(ξ (θ ); θ t ( y; θ )). (E.5) We assume that the limit defining r( y; θ ) exists and is finite, and that it is differentiable with respect to y. Clearly u (ξt, t; θ t ) = (∂ θ t /∂ y) −1 . Since r( y; θ ∞ ) is finite, and since, by the assumption that (v eff ) (ξ; θ ) = 0, the functionũ(ξ; θ ) has simple zeroes in θ , then θ t approaches the zero θ ∞ with corrections of order t −1 : Let us assume that the corrections o(t −1 ) are smooth enough in y. Then the term displayed gives the correct variation of θ t ( y; θ ) with respect to y at fixed θ to leading order in t −1 : where we used (D.14). Now consider the integral equation (2.23). Subtracting that for tũ(ξ; θ t ) from that for u(ξt, t; θ t ), and using the fact that the states n(z, t; θ t ) and n(z/t; θ t ) have the same asymptotic at large distances, we get We may take the large-t limit. Assuming that the initial state densities are continuous in rapidity, we therefore find (changing integration variable on the left-hand side) Here we have been careful not to simply replace, on the left-hand side, the integrand by its limit. The physical meaning of the left-hand side of (E.9) is as follows. First observe that thanks to the conservation equation ∂ t ρ s (x, t; θ ) + ∂ x (v eff (x, t; θ )ρ s (x, t; θ )) = 0 [7,8]. This represents the inflow and outflow at the asymptotic boundaries of the system. The same t-dependence occur for ∞ −∞ dz ρ s (z/t; θ ), as the partitioning protocol is based on the same states at its asymptotic boundaries. Therefore ∞ −∞ dz ρ s (z, t; θ ) − ρ s (z/t; θ ) does not depend on t. Taking t → 0, we see that it is simply equal to the total difference between the initial state density, which effectively regularizes the partitioning protocol, and the partitioned initial state with a discontinuity at the origin; this difference is finite by (5.9). Taking the large-t limit, it is clear that we cannot simply take the limit on the integrand, as we would get zero. As time evolves, the difference between the two initial conditions is redistributed in space: the difference between the integrands goes to zero, but the integrated difference does not. The limit in (5.12) measures how much of the initial state density difference has been transferred to the left of the ray ξ. This quantity should become constant in time: we would expect it to be the portion of the initial, finite difference carried to the left of ξ by all quasi-particles of effective velocities allowing them to cross the ray. It is a nontrivial quantity, as it depends on the details of the initial condition n 0 (x; θ ).
As explained in the main text, this the implies (5.14), which states monotonicity of r( y; θ ); thus the assumption is indeed consistent.

E.3 Derivation of the main result
There are two contributions to A i j (ξ; y): from the direct and the indirect propagators. Consider first that from the direct propagator, the first term on the right-hand side of (3.25).
We need to evaluate the long-time asymptotic of the spectral derivative of u(x, t; θ ). By (D.1), u(x, t; θ ) grows with t as u(ξt, t; θ ) ∼ tũ(ξ; θ ) (t → ∞) (E. 18) for any θ such thatũ(ξ; θ ) = 0. This simply means that the point y = u(x, t; θ ) from which a quasi-particle reaches x after a long time t is generically very far away from the origin. Thus the spectral derivative should also grow with t. However, we cannot simply take the θ derivative of (E.18) in order to obtain the large-t asymptotic of u (ξt, t; θ ). This is because the finite correction to (E.18) also has a very large derivative. This finite correction occurs when the spectral parameter θ is very near to that of a quasi-particle traveling along the ray ξ in the partitioning protocol; that is, very near to satisfying v eff (ξ; θ ) = ξ, equivalentlyũ(ξ; θ ) = 0, or according to the notations introduced, θ ∈ θ (ξ; 0) = θ (ξ). If θ approaches such a point as time grows, the point y may stay finite. But at long times, a small change of θ (of order 1/t) will occasion a large change of y (of order 1), because the trajectory depends on the precise structure of the state in finite regions around the origin, and a finite region is spanned by a small change of θ . This contribution is in fact immediate to evaluate form the result of the Appendix E.2. Indeed, from (E.7) we have u (ξt, t; θ t ( y; θ )) = −t V (θ ) (v eff ) (ξ; θ ) ∂ r( y; θ )/∂ y + o(t), (E. 19) where ξ = ξ (θ ).
In order to go further, we need to understand how other functions behave when evaluated at (ξt, t; θ t ( y; θ )). Clearly, the occupation function n t (ξt; θ t ( y, θ )) does not tend to its partitioning value n(ξ; θ ) = n sgn(ũ(ξ;θ )) 0 (θ ), but rather equals n 0 ( y; θ ). However, the principle we will use below is that any dressed quantity, h dr (ξt, t; θ t ( y; θ )), tends to its partitioning value h dr (ξ; θ ) at large t. This is because dressing involves spectral integrals of n t (ξt; α), and the set of values of α around θ t ( y; θ ) for which n t (ξt; α) is significantly different from its partitioning value becomes of measure zero, at large t, under the dα measure. This is the same effect as that explained in Appendix E.2, and the above principle was used there for the state density ρ s = (p ) dr /(2π).
Therefore, combining (E.19) with (5.14), the first term in (3.25) gives the following contribution to A i j (ξ; y): Remark that this can be seen as coming from the integral Consider now the contribution from the indirect propagator: the second term in (3.25). In order to have an intuition of its contribution, recall that the effective acceleration a eff [n 0 ] (z; θ ) is zero, except for z in a region around the origin which can roughly be considered as finite.
With an argument similar to that made above and in Appendix E.2, since u(ξt, t; θ ) diverges proportionally to t as per (E.18), it is clear that a eff [n 0 ] (u(ξt, t; θ ); θ ) vanishes for almost all values of θ , except for a small region, whose extent decreases as t −1 , around the zeroes of u(ξ; θ ). Therefore, any integral of the form S dθ a eff [n 0 ] (u(ξt, t; θ ); θ ) g(θ ), for bounded spectral function g(θ ), decreases proportionally to t −1 and is supported on the point set θ (ξ). Thanks to (3.21), the indirect propagator has an overall factor a eff [n 0 ] (u(ξt, t; θ ); θ ), wherefore the indirect propagator contribution -the second term in (3.25) -is an integral of the above form. Thus, in order to obtain the leading t −1 decay of the correlation function, we only need to keep terms that stay finite at large t on the right-hand side of the integral equation (3.21).
In order to determine the finite contribution on the right-hand side of (3.21), consider first the source term (3.22). It has itself two contributions. For the first, we do the change of variable z = ηt to write it as The second term in the source (3.22) is clearly finite. Next, consider the second term in (3.21), the integral of a star-dressed quantity involving the indirect propagator itself. From the definition (3.9), a star-dressed quantity can be written as a series of terms each involving at least one spectral integral. Since, as argued above, the indirect propagator has an overall factor a eff [n 0 ] (u(ξt, t; θ ); θ ), and spectral integrals involving such factors decrease as t −1 , we conclude that the second term in (3.21) also decreases as t −1 . Therefore, for the purpose of evaluating the leading decaying term of the correlation function, we may make the replacement ∆ ( y,0)→(x,t) g (θ ) → 2πa eff − Θ(u(ξt, t; θ ) − y) ρ s ( y, 0) f ( y, 0)g * dr ( y, 0; θ ) .
We therefore get the following contribution to A i j (ξ; y): We can then sum this contribution with (E.20) to get the full coefficient. Before going further, however, let us note the apparent lack of space parity symmetry in the expression (E.28), both in the integral I( y; γ), and in the integral ξ −∞ dη. We do not assume the model to be parity symmetric, however what we note here is that the expression treats left and right regions of space differently, independently of the properties of the model. This is due to the lack of a manifest parity symmetry in the solution by characteristics (2.23), where the integral is chosen to start at a left asymptotic stationary point. As mentioned in [33], this is a conventional choice, and a similar formula can be obtained by integrating towards a right asymptotic stationary point instead. It is not too difficult to obtain the result with this different choice: If the asymptotics n ± 0 (θ ) are different, then the expression within the large parentheses cannot be zero: as a function of y, the first terms has a jump at y = y (α), while the second term is continuous. Therefore we conclude that I(γ) = 0. (E.32) Since I(γ) is continuous as a function of the asymptotics n ± 0 (θ ), the case of equal asymptotics is obtained by taking the limit, thus also giving 0.