Euler-scale dynamical correlations in integrable systems with fluid motion

We devise an iterative scheme for numerically calculating dynamical two-point correlation functions in integrable many-body systems, in the Eulerian scaling limit. Expressions for these were originally derived in Ref. [1] by combining the fluctuation-dissipation principle with generalized hydrodynamics. Crucially, the scheme is able to address non-stationarity, inhomogeneous situations, when motion occurs at the Euler scale of hydrodynamics. Using our scheme, we study the spreading of correlations in several integrable models from inhomogeneous initial states. For the classical hard rod model we compare our results with Monte-Carlo simulations and observe excellent agreement at long time scales, thus providing the first demonstration of validity for the expressions derived in Ref. [1]. We also observe the onset of the Euler-scale limit for the dynamical correlations.

In many-body systems, dynamical correlation functions -measuring the correlations between observables at different space-time points -are of particular interest. They provide a powerful way of characterizing quantum and classical systems, encoding the emergent degrees of freedom and their dynamics. Unfortunately, they are difficult both to calculate theoretically from first principles, and to extract experimentally [10,12,53]. The frameworks of linear and nonlinear fluctuating hydrodynamics [54], and the theory of hydrodynamic projections [44,55], provide powerful methods able to access the large-scale behaviours of dynamical correlation functions. By focusing on the Euler-scale hydrodynamic equations and the propagation of long-lived modes, they allow one to extract exact asymptotic expressions for correlation functions along the propagation of such modes. However, these standard methods, as developed until now, are limited to correlation functions in homogeneous and stationary states: they are based on (non)linear response mechanisms with respect to such states. Yet, out of equilibrium, the state of the many-body system is often inhomogeneous and non-stationary. In fact, one expects that, at large enough times, before thermodynamic stationarity is reached, local relaxation occurs and there is a long period of large-scale motion, well described by hydrodynamics. A natural question arises as to the correlations within such states.
The theories of inhomogeneous Luttinger liquids [56][57][58][59][60][61][62] and, more recently, quantum GHD [48] are able to account for quantum critical correlations in such cases, which dominate at zero temperature or more generally in zero-entropy states (not necessarily of Gibbs form). At finite entropy, the dominant correlations are instead due to classical fluctuations 1 , both in quantum and classical systems. For these, a combination of the fluctuation-dissipation principle and GHD lead to a recursive procedure for generating n-point correlation functions at the Euler-scale [1]. Crucially, the method gives results valid also when the system displays nontrivial Euler-scale motion. Given the universality of GHD, the results obtained should be applicable to any model obeying the hydrodynamic equations. Through the developed technique, exact dynamical two-point correlations for a large class of local fields were derived. However, numerical checks were made only in stationary states [43,63]. The much more nontrivial part of the theory of Ref. [1] is that dealing with correlations within fluids with Euler-scale motion, where correlations depend on the full initial profile of the GHD fluid. The results take the form of a quite involved set of nonlinear integral equations, which were never solved numerically, let alone verified against microscopic simulations.
In this paper, we develop a numerical scheme for calculating dynamical two-point correlations in Euler-scale GHD, in the full generality of the theory of Ref. [1]. We use the scheme to describe the propagation of correlations from inhomogeneous states in the Lieb-Liniger model, the relativistic sinh-Gordon model, and the classical hard rods gas. To confirm the validity of our method, we compare our results for the hard rod model to Monte-Carlo simulations and examine the time required to reach the Euler-limit of correlations.
The paper is organized as follows. In Section 2, we briefly review the basics of GHD, while the correlation functions originally derived in Ref. [1] are discussed in Section 3. Next in Section 4, we calculate and analyze the spreading of correlations in three distinct systems: The homogeneous Lieb-Liniger model, bump releases in the classical hard rod model, and the partitioning protocol in the Lieb-Liniger model. Finally, we conclude in Section 5.

Summary of GHD
Although the theory of generalized hydrodynamics is still relatively young, several works have provided reviews, see for example Refs. [64,65] and the lecture notes [66]. In this section we simply reiterate some of the central concepts of GHD, which later will be relevant for computing dynamical correlation functions. Note that various notational conventions have been used in the literature; here we will follow that of [65].
Generalized hydrodynamics describe the transport properties of integrable system, which due to their infinite set of conservation laws exhibit dynamics very different from conventional hydrodynamics [13,32,33,67,68]. The theory employs the language of the thermodynamic Bethe ansatz (TBA), where the full thermodynamics of an integrable system is encoded in a root density, ρ p (λ). The root density can be interpreted as a density of quasi-particles parametrized by the rapidity λ [19,22,69,70]. Equivalently, the system can in this stage be described by the generalized Gibbs ensemble (GGE) [15,16], which is completely specified by 1 Recall, though, that correlations due to classical fluctuations, that is, to the statistical distribution of states, keep a memory of quantum effects in quantum systems: the statistics of the underlying degrees of freedom affect the distribution function. the knowledge of the so-called pseudoenergy function ε(λ), that is in turn determined by the non-linear integral equation where T is the differential two-body scattering phase given by T (λ, λ ) = ∂ λ Θ(λ, λ )/2π and Θ(λ, λ ) is the scattering phase; F(ε) is named free energy function and it depends on the statistics of the quasiparticles of the theory. For fermions F(ε) = −log(1 + e −ε ) (as in the case of the Lieb-Liniger model and the sinh-Gordon) and F(ε) = −e −ε for classical particles (as for the hard rods gas); w(λ) is named source term and it is given by where h i (λ) is the single particle eigenvalue of the conserved charge Q i and β i the corresponding Lagrange parameter. The thermodynamic Bethe ansatz (and GGE) is in principle only valid on stationary states, however, the root density can be propagated according to Eulerian hydrodynamical equations, whereby dynamics of the system can be treated. This realization sparked the theory of generalized hydrodynamics, which considers a system consisting of locally homogeneous space-time fluid-cells exhibiting only small variations to neighbouring cells. Within these cells the microscopic dynamics establishes local equilibrium at time-scales faster than the global dynamics, whereby the system remains in a quasi-stationary state. Thus, the thermodynamic Bethe ansatz still applies, although the root density is now both weakly time and space dependent, ρ p (λ) = ρ p (x, t; λ) [32,33].
In the absence of any inhomogeneous external fields, the dynamics of the system are given by a simple, Eulerian fluid equation [32,33] where ϑ t (x; λ) ≡ ϑ(x, t; λ) is dubbed the filling function and, like the root density, encodes the full thermodynamics of the system. It is related to the root density via the non-linear relation The local velocity field v eff (λ) in Eq. (3) has been first defined in Ref. [71] and it specifies the propagation velocity of the quasi-particles where (λ) and p(λ) are the one-particle energy and momentum respectively. Here, the superscript dr denotes that the quantities have been dressed, encoding a modulation of the quantity induced by interactions between the quasi-particles. The effective velocity encodes the Wigner delay time, which is associated with the resulting phase shifts of elastic collisions in integrable systems [72]. In the absence of interactions the particles propagate with the group velocity v gr (λ) = ∂ λ /∂ λ p. The dressing of some quantity h, e.g. the single particle eigenvalue h i (λ) in Eq. (2) of some conserved charge Q i , is defined through the integral equation The hydrodynamic principle states that the system during its evolution via Eq. (3) remains in a quasi-stationary state. Thus, the entire TBA framework can be applied at any point in time. For instance, given the root density one can calculate the density of the i'th conserved charge and its associated current via the relations where h i is the one-particle eigenvalue of the charge Q i . The exact average currents fall outside of the historically developed TBA machinery, but are derived in [32,[73][74][75][76][77][78][79] Lastly, Eq. (3) admits a solution by a characteristic function, U(x, t; λ), encoding the inverse trajectories of the quasi-particles [37,65,80]. Given the characteristic, one can directly relate the evolved state to the initial state via the relation The characteristics follow the same hydrodynamical equation as the filling function whereby they can be efficiently computed alongside the filling function. The characteristic function also solve a system of nonlinear integral equations [80] where time enters as a fixed parameter. This forms the basis for the recursive method leading to the exact Euler-scale dynamical correlation functions in GHD [1], which we make use of below. Note that some TBAs admit multiple root densities (that is, multiple quasi-particle species). In this case all integrals over the rapidity transform into integrals over the rapidity of each root density, where the contribution from each root density is added together.

Exact Euler-scale dynamical two-point correlations
Euler-scale hydrodynamics describes the flow of conserved charges between locally homogeneous space-time fluid cells, whose variations to neighbouring cells are sufficiently small. Similarly, an Eulerian scaling limit for correlations exists, which is defined as the large-scale limit of connected correlation functions [1,55]. Thus, the Eulerian scaling limit probes only long wavelengths, which is consistent with the small variations between fluid-cells assumed in GHD. The precise definition of the Eulerian scaling limit for correlation functions needs some care, and in some cases nontrivial averages over fluid cells are necessary [81]; in the model analysed below, a simple averaging in space is sufficient, as described in Subsection. 4.2.
In Ref. [1], a recursive procedure for generating such correlations was obtained based on linear responses of ϑ t (x; λ) to variations in the initial state. In this method, the Eulerian dynamical correlation functions are obtained from a linear response analysis, a generalization of the fluctuation-dissipation theorem. Importantly, due to the exact solution to the problem of characteristics [80], in integrable systems, this allows one to extend the method to situations with large-scale, inhomogeneous initial states.
A main feature of GHD is the quasi-stationary fluid state, implying that at every slice of time t, the full thermodynamics of the system is determined by a weakly spatially inhomogeneous root density. Equivalently, one can consider each fluid cell being described by a separate GGE. Thus, all equal-time, space-separated connected correlation functions vanish (that is, decay exponentially fast with a microscopic-scale correlation length). However, over time the propagation of quasi-particles causes quantities in separated fluid cells to become correlated in a non-trivial manner. Hence, dynamical correlations at the Euler scale can be viewed as initial delta-functions correlations, which over time ballistically spread and propagate throughout the system. This intuitive picture is reciprocated in the expression for the exact two-point Euler-scale correlation function, wherein the propagator Γ (y,0)→(x,t) (λ, λ ) [1] encodes the ballistic propagation of local quantities between two points in space-time. This, in general, extends to inhomogeneous states the concept from hydrodynamic projections [44].
One of the main results of Ref. [1], is the derivation of an exact formula for two-point correlations of generic local observables Here, f (λ) is the statistical factor of the model. For models with fermionic quasi-particle statistics f = 1 − ϑ(λ) (the case of Lieb-Liniger and sinh-Gordon), while for classical particle models f = 1. Meanwhile, the square brackets in Eq. (11) denote the contraction Lastly, the field V O in Eq. (11) is the one-particle-hole form factor of the operator O. It is a functional of the filling function such that The form factors must be worked out for every operator individually. For charge densities and associated currents they are very simple [44] however, other observables such as vertex operators in the sinh-Gordon model have form factors with more complicated expressions (see Refs. [1,82] and Appendix C). The propagator, Γ (y,0)→(x,t) (λ, λ ), describes how the local quantity V O (y, 0; λ ) travels through the system on a given trajectory, until it reaches the location x at time t. The propagator itself can be split into two terms where each term has a clear physical interpretation. The first term denotes the direct propagation of the quantity carried by the quasi-particles with the inverse trajectory U(x, t; λ). Thus, of all the quasi-particles found at (x, t), only those with the right rapidity have arrived from the point (y, 0). Meanwhile, the second term is dubbed the indirect propagator, as it describes modifications to the correlations due to perturbations of the quasi-particle trajectories from local inhomogeneities at (y, 0). Hence, all rapidities can in principle contribute to the indirect correlations.
Inserting Eq. (15) into the two-point correlation formula of Eq. (11) yields where the set λ (x, t; y) = {λ : U(x, t; λ) = y} contains only the rapidities of quasi-particles directly propagating the correlations. While the direct correlations given by the first term of (16) are relatively simple to evaluate, the indirect term poses more of a challenge. The indirect propagator follows the linear integral equation where the field D 0 encodes the degree of inhomogeneity of the initial state (it is the "effective acceleration" [64]) and the so-called source term reads Note, if the state is homogeneous D 0 vanishes, thus eliminating any indirect correlations. In the equations above h * dr (λ) = h dr (λ) − h(λ) and Θ(x) is the Heaviside function, while x 0 is an asymptotically stationary point, which must be chosen such that ϑ s (x; λ) = ϑ 0 (x; λ) for x < x 0 and s ∈ [0, t] [80]. Thus, x 0 denotes the boundary for which disturbances of correlations have yet to spread within the time t. One could think of x 0 = −∞, although for numerical simulations it is set as the first spatial gridpoint, which much be chosen sufficiently far away from the point y or any inhomogeneities. Solving Eq. (16) requires mostly quantities already available from the TBA and GHD frameworks, however, currently no numerical solution of Eq. (17) and (19) have been shown. In Appendix A we report an iterative scheme for calculating the indirect propagator.
Finally, one should note that in Euler-scale GHD, one does not fully specify the observables whose correlations are evaluated: only partial information is given. For instance, conserved densities are ambiguous. Indeed, the GGE density matrix enables the exact calculation of thermodynamic averages, but contains information only of the conserved charges, Q i [15]. Thus, the conserved charge densities, q i , from which the Euler-scale two-point correlation functions are derived, are defined only up to a total spatial derivative. However, in the Eulerian scaling limit any derivative corrections to q i are expected to be vanishing small, since the large-scale limit only probes long wavelengths [1].

Numerical calculation of correlations
In this section we calculate dynamical, Euler-scale two-point correlation functions by numerically solving Eq. (16). To demonstrate properties of correlations at the Euler scale we examine three different scenarios, whose hydrodynamical properties have already been well studied: First, we calculate the spreading of correlations in a homogeneous system. The absence of inhomogeneities drastically simplifies the problem, as the indirect contribution to the correlations vanish. Next, we study how correlations spread during a bump release protocol in the classical hard rod model. As this is a classical model, we can simulate it via Monte-Carlo methods and microscopically measure the spreading of correlation functions, thus giving us the opportunity to test the equations for Euler-scale correlations. Lastly, we examine the iconic partitioning protocol, where two homogeneous, semi-infinite systems initially at different temperatures are joined together at t = 0. Although partial analytic predictions for the correlations in such a setup were made in Ref. [1], these in fact contain inaccuracies. Our numerical analysis unveils the full dynamics even at short time scales.
The exact numerical procedure for solving Eqs. (17) and (19) is fully detailed in Appendix A. Aside from the propagator, Γ (y,0)→(x,t) (λ, λ ), all other quantities of Eq. (16) are readily available via iFluid, an open-source framework for GHD calculations [65]. As part of this work, the code for calculating the propagator and two-point correlations have been integrated as a standalone module in the framework [83].

Homogeneous state
The spreading of Euler-scale correlations in a homogeneous system, ϑ t (x; λ) = ϑ(λ), is particularly simple as D 0 = 0, causing the indirect propagator (17) to vanish. Furthermore, the velocity of the quasi-particles is spatially independent, whereby the characteristic solution to Eq. (3) becomes U = x − v eff (λ)t. Therefore, the the full propagator (15) reduces to and the dynamic two-point correlation function of the zeroth charge density, q 0 = n, for y = 0 becomes In Eq. (21), λ (ξ) is the set of solutions to the equation v eff (λ) = ξ = x/t. Thus, the correlations spread at the same velocity as the quasi-particles move, while they diminish over time as t −1 . This formula was obtained in [44], and follows from a direct application of hydrodynamic projection methods. The decay in t −1 is a consequence of the continuum of hydrodynamic modes (parametrised by λ) on which projection occurs -and thus this is a special property found in integrable models. Note, in models like the Lieb-Liniger and sinh-Gordon model, v eff (λ) is a monotonically increasing function of λ. Hence, for any combination (x, t) the set λ (ξ) will contain only one element. Computing dynamical two-point correlation functions via Eq. (21) is remarkably straightforward, as the expression can be evaluated using only information available from the TBA without performing any hydrodynamical evolution of the system. Simulation parameters can be found in the main text.
In Fig. 1, density-density correlations calculated via the full formula (16) and simplified formula (21) are compared. The simulation was carried out for the Lieb-Liniger model with inverse temperature β = 1, interaction strength c = 1, and chemical potential tuned to a linear density of n(x, t) = 0.5. Fig. 1(a) depicts n(x, t)n(0, 0) Eul and illustrates the aforementioned interpretation of the Euler-scale dynamic correlations; an initial delta function which spreads ballistically throughout the system. Since the quasi-particles move with the same velocity regardless of position and time, the solutions of the hydrodynamic equation (3) are constant on the ray ξ = x/t, even at short timescales. This is exemplified in Fig. 1(b), where the two-point correlation function scaled by the time, t, is plotted as function of the ray, ξ. Here, correlations calculated via the full and the simplified formula overlap perfectly, as one would expect. The shape of the two-point correlation profile features a dip towards the center originating from the statistical factor f (λ). In the TBA of the repulsive Lieb-Liniger model, the quasi-particles are fermions (despite the model describing a Bose gas). Hence, the statistical factor reads f (λ) = 1 − ϑ(λ) and the filling function is capped at one. For sufficiently cold temperatures, a Fermi sea of quasi-particles can form at low rapidities and create a barrier for other quasi-particles trying to pass through. While the system studied here is not cold enough to form a full Fermi sea, its filling function at low rapidities is still somewhat close to unit. Therefore, the propagation of correlations at lower rapidities (and by extension low values of ξ) is limited causing the dip visible in Fig. 1(b).

Bump release and Monte-Carlo comparison
Next, we study the spreading of correlations in the classical hard rod model. This model describes classical rods of length a that propagate freely except for elastic collisions, at which rods exchange their velocities. The TBA functions describe the velocity tracers, the tracers of rods at a given velocity λ. These propagate with linear trajectories interrupted by actual jumps occurring when collisions happen. The main ingredients of the TBA description of the model, first developed in Refs. [25,44], are reported in the Appendix C.3. Fundamentally, thanks to its classical properties, we can directly measure the spreading of connected correlation functions via classical Monte-Carlo simulations and therefore compare with the Euler-scale formulas.
For this demonstration, we turn to another well-studied protocol, namely the release of a density bump, initially located around x = 0, created by an inhomogeneous temperature profile. In addition, we also consider the more intricate case of the release of two bumps that initially do not overlap on top of a thermal background. Both setups are akin to what was studied experimentally in Ref. [13]. In the Appendix D further results for the release of two density bumps are presented. The initial state is a thermal GGE identified by the source term w (th) (x, λ) = β(x)λ 2 /2, cf. Eqs. (1), (2) and h(λ) = λ 2 /2 the single particle energy eigenvalue for a Galilean model in Eq. (6), with β(x) for the two bumps problem while for the single bump case a single Gaussian profile is considered z is the parameter controlling the smoothness of the bump space dependence, ±x 0 the bumps positions (for simplicity we take them symmetric with respect to the origin) and β as , β in are the thermal background and the bump inverse temperatures respectively. The thermal root density ρ (th) p of the initial state at time t = 0 reads where d(β) = 1/ √ 2πβ, whereby the initial linear density of particles n(x, 0) is Here W (z) is the Lambert W function on its principal branch. In the Monte-Carlo simulations, rods are at the initial time t = 0 distributed in space according to Eq. (25) starting from some initial point −L (L > 0), while the velocity of each rod is drawn from a Gaussian distribution with variance 1/β(x) dependent on the point x where the rod is initially located, according to Eqs. (22) (or (23)) and (24). From this initial condition, we then run the deterministic classical dynamics of the hard rods gas. For each sample of the initial condition, the linear particle density n(x, t) (h(λ) = 1 in Eq. (6)), and the density-density connected correlation function t n(x, t)n(0, 0) c multiplied by time t are acquired by counting for each space point x the number of particles in an interval (x − l/2, x + l/2) of length l both at time 0 and after some time t. The average of the aforementioned quantities with respect to many independent realizations of the initial rods' positions and velocities is eventually computed. We stress that only the initial configuration of the particles is random, while the dynamics are completely deterministic. More details about the Monte-Carlo simulations are in Appendix D.
The parameter z has to be chosen big enough such that ρ (th) p (x, 0, λ) is smooth and a sufficiently large number of rods are contained within the bump. In this way one can then expect that the root density ρ   (16)).
β as and is needed to avoid to consider space regions with no particles, which could cause non Eulerian-effects. In particular, β as is taken larger than β in in order for the background density to be smaller than the bump density. Furthermore, since the variance of the rods velocity distribution is 1/β(x), particles from the background density intervals move slower than the ones initially located in the bumps and the dynamics is therefore characterized by the propagation of the particles from the hot high-density bump regions to the cold low-density background. In particular, for short times each of the two density peaks evolves independently of the other, while for large enough times the density in the central background region, around x = 0, increases as a consequence of the arrival of the rods from both the bumps, thereby inducing correlations among the particles coming from the left and the right density peak.  Fig. 2(b) and t = 15, 30 in Fig. 2 (d)) discrepancies between the Monte-Carlo simulations and Euler-scale results are evident. These differences are absent for longer time scales (t = 30, 45 in Fig. 2(b) and t = 70, 90 Fig. 2(d)) so that correlations are well reproduced by their Euler scale limit. The latter aspect is quantified by looking at the relative distance σ between the results of the two methods which is reported in Fig. 2(b),(d).
A subtle aspect of Eq. (16) is the presence of the indirect propagator (17). As mentioned, the direct propagator, the first term in Eq. (16), represents the direct contribution of the normal modes (the quasi-particles), where correlations are due to the direct transport of quasi-particle along their trajectories within the inhomogeneous, non-stationary state. The indirect propagator is a correction to this, and is due to nonlinearity of the GHD equations: in a linear-response picture of the correlation function, it encodes the effects of the local disturbance of normal mode λ on normal mode λ . In our numerical analysis, we observe that this correction is extremely small, the dominant part of the correlation function coming from direct propagation. However, the correction is nonzero, and, as we report in Fig. 2(b), neglecting it renders the agreement with the simulation slightly worse. The subtle effect of indirect propagation is therefore explicitly observed.
We stress that this is the first comparison against numerical simulations of the formulas for the inhomogeneous Euler-scale correlation functions in Eqs. (16), (17) and (19) of Section 3. In the simpler homogeneous thermal framework, Euler-scale correlation functions have been compared in Ref. [81] against Monte-Carlo simulations for the classical sinh-Gordon field theory. In the latter case, results of the simulations oscillate at all times around the GHD predictions and fluid cell averaging is necessary in order to integrate them out. In the present study, on the contrary, the agreement between the classical simulations of the dynamics and the hydrodynamic expression of correlation functions become evident at larger times without the need of any further averaging procedure.

Partitioning protocol
Finally, we turn our attention to the well-known partitioning protocol [32,33], where two homogeneous, semi-infinite systems are stitched together at the point (x = 0, t = 0). The two subsystems have different initial root densities causing a flow of charges between the two subsystems once they are joined together. In this example we study a partitioning protocol in the Lieb-Liniger model, where two subsystems of different temperature β L = 1 and β R = 0.5, but equal linear density n L = n R = 0.5 and interaction strength c = 1, are merged. As shown in [32], the temperature difference alone causes a net flow from the hot side (right) towards the cold (left), as quasi-particles in the hot side generally travel faster.
The standard partitioning protocol features an abrupt transition between the two subsystems, however, this setup is not suitable for numerical calculation of correlations, as the initial inhomogeneity D 0 of Eq. (18) is evaluated via finite difference. Instead, we employ a softened transition achieved via a steep hyperbolic tangent temperature profile. Meanwhile, the chemical potential was adjusted in order to maintain a constant linear density across the system.
In Fig. 3 we have plotted several quantities showing the propagation of density-density correlations, n(x, t)n(y, 0) Eul . Subfigures 3(a) and 3(b) depict the spreading of the direct and indirect correlations for y = 0 respectively. Starting with the direct correlations, they appear very similar to the correlations in the homogeneous setup showcased earlier. This is somewhat expected, as the partitioning setup is (initially) piece-wise homogeneous with linear densities equal of the system in section 4.1. Upon closer inspection of Fig. 3(a), one might notice slightly higher correlations at the negative side (seen more clearly in Fig. 3(h)). This asymmetry reflects the net flow of quasi-particle from right to left, which is further exemplified in Fig. 3(g) showing the formation of the distinct, self-similar density profile as function of the ray, ξ. Moving on to the indirect correlations, we observe that the indirect correlations initially are antisymmetric around x = 0. As time passes, the indirect correlations become more asymmetric due to the flow of particles.
The partitioning protocol is interesting from the point of correlations, since the inhomogeneities are very localized around x = 0, where the subsystems are joined. Therefore, it is interesting to vary y such that it is not necessarily centered on the inhomogeneity. Subfigures 3(c-f) display the time-scaled correlation matrices, t n(x, t)n(y, 0) Eul . Here, one clearly sees how the correlations start as delta functions, whereafter the propagate ballistically throughout the system. As the quasi-particles from the hot subsystem in general move faster, so do the correlations in that side propagate more rapidly. We see this in the correlation matrices, where one half of the correlations extend farther. Furthermore, towards the edges of the correlation matrices the correlations appear homogeneous, whereas around x, y ≈ 0 a transition occurs. These three regions; the left side, the center, and the right side, are further explored in Subfig. 3(h), where n(x, 1)n(y, 0) Eul is plotted for y = −4, 0, 4. The y = 0 profile we have already discussed: it is skewed toward the left due to the particle flow. Meanwhile, the remaining two profile are the result of placing the point y within the homogeneous subsystems. The left subsystem exhibit a correlation profile very similar to the homogeneous system in section 4.1, as the two systems have identical temperatures. Again, the visible dip in correlations in the center of the profile is due to the high filling factor at lower rapidities present in the colder system. Conversely, the right profile exhibit no dip, as the subsystem is too hot to form any Fermi-sea-like quasi-particle distribution, whereby correlations can propagate freely even as low rapidity.

Conclusion
In this paper we have developed an iterative scheme (cf. Section 4 and Appendix A for further details) for finding the correlation propagator necessary for the calculation of exact, dynamical two-body correlation functions in the Eulerian limit (see Section 3 and Eqs (16), (17) and (19) therein). We have applied our scheme to three different setups, whose transport properties have already been well-studied, namely a homogeneous system in Section 4.1, a bump release in Section 4.2, and a partitioning protocol in Section 4.3. Furthermore, the universality of GHD enables our scheme to be applied to most integrable models. Thus, we have studied the spreading of correlations in the Lieb-Liniger model, the classical hard rod model, and the relativistic sinh-Gordon model (cf. Appendix B). In Section 4.2, by comparing for the classical hard rod model against the results obtained via Monte-Carlo simulations (the results are presented in Fig. 2), we have provided the first demonstration of the validity of the formulas derived in Ref. [1] for non stationary and inhomogeneous states. Crucially, we succeeded in explicitly confirming the subtle effect of indirect propagation of correlationscorrelations due to the nonlinearity of GHD, and not directly interpreted as coming from the propagation of normal modes along their curved trajectories in the moving GHD fluid.
From this comparison we are able to observe the onset of the Eulerian limit at longer time scales, while for short times deviations between the classical microscopic simulations and the generalized hydrodynamic predictions for correlation functions are observed. We expect that these discrepancies at smaller time scales can be accounted by considering diffusive terms into the GHD equation (3). Although the effect of diffusive corrections on one-point functions is by now well understood [38][39][40], for dynamical two point correlators in inhomogeneous and non-stationary states, instead, no analytical result is currently available. It would be interesting to investigate this point further in the future. Lastly, our results also enable us to analyze how statistical and dispersion properties of the various models affect the spreading of correlations (cf. Appendix B).
Finally, we point out that our method not only enables future studies of correlation functions, but it also allows one to extend the analysis of Refs. [63,84] for the calculation of the full counting statistics of the time integrated current of ballistically transported conserved quantities to the more complex and interesting case of inhomogeneous states. The authors plan to carry out this analysis in a future publication.
We have made the implementation of our numerical scheme public as a module of the iFluid framework for GHD calculations [65,83].
In order to solve Eqs. (17) and (19) we employ an iterative scheme. First, let the system be spatially discretized on a grid x = {x 0 , x 1 , . . . , x N } where the first gridpoint, x 0 , fulfills the same requirements as the lower limit of the integrals in Eqs. (17) and (19), namely that the filling function at values x < x 0 remain constant for the entire duration. Introducing the grid spacingx n = x n+1 − x n and the contraction W (y,0)→(xn,t) V O (λ) ≡ W (n) (λ), we can write Eq. (19) for the source term as follows Separating the two terms, as done in the final line, highlights the iterative nature of the equation. Starting from W (0) (λ) at the very left side of the grid, one iterates over the spacial grid while updating the source term at each iteration. Using the very same approach, the indirect propagator of Eq. (17) can be calculated. Let ∆ (y,0)→(xn,t) V O (λ) ≡ ∆ (n) (λ) and D 0 (U(x n , t; λ) ≡ D (n) (λ) for a lighter notation. Then, where we with a slight abuse of notation have introduced the term W (2) (λ) only contributes when U(x n , t; λ) > y, i.e. when any quasi-particles found at point (x n , t) originated to the right of y. By the definition of x 0 , this term should vanish for n ≤ 0. Meanwhile, the first term of (27) is dependent on all previous terms, however, only quasi-particles with rapidities within the root set λ * (x m , t; y) contribute. In some non-relativistic models there is no upper limit to the quasi-particle velocity, like the Lieb-Liniger model where the group velocity is linearly proportional to the rapidity. Since the quasi-particle move according to a purely Eulerian equation with no inhomogeneous couplings, the overall rapidity distribution of the system is conserved throughout evolution. Thus, the grid must be chosen large enough such that no quasi-particle of the initial state has a large enough rapidity to reach x 0 within the time t. Finally, the name "source term" for W (n) (λ) becomes apparent when considering Eq. (28). Given the definition of x 0 , the source term at the leftmost grid point must vanish, W (0) (λ) = 0. Thus, the solution to Eq. (28) is likewise zero ∆ (0) (λ) = 0. As we iterate over x, we eventually reach a point where quasi-particle starting at y appear, whereby the source term becomes non-zero. This is turn results in a non-vanishing indirect propagator.
In order to solve (28) numerically, we have to discretize the rapidity as well. First, we rearrange the equation as follows Next, we discretize the quantities following the notation of [65], where rapidity and type indices are lower and upper ones, respectively. In this notation the dressing operation of Eq. (6) reads h dr = U −1 h, where the matrix is U kl ij = δ ij δ kl +w l j T kl ij ϑ l j . Since all the models treated in this paper only have a single quasi-particle type, we omit the type index for readability.
Hence, denoting g (n) (λ) ≡ ρ s (x n , t; λ)f (x n , t; λ) the equation above can be written as where X (n) depends on the previous iterations. Thus, the calculation of the indirect propagator at the n'th grid point reduces to solving a linear matrix equation.

B Comparing light cones of different models
In this section we illustrate the difference in the spreading of correlations between a relativistic quantum field theory (sinh-Gordon) and the non-relativistic Lieb-Liniger model. Hence, Fig. 4 displays the correlations from a bump release in the two models. Both bumps were realized using inhomogeneous chemical potentials: For the Lieb-Liniger model, the inverse temperature is β = 0.25, the coupling is c = 1, and the chemical potential µ(x) = 2 − 2x 2 . For the sinh-Gordon model the inverse temperature is β = 0.25, and we have α = 0.0369, m = 0.9989, and µ(x) = 2 − 2x 2 .
In the Lieb-Liniger model, the quasi-particle group velocity is directly proportional to its rapidity. Thus, there is no maximum velocity, and the emerging light cone of the direct correlations has a smooth edge. Meanwhile, the group velocity in the relativistic sinh-Gordon model is bounded, as it scales a v gr (λ) ∼ tanh λ. Thus, the light cone of its direct correlations has a characteristic sharp edge.
Interestingly, the indirect correlations of the two models are fairly similar, and do not reflect the quasi-particle velocities to the same extend. Instead, the indirect correlations are mainly determined by the inhomogeneity of the system, which is fairly similar in the two cases (both are bump releases). Unlike the direct correlations, the indirect correlations do not decrease monotonically but in fact increase at first. We can understand this from the definition of the indirect correlations, namely how they are a consequence of the change in quasi-particle trajectories due to inhomogeneity. Over time, more and more particle will cross the point y = 0, thus increasing the indirect correlations. Meanwhile, as the correlations disperse, they start trailing of as ∼ t −1 . These two competing effect produces the light cones observed in Fig. 4.

C TBAs of studied models
The models used for studying the spreading of correlations are all very well known. We here briefly summarize their thermodynamic Bethe ansatz. A brilliant feature of GHD is the universality of the equations. Thus, each model needs only to provide a few specific quantities: the one-particle energy (λ), the one-particle momentum p(λ), the differential two-body scattering phase T (λ, λ ), and the statistical factor f (λ).

C.1 Lieb-Liniger model
The Lieb-Linger model describes a one-dimensional Bose gas with contact interactions governed by the Hamiltonian [17,18] whereψ † (x),ψ(x) are the bosonic fields, while c is the interaction strength, µ is the chemical potential, andh = 2m = 1. The TBA detailed here is only valid for repulsive interactions c > 0. Thus, the three main functions (single-particle energy, momentum and scattering) required for solving the GHD equations read The quasi-particles in the Lieb-Liniger TBA are fermions, whereby their statistical factor reads f (λ) = 1 − ϑ(λ) .
Like the Lieb-Liniger model, the thermodynamic Bethe ansatz is determined by only a single root density. The TBA functions of the model read where we have added a chemical potential, µ, to the energy function. Once again, the quasiparticles of the sinh-Gordon model are fermions, whereby Additionally, the expectation value of the vertex operator is given by [82,86,87] where and The one-particle form factor of the vertex operator used for calculating the correlations then reads

C.3 Classical hard rod model
The notion of integrability is not limited to quantum models but applies to some classical models as well. One of these models consists of hard rods of length a on a one dimensional line [25,26,88]. Similar to models described above, a TBA description exists for the hard rod model, where the relevant quantities are For classical models, such as the hard rod model, the statistical factor is merely The filling function for a thermal state can be easily written in terms of the source term w (th) (λ) = βλ 2 /2 as outlined in [63] ϑ (th) (λ) = e −ε (th) (λ) = e −w (th) (λ)−W (a d(β)) , with the thermal pseudoenergy ε (th) (λ) being the solution of Eq. (1) with the source term w (th) (λ) ε (th) (λ) = w (th) (λ) + W (a d(β)), while d(β) has been defined after (24) and W (a d(β)) is the Lambert-W function, which is defined as the solution of the equation Similarly ρ (th) s (λ) = 1 2π(1 + W (a d(β))) , whence, together with (4) the expression for the thermal root density ρ (th) p (λ) in (24) immediately follows. From the latter the thermal linear density distribution of rods is constructed in the Monte Carlo simulations as explained in Section 4.2 of main text and in the following paragraph.

D Details of the Monte-Carlo simulations for the hard rods gas
We present here additional details about the Monte-Carlo simulations presented in Section 4.2. In our simulations we fix the initial point −L (L > 0) whence rods are distributed in space and the number N of particles. The initial position L M of the rightmost rod is therefore fluctuating for each different realization of the initial rods' configuration. The number N of particles is chosen such that it is larger than the average number N of rods contained in the interval [−L, L] (we take it symmetric for simplicity) where we want to compute the dynamics of the density n(x, t) : The expression of the initial linear density n(x, 0) used in the simulations is given in (25). As a consequence L M > L. The simulations are performed in infinite volume, however, the initial rods' distribution is zero outside the interval [−L, L M ] and there are two depletion zones that move inwards as time elapses in proximity of which the GHD solution does not hold anymore. Velocities are drawn from a Gaussian distribution with variance 1/β(x) according to Eqs. (22) and (25). Notice that one can account for boosted thermal distributions by replacing λ → λ − µ in Eq. (25). In all the simulations presented in the manuscript we have set for simplicity µ = 0. The density and the two-point correlation function are computed by averaging over the number M of independent sampled initial conditions n(x, t)n(x, 0) c where N i (x, t) and N i (0, 0) denote the number of rods at time t in the interval (x−l/2, x+l/2) and at time 0 in (−l/2, l/2) respectively, for the i = 1, 2...M realization of the initial rods' positions and velocities. The results obtained for a double thermal bump release on top of a constant thermal background for rod length a = 0.1 and inverse temperature profile β(x) as per Eqs. (22) and (25) are further reported in Fig. 5 for completeness. The parameters are as follows: x 0 = 300, β as = 10, z = 120, β in = 0.4 and a = 0.1. The number of rods used in the Monte-Carlo simulations is N = 270, L = 660 and l = 10. For t = 15 and t = 30 we use 2 · 10 6 samples, while for t = 70 and 90, since the noise in the simulations increases, the sampling is enlarged to 7 · 10 6 and 8 · 10 6 samples respectively.
Similarly to the cases analyzed in the main text, for the density dynamics n(x, t) an excellent agreement with the GHD results is obtained for all the times shown (t = 15, 30, 70 and 90). As long as correlation functions are concerned, instead, for t = 15, 30 a deviation with the Euler scale expression is present, for longer times t = 70, 90 an excellent agreement is again attained. This aspect is witnessed by the relative distance σ between the two methods, that for t = 15 is significantly larger than the one for the values t = 30, 70 and 90. In the latter cases, σ is solely determined by the noise in the Monte-Carlo sampling