Euler-scale dynamical fluctuations in non-equilibrium interacting integrable systems

We derive an exact formula for the scaled cumulant generating function of the time-integrated current associated to an arbitrary ballistically transported conserved charge. Our results rely on the Euler-scale description of interacting, many-body, integrable models out of equilibrium given by the generalized hydrodynamics, and on the large deviation theory. Crucially, our findings extend previous studies by accounting for inhomogeneous and dynamical initial states in interacting systems. We present exact expressions for the first three cumulants of the time-integrated current. Considering the non-interacting limit of our general expression for the scaled cumulant generating function, we further show that for the partitioning protocol initial state our result coincides with previous results of the literature. Given the universality of the generalized hydrodynamics, the expression obtained for the scaled cumulant generating function is applicable to any interacting integrable model obeying the hydrodynamic equations, both classical and quantum.


Introduction
One-dimensional isolated many-body quantum systems have recently received a great amount of attention thanks to ground-breaking advances in cold-atomic experiments, see e.g., Refs. [1,2] for excellent reviews on the subject. In one dimension powerful analytical techniques are available for integrable models, which possess an extensive number of conservation laws, see e.g., the collection of reviews in Ref. [3]. In particular, it is by now well established, that integrable systems do not relax to a thermal steady state after an homogeneous quantum quench [4][5][6], where a global parameter of the Hamiltonian is abruptly changed and the system undergoes unitary evolution from a translationally invariant initial state. Relaxation occurs to the so-called generalized Gibbs ensemble (GGE), first proposed in Ref. [7] with Refs. [6,8] reviews on the subject, which accounts for the presence of the infinite number of conservation laws in the thermodynamic limit.
Out of equilibrium, the initial state of the system is, however, often inhomogeneous and non-stationary. The study of the dynamics of interacting integrable models from inhomogeneous states is much more difficult than the one in homogeneous quantum quenches. In this context, the introduction of the generalized hydrodynamics (GHD), in the independent works in Ref. [9] and Ref. [10], has been a groundbreaking advancement. This theory is an extension of hydrodynamics to interacting integrable systems: while Euler hydrodynamics describes fluids with a finite number of conservation laws, the generalized hydrodynamics is concerned with integrable systems, admitting an infinite number of conserved quantities. As a consequence, the theory of generalized hydrodynamics is based on the GGE, instead of the canonical Gibbs ensemble. As GHD is a hydrodynamic theory, at its leading order, it applies in the so-called Euler-scaling limit, where the scale of space and time at which variations in the state of the system are observed, is taken to infinity. Under this assumption, at any space-time point (x, t) the many-body system can be considered to locally relax to a GGE which depends on (x, t), in the spirit of a local-equilibrium approximation [11]. The dynamics is ruled by hydrodynamic equations having a form analogous to the Euler equations. The input of the method is the thermodynamic Bethe ansatz (TBA), see, e.g., Refs. [12,13]. Despite the name, the TBA is not restricted to quantum models: it applies also to classical systems, such as the hard-rod gas [14], the classical sinh-Gordon [15] and the Toda chain [16,17]. The initial application of the GHD formalism in Refs. [9,10], and later in Refs. [14,15,[17][18][19][20][21][22][23][24][25], has been to the study of ballistic transport from the partitioning protocol initial inhomogeneous state, but the versatility of GHD allows to study various inhomogeneous setups such as the effect of confining potentials [26,27], bump-release protocols [28,29], correlation functions [30][31][32] and entanglement spreading [33][34][35][36][37]. Remarkably, the formalism can also be used to study nonintegrable systems, provided the integrability breaking is weak enough so that on large enough length and time scales, the conserved charges may be assumed not to be broken [26,[38][39][40][41][42][43][44][45][46][47].
We also note that GHD has been experimentally verified in Ref. [48].
A complete understanding of the non-equilibrium processes where transport is present requires, however, the knowledge not only of the mean values of the ballistically transported conserved densities and the associated currents -the natural quantities at the basis of hydrodynamics -but also of the corresponding fluctuations. A relevant framework to account for fluctuations is given by the large deviation theory, see e.g., Refs. [49][50][51]. This framework builds on the computation of the large deviation function, which provides the statistics of rare but significant fluctuations. It generalizes to non-equilibrium conditions the concepts of entropy, thermodynamic phases and phase transitions, that are canonical quantities for describing the equilibrium statistical mechanics of many-body systems. The large deviation function can be computed from the associated scaled cumulant generating function (SCGF) [49], which is the non-equilibrium analogue of the free energy. In the context of isolated systems, until recently only few results for transport fluctuations were present, mostly regarding one-dimensional critical systems [52][53][54] and non-interacting models, such as the Levitov-Lesovik formula for free fermions [55][56][57][58][59][60][61][62][63][64], the free Klein-Gordon field theory [65] and harmonic chains [66]. Only recently, results became available in interacting integrable systems for the statistics of various observables not related to transport, see e.g., Refs. [67][68][69][70], and most importantly for the present paper, for ballistic transport [71,72]. In the latter, an exact expression for the current fluctuations at the Euler scale has been derived combining GHD and the large-deviation theory.
The results described above apply to homogeneous and stationary GGE states -even though these states can be obtained from inhomogeneous settings, such as the partitioning protocol for non-equilibrium steady states. At present, no general formula is available for the Euler-scale fluctuations of current flows within states affected by a nontrivial large-scale dynamics. The work [73] provided an important first step, obtaining, for non-interacting fermionic and bosonic systems, a formula for the scaled cumulant generating function describing the energy current fluctuations far from the connection point in the partitioning protocol, where a nontrivial Euler-scale dynamics occurs. The derivation is based on stationary phase methods, which are characteristic of free models, and therefore it is difficult to immediately generalize it to interacting integrable systems. A full understanding of the Euler-scale fluctuations of current flows in interacting integrable systems in states with large-scale inhomogeneity and dynamics is therefore still missing.
In this manuscript we aim at filling this gap, by providing an exact expression, in the Eulerscaling limit, for the SCGF of the time-integrated current associated to an arbitrary ballistically transported conserved charge. Our results apply to a broad class of inhomogeneous and dynamical states, which include, e.g., the initial state of the partitioning protocol, as shown in Fig. 1. The calculation of the SCGF is based on the biasing of the measure of the initial state by the exponential of the time-integrated current, in a similar way as the procedure followed in Refs. [71,72]. The biasing of the measure is shown to be fixed by the knowledge of two-point correlation functions in the inhomogeneous initial state, derived in Ref. [30]. Our derivation accordingly shows that the study of two-point correlation functions in GHD [30][31][32] is tightly related to the large-deviation theory of current fluctuations. The biasing technique developed in [72], generalised here to inhomogeneous, dynamical situations, is in principle applicable to a wide class of hydrodynamic systems, integrable or not. In this paper, however, we concentrate on integrable systems, where all the necessary technical tools are available for the technique to lead to calculable results.
We provide exact expressions for the first three cumulants of the time-integrated current and we further show that in the non-interacting limit, for the partitioning protocol initial state, our general formula for the SCGF reduces to the result of Ref. [73]. Our study of the SCGF thereby deeply generalizes the analysis of Refs. [71,72] and of Ref. [73] by accounting for inhomogeneous situations and interactions, respectively. The general expression presented for Figure 1: Schematic representation of the fluctuations of the time-integrated current in the Euler-scaling limit. In the Figure we consider, as an example, the partitioning protocol inhomogeneous state, where two GGEs, denoted as ϑ l (in red at x < 0) and ϑ r (in blue at x > 0), are joined at x = 0 at the time t = 0. The ensuing dynamics is described in terms of a continuum of states, constant on the light rays parametrized by ξ = x/t, which interpolate between ϑ l and ϑ r . We consider the time integral ∆q i * (z x, z t) = z t 0 dτ j i * (z x, τ) of the current j i * , associated to the conserved charge q i * , flowing across the point z x over the time interval (0, z t). The fluctuations of ∆q i * (z x, z t) in the Euler-scaling limit are encoded in the SCGF G(s, x, t) = lim z→∞ 1/(z t) ln 〈 exp(s ∆q i * (z x, z t))〉 inh,z . In the case depicted in the Figure, the average 〈. . . 〉 inh,z is taken over the partitioning protocol inhomogeneous state. We emphasize, however, that the formula for G(s, x, t) derived in the manuscript applies more generally to a wide class of inhomogeneous and dynamical states varying at large scales z.
The rest of the paper is organized as follows. In Sec. 2 we review the basic concepts underlying the generalized hydrodynamics theory of integrable systems that will be needed for the presentation of our results. In Sec. 3 we briefly review the large-deviation theory for the rare fluctuations of ballistic current flows. Sec. 4 contains our main result for the exact expression of the SCGF in inhogeneous and dynamical initial states together with the applications of the result to the calculation of the cumulants and to the non-interacting limit. In Sec. 5 we draw our final conclusions. Some technical aspects are reported in the Appendix A.

Introduction to the generalized hydrodynamics
In this Section we review integrable models and their generalized hydrodynamics theory. The results that will follow in Sec. 4 apply to this class of systems. Reference [26] and the lecture notes in Ref. [76] provide a very good review on the subject. In Subsec. 2.1 we recall the basic identities regarding the description of the GGE in the thermodynamic Bethe ansatz formalism.
In Subsec. 2.2, we first define the Euler-scaling limit for a generic local observable, and then we report the main GHD hydrodynamic equations. In Sec. 2.3 we briefly review the expressions of the dynamical two-point functions in inhomogeneous and non-stationary states, derived in Ref. [30] and later numerically evaluated in Ref. [37], which constitutes a central building block of the results presented in Sec. 4.

The GGE and its thermodynamic Bethe ansatz description
We consider interacting integrable models in one spatial dimension, where, in the thermodynamic limit, an infinite number of local (or quasi-local) conserved charges {Q i } exists, with i, j = 1, 2 . . . ∞. For Hamiltonian time evolution one of the conserved charges is the Hamiltonian of the system H; henceforth we further assume that the conserved charges commute among themselves [Q i , Q j ] = 0, for any i and j. The charges can be written as where the density q i associated to Q i obeys a continuity equation with the current density j i , corresponding to the current J i = dx j i (x), written in the form 1 We will indicate by q 0 , q 1 and q 2 , respectively, the number, momentum and energy charges.
In non-integrable models the latter are the only conserved quantities. The thermodynamic properties of integrable systems are described by the so-called generalized Gibbs ensemble (GGE) [6][7][8], which generalizes the canonical Gibbs ensemble to the case where an infinite number of conserved charges is present and Z GG E is the normalization constant. The parameters {β i } are Lagrange multipliers whose knowledge uniquely determines the GGE. Their values are fixed by the initial values of the conserved charges {Q i }. In this manuscript we are using the superscript notation {β i } following the notation of Refs. [71,72], but the superscript i should not be confused with an exponent. Notice that the present discussion applies both to classical and quantum models. In the former case ρ GG E is a statistical distribution in the phase-space, while in the latter it is a density matrix. The thermodynamic Bethe ansatz (TBA) formalism, see, e.g., Refs. [12,77], provides an efficient method to represent the GGE. The system is described in terms of quasi-particles with rapidity λ, where each quasi-particle carries a bare charge h i (λ) corresponding to the single-particle eigenvalue of Q i in Eq. (1). The statistics of the quasi-particles is encoded into the free energy function F ( ). The latter is given by F ( ) = −ln(1 + e − ) for fermions, F ( ) = ln(1 − e − ), for bosons, and F ( ) = −e − for classical particles. The GGE in Eq. (2) is then fixed by the pseudo-energy function (λ), which, in turn, is determined by solving the non-linear integral equation where T (λ, µ) = Θ(λ, µ)/2π and Θ(λ, µ) is the differential two-body scattering phase. The term w(λ) on the right hand side of Eq. (3) is named source term and it is defined as where the {β i } are the parameters of the GGE in Eq. (2). For a Galilean model, for example, h 0 (λ) = 1, h 1 (λ) = λ and h 2 (λ) = λ 2 /2. The only model-specic quantities in Eqs. (3) and (4) are the single-particle eigenvalues h i (λ), the integral kernel T (λ, µ), the free energy function F ( ), and the spectral space dµ. The latter, in the cases where multiple quasi-particle species n are present, has to be enlarged to n dµ. In this work we will drop the summation over the different quasi-particle species for brevity, but the results of Sec. 4 straightforwardly generalize to those cases. From Eq. (3) it is customary to introduce the mode-occupation (or filling) function ϑ(λ), the root density ρ p (λ) and the state density ρ s (λ) as where the dressing of a generic function h dr (λ) of λ is defined by the linear integral equation The GGE may be identified in various equivalent ways. The root density ρ p (λ) completely specifies it, as well as the mode occupation function ϑ(λ); since the knowledge of either determines the average values of all the conserved charges {Q i } of the model. As a consequence, in the following, we will the denote the average of a local observable O(x, t) at point x and time t over the GGE in Eq. (2) as where the second equality follows from the fact that the GGE is homogeneous and stationary. The dependence of ϑ on λ is omitted in this notation. In Eq. (7) we have used for convenience the quantum-mechanical trace notation, for classical systems it has to be replaced by an integral in phase space.

The Euler-scaling limit and the GHD equations
The TBA formalism and the GGE describe integrable systems in homogeneous and stationary states. The generalized hydrodynamics (GHD) theory accounts for situations where the state of the many-body system is inhomogeneous and non-stationary. In order to introduce this formalism, we first define the Euler-scaling limit, as the latter will be a central concept that will be used throughout the rest of the manuscript. Let us consider the case in which an integrable system is initialized in some generic state ρ 0 which is inhomogeneous and nonstationary. The average of a local operator O(x, t) over a generic inhomogeneous state ρ 0 will be denoted as 〈O(x, t)〉. Upon considering inhomogeneities which slowly vary in space-time, one can use the local relaxation assumption or "local maximization of entropy principle" (an excellent book on this is Ref. [11]), which asserts that the system locally relaxes to a GGE: The space-time dependence is recovered, in the spirit of a local-density approximation, by promoting the GGE Lagrange multipliers {β i } → {β i (x, t)}, and, equivalently, the mode occupation function ϑ → ϑ(x, t), to be dependent on the space-time point (x, t) where the operator is located. Equation (8) becomes exact in the limit of infinite length scale of the variation of the inhomogeneity. In the case of large, but finite, length scales for the variation in space and time of the inhomogeneity, instead, Eq. (8) is only approximate. The limit where Eq. (8) becomes exact is usually named "Euler-scaling limit" or simply "Euler scale". Henceforth we will use both names interchangeably. For one-point functions at the space-time point (x, t), we will use the subscript ϑ(x, t) to denote the averages over the local homogeneous GGE in Eq.
To make the discussion more concrete, here we introduce a particular class of inhomogeneous and dynamical initial states, which generalize the GGE in Eq. (2), which are given by In this expression z is a scale parameter which has to be large enough such that the resulting Lagrange parameters {β i (x/z)} are functions of the position x that vary only on large scales. If the multipliers {β i } do not depend on space, the state ρ 0 is a GGE as in Eq. (2). For this reason, we will refer henceforth to ρ 0 as inhomogeneous GGE.
In this manuscript, we will denote 〈. . . 〉 inh,z the averages over ρ 0 in Eq. (9). In the case of ρ 0 in Eq. (9) the Euler-scaling limit can be achieved by taking the limit where the scale of the inhomogeneity is sent to infinity z → ∞, simultaneously with the space-time observation point (x, t) of the operator O. In formulas A comment is in order here: in fact, sharp jumps in the state are also allowed and describable by Euler hydrodynamics. These are interpreted as "weak solutions" to the hydrodynamic equations. Therefore, β i (x, 0) do not need to be smooth functions: jumps are allowed, and the limit in Eq. (10) is still described by Euler hydrodynamics. A particular initial inhomogeneity of the form of ρ 0 in Eq. (9), which we will discuss, is the partitioning protocol state, where indeed a sharp jump occurs. In this state two homogeneous GGEs, conventionally denoted as right (r at x > 0) and left (l at x < 0) are joined at the point x = 0 at time t = 0, see Fig. 1. Accordingly, the generalized inverse temperatures are chosen as where Θ(x) is the Heaviside step function. The initial filling function ϑ(x, 0, λ) corresponding to Eq. (11) is given by where ϑ l,r (λ) are fixed by the boundary conditions one imposes on the left {β i l } and right {β i r } reservoirs. For our results in Sec. 4 we will consider initial inhomogeneous and dynamical states of the form of Eq. (9) in the Eulerian limit z → ∞.
Once the Euler-scaling limit in Eqs. (8) and (10) is assumed, it is rather simple to derive the hydrodynamic equations ruling the evolution of the system. Starting from Eq. (1), one has The mean value of the charge density in Eq. (1) is given by the TBA description where the dressing operation has been defined in Eqs. (5) and (6). In particular, all dressed quantities become functions of space and time, being determined by ϑ(x, t) via Eq. (6). The expression of the current in Eq. (1) has been first discovered in Refs. [9,10] and it reads as where h 2 (λ) = E(λ) is the bare single-particle energy eigenvalue and the effective velocity v eff is given by where in the last equality we used that for non-relativistic models h 1 (λ) = P(λ) = λ (even though the second relation in Eq. (16) holds also for more general parametrizations of the momentum as a function of the rapidity). The effective velocity v eff is a generalization of the group velocity v g , which takes into account the interactions among the quasi-particles via the dressing operation. v eff has been first defined in Ref. [78] in the context of the lightcone spreading of correlations in homogeneous quantum quenches. We stress that Eq. (15) goes beyond the results available from the TBA and it has been first proved in relativistically invariant quantum field theories in Ref.
[10] and numerically tested in quantum spin chains in Ref. [9]. Later, it has been proved in a variety of contexts including quantum spin chains and classical models [79][80][81][82][83][84][85]. Inserting Eqs. (14) and (15) into the continuity equation (13), and assuming the completeness of the set of local (or quasi-local) charges h i (λ), the final GHD equation for the filling function, see Refs. [9,10] for the original derivation, turns out to be Equation (17) express the fact that the occupation function ϑ(x, t, λ) represents the normal modes of the hydrodynamics, which, in integrable systems, form a continuum labelled by the rapidity λ. At the Euler scale, the normal modes are convectively transported with velocity v eff (x, t, λ). Generalizations of Eq. (17) to account for the presence of trapping potentials [26], diffusive corrections [43,[86][87][88][89], space-time variations of the interaction terms of the Hamiltonian [38] and Markovian coupling to an external bath [40] have been further developed. For the derivation of our results in Sec. 4 the form in Eq. (17) will be sufficient. As a final piece of introduction, we mention that Eq. (17) admits a solution by the characteristic function, U(x, t, λ, t 0 ), encoding the position at time t 0 of the characteristic curve x(U , t, λ, t 0 ) of the quasi-particle with rapidity λ that at time t is in x (therefore U(x, t 0 , λ, t 0 ) = x). The characteristic curve x(U , t, λ, t 0 ) is defined as the curve tangent to the effective velocity v eff in Eq. (16). The function U(x, t, λ, t 0 ) is defined by inverting x(U , t, λ, t 0 ) with the respect to the initial position U, see Refs. [29,36] for a detailed discussion. From Eq. (17), it is simple to see that the filling function ϑ( Remarkably, in Ref. [29], it has been proved that U(x, t, λ, t 0 ) can be determined by solving the following integral equation where x 0 is an asymptotically stationary point (see the discussion after Eq. (83) in Appendix A). Equations (18) and (19) provide an efficient way to solve the main GHD equation (17) initial value problem, as shown in Ref. [29]. The integral equation (19) for U(x, t, λ, t 0 ), furthermore, is fundamental for the derivation of the exact expression of the Euler-scale dynamical correlation functions, which will be the object of the next Subsection.

Euler-scale dynamical correlation functions
The discussion of Subsec. 2.2 focused on the Euler scaling of one-point functions (mean values) of local observables. Since for weak inhomogeneities the system is composed by locally homogeneous space-time fluid cells, one-point functions at the space-time point (x, t) can be computed based on the knowledge of their expression 〈O(0, 0)〉 ϑ(x,t) in the local homogeneous GGE in Eq. (2) at (x, t) according to Eq. (10). Euler-scale connected correlation functions are, instead, harder to compute as they depend on the whole inhomogeneous state ρ 0 characterizing the initial state of the system, and not only on the homogeneous GGE on a specific fluid cell. Considering two-point functions at the space-time points (0, 0) and (x, t), the approach for their calculation [11,32,90] consists in studying the linear response of the system at the space-time point (x, t) to a local perturbation at the point (0, 0). Correlations, at the Euler scale, are then determined by the stable normal modes propagating between the two space-time points. In this way, exact expressions for the two-point functions of generic local observables in the large-scale limit can be obtained. However, the expressions obtained in this way in Refs. [11,32,90] assume that background fluid state at the initial space-time point (0, 0) is homogeneous and stationary. In Ref. [30], this approach has been extended to account for two-point correlation functions from inhomogeneous initial states of the form in Eq. (9), under homogeneous evolution (for a dynamics with an homogeneous and time-independent Hamiltonian -that is, without spatially or temporally dependent external force fields). We report just the main ideas behind the derivation together with the final expression of the the dynamical two-point function, which will be essential for the derivation of the scaled cumulant generating function detailed in Sec. 4. Consider a small perturbation of one Lagrange parameter β j ( y, 0) → β j ( y, 0)+δβ j ( y, 0) in the initial state ρ 0 in Eq. (9). The response of the system to the small perturbation of β j ( y, 0) is related to the connected correlation function involving the associated density q j ( y, 0). Consider the average of some other density q i (x, t) over ρ 0 in Eq. (9): the functional derivative of this average with respect to β i ( y, 0) is with the Euler-scaling limit for the connected correlator defined as the subscript ϑ 0 denotes the filling function which characterizes globally, as a function of space, the inhomogeneous state ρ 0 in Eq. (9) at the time slice t = 0 in the Euler-scaling limit z → ∞. For the partitioning protocol initial state, for example, ϑ 0 is given in Eq. (12). This notation, where time appears as a lower index, stresses the fact that two-point correlation functions depend on the whole initial inhomogeneous state ϑ 0 of the system and not only on the homogeneous GGE ϑ(x, t) at a specific space-time fluid cell (x, t). For one-point functions, instead, we have denoted with 〈q i (x, t)〉 Eul ϑ 0 = 〈q i (0, 0)〉 ϑ(x,t) the average over ρ 0 in Eq. (9) in the limit z → ∞ according to Eq. (10) and in agreement with the notation introduced in Subsecs. 2.1 and 2.2. [One-point functions of the charge densities and of the associated currents are given in Eqs. (14) and (15), respectively.] We emphasize that, for some specific models, the Euler-scale limit of two-point correlation functions requires additional care and q i (z x, z t), q j (z y, 0) in the r.h.s. of Eq. (21) have to be averaged over space-time fluid-cells, as explained in Ref. [30] and shown in Ref. [15] for the classical sinh-Gordon field theory. For the classical hard-rod gas considered in Ref. [37], instead, fluid-cell averaging is not necessary and the Euler-scale predictions for the two-point function are simply obtained by averaging in space. Equation (20) represents an extension of the fluctuation-dissipation theorem of statistical mechanics [91] to the framework of GHD, where infinitely many conserved charges are present. From Eq. (20) an expression for the two-point connected correlation function can then be derived exploiting the fact that the expression of the onepoint function on the l.h.s is known in Eq. (14), as shown in Ref. [30].
These results for correlation functions of conserved densities can then be extended to the connected two-point function of any pair of local observables O(x, t)O ( y, 0) c,Eul ϑ 0 , via hydrodynamic projections [32]. The main idea of the latter is that at the Euler scale correlations are determined by the ballistically propagating conserved charges. Correlations can then be determined by "expanding" the operator O(x, t) of interest in the basis of the local (or quasi-local) conserved charges where the scalar product 〈q k |O〉 ϑ(x,t) is computed at the space-time point (x, t) of the observable as . (23) f (x, t, λ) is dubbed statistical factor of the model. For models with fermionic quasi-particle (23) is named one-particle-hole form factor of the operator O, it a functional of the filling function defined such that Eq. (23) holds. It must be worked out for every operator individually. For charge densities and the associated currents they are simply related to dressed single-particle eigenvalues h dr i in Eq. (6) [32] as (23) is the inverse of the correlation matrix (cf. Eq. (23) and Eq. (24) for V q i ) The factor C −1 ϑ(x,t) is introduced in Eq. (22) because the densities {q i } are not orthonormal under the scalar product (C ϑ(x,t) ) i j = 〈q i |q j 〉 ϑ(x,t) given in Eqs. (22) and (23). A detailed justification of the hydrodynamic projection identity in Eq. (22) is provided in Refs. [32,92].
For our purpose, it is sufficient to say that an exact formula for two-point correlations of can be derived using Eq. (22) and the result for the two-point correlator of charge densities obtained from Eq. (20). The formula reads as The square brackets in Eq. (26) denote the integral-operator notation The propagator, Γ ( y,0)→(x,t) (λ, µ), describes how the local perturbation at the space-time point ( y, 0) is transported by the normal modes of the hydrodynamics on given characteristic curve U(x, y, λ, 0), until it reaches the space-time point (x, t). In particular, Γ ( y,0)→(x,t) can be split into two parts The first term is dubbed direct propagator, it describes the normal-modes propagation along the curved characteristic curve U(x, t, λ, 0) (see Eqs. (18) and (19) in Subsec. 2.2) within the inhomogeneous fluid background. Therefore, only the quasi-particles with the suitable rapidity λ to propagate from ( y, 0) to (x, t) enter in this term. On the other hand, ∆ ( y,0)→(x,t) is named indirect propagator and it encodes a more subtle effect. ∆ ( y,0)→(x,t) describes the perturbation to the trajectory of the rapidity λ due to the interaction with other rapidities λ . Quasi-particles with arbitrary rapidities λ , not necessarily connecting the two observation points ( y, 0) and (x, t), are therefore involved into the definition of the indirect propagator. Inserting Eq. (28) into Eq. (26) an analogous decomposition for the two-point correlation formula applies where the term on the first line of the r.h.s of Eq. (29) is named direct correlator, while the term on the second line is dubbed indirect correlator. The set λ (x, t, y, 0) = {λ : U(x, t, λ, 0) = y} over which the sum on the first line runs expresses the fact that direct correlations are determined solely by the propagation of the quasi-particles with the right rapidity λ to connect the space-time point ( y, 0) with (x, t).
Note that in Eqs. (26)- (29) we have placed for convenience the operator O ( y, 0) at the space-time point ( y, 0). However, the very same formulas apply for the operator O ( y, t 0 ) at a generic time t 0 = 0 upon replacing U(x, t, λ, 0) with U(x, t, λ, t 0 ), according to Eqs. (18) and (19), and similarly for all TBA quantities evaluated at the time 0. In particular, Eqs. (26)- (29) apply also in the case t < 0, where the operator O(x, t) is computed at a time earlier than O ( y, 0), and time ordering is not needed. This is a consequence of the fact that the Euler-scale equation (17), on which Eqs. (26)- (29) are based, is time reversible. The indirect propagator ∆ ( y,0)→(x,t) can be obtained by solving a linear integral equation, which we report in the Appendix A for completeness.
Only very recently, in Ref. [37], ∆ ( y,0)→(x,t) and the formula in Eq. (29) have been numerically evaluated for various interacting integrable models (Lieb-Liniger, sinh-Gordon and the classical hard-rod gas) in different non-equilibrium protocols involving inhomogeneous initial states. Moreover, in Ref. [37], the first numerical demonstration of the validity of the Eulerscale formulas in Eq. (29) has been given by comparing them with Monte Carlo simulations of the microscopic hard-rod dynamics. An excellent agreement has been, in particular, found for long times and large values of the scale parameter z of the inhomogeneity in Eq. (9), consistently with the expectation that the Euler-scale predictions apply to weakly inhomogeneous settings.

Large deviation principle and the scaled cumulant generating function
Considering the continuity equation in Eq. (1), each conserved density satisfies, for a particular density q i * denoted by the subscript i * we define the time-integrated current q i * could denote, for example, the energy, the particle or any other density of the model. The GHD equation in Eq. (17) describes ballistic transport of the hydrodynamic modes. Ballistic transport, indeed, is relevant in integrable models as a consequence of the presence of an infinite number of conserved charges. For ballistic motion, the time-integrated current is expected to depend extensively on the time t as ∆q i * (x, t) ∝ t for large times t, and it is therefore convenient to focus on the intensive variable J i * = ∆q i * (x, t)/t. According to the large deviation theory, see, e.g., Refs. [49][50][51], the probability density function of J i * for large times t peaks exponentially around the mean, and most-probable, value 〈J i * 〉 as In Eq. (31) the notation " ", taken from Ref. [49], denotes equality of the right and the left hand side in logarithmic scale in the long-time limit t → ∞. Equation (31)  Here, as we are interested in the Euler-scaling limit, the SCGF is defined as The average is taken over the rescaled inhomogeneous state ρ 0 in Eq. (9). Equation (32) serves to define the meaning of the probability p(J i * ) in (31), which depends on the initial state, and therefore in particular on z; it is the scaling in z as per Eq. (32) that gives rise to the asymptotic equality represented by " " in Eq. (31). From the knowledge of G(λ, x, t) in Eq. (32), one can obtain by Taylor expansion the cumulants {c k } of the time-integrated current, which are defined as The validity of the large deviation principle in Eq. (31) implies that all that the connected correlation functions {〈[∆q i * (z x, z t)] k 〉 c inh,z } scale as z, with the cumulants {c k } being therefore finite and the series expansion in Eq. (32) being valid in some interval of s around the origin. The main result of this manuscript is the derivation of an exact expression for G(s, x, t) in Eq. (32) valid for interacting integrable models, both classical and quantum, which are described by GHD.

SCGF for homogeneous GGEs: review of the result
In Refs. [71,72] a general approach has been developed to compute G(s, x, t) in the case where the average in Eq. (32) is taken w.r.t. the homogeneous GGE ϑ as in Eqs. (2) and (7) (and not w.r.t. the inhomogeneous state ρ 0 in Eq. (9), as it will be the case in Sec. 4). In this case, because of the homogeneity of the GGE, G(s, x, t) = G(s). Let us review this approach.
The analysis of Refs. [71,72] relies on the GHD description of the Euler-scale hydrodynamics. In particular, we need to introduce the so-called flux jacobian matrix, with i the row index and j the column one, which will play an important role in the analysis of this Subsection. This matrix describes how the average currents depends on the average densities and therefore it depends on the equation of state of the model. In the case of integrable systems, where the averages in Eq. (35) are computed over the homogeneous GGE ϑ in Eqs. (2) and (7) as per Eqs. (14) and (15), an expression for the matrix elements of A j i can be given by means of the hydrodynamic projection formalism, as shown in Ref. [32]. In the basis of the single-particle eigenvalues with the effective velocity v eff (λ) given in Eq. (16). The function h i dr (λ) denotes the orthonormal conjugate of h dr i (λ), which satisfies the orthogonality and completeness relations respectively. Notice that the matrix A j i in integrable systems is infinite-dimensional, yet it can be formally defined. From Eqs. (36) and (37) which shows that the flux jacobian has a continuous spectrum indexed by the rapidity λ with eigenvalue the effective velocity v eff (λ).
In order to compute G(s) one then biases the GGE measure in Eq. (2) by multiplying it by the exponential of the time-integrated current appearing in Eq. (32). Averages over this tilted measure become dependent on the parameter s conjugate to the time-integrated current 〈O〉 ϑ(s) . Note that we are here extending the notation introduced in Eq. (7) for the average over the homogeneous GGE by promoting ϑ(s) to be dependent on s (in addition to the rapidity λ whose dependence is not reported for brevity). In Subsec. 4.2 we will explain this procedure in more details, for the moment we just report the main result for the flow equation from Ref. [71,72].
The flow equation describes how the homogeneous GGE ϑ in Eq. (2) is modified by the insertion of the exponential of the time-integrated current. Fundamentally, the s modified state ϑ(s) is still an homogeneous GGE, yet with Lagrange parameters {β n } dependent on s. This dependence is captured by the flow equation, which can be written as where the sign of the flux jacobian is obtained by diagonalizing the latter and by taking the sign of its eigenvalues. For interacting integrable systems, from Eq. (36), this implies By inserting Eq. (39) into Eq. (3) and by taking the derivative with respect to s, the flow equation can be recast in a form where the dependence of the pseudo-energy (λ; s) on s is exposed Notice that, as a consequence of Eq. (39), all the dressed quantities, which depend on the GGE state, acquire an additional dependence on the parameter s. Exploiting Eqs. (39), and equivalently (41), the SCGF G(s) can be eventually computed as where we stress that the current expectation on the r.h.s. is taken over the homogeneous GGE and it is thereby given by Eq. (15). In Eq. (42) we have extended the notation introduced after  41) to fix the dependence of the state on s, and then the result must be plugged into Eq. (42) to get G(s). This procedure has been carried on in Ref. [71] for the classical hard-rod gas and for the quantum Lieb-Liniger model, where G(s) has been computed in the homogeneous steady state developing at long times from the partitioning protocol initial state in Eqs. (11) and (12). In the classical hard-rod fluid, moreover, the cumulants obtained from the series expansion of G(s) according to Eq. (33) have been compared against Monte-Carlo simulations, finding an excellent agreement and thus corroborating the validity of the approach. In concluding this Subsection, we mention that the flow equation in Eq. (39) can be proved, as shown in Ref. [72], solely on the basis of linear fluctuating hydrodynamics and hydrodynamics projections [11,90] without the need of using any tool coming from integrability and Bethe ansatz. Only in deriving the form in Eq. (41) the TBA machinery is exploited. The expression for G(s) in Eq. (42) together with Eq. (39) is therefore expected to apply more generically to any system, integrable or not integrable, displaying ballistic transport.

The Euler-scale SCGF for inhomogeneous states
This Section contains the main result of this work, which is the exact expression of G(s, x, t) in the Euler-scaling limit z → ∞ as per Eq. (32), with the initial state ρ 0 the inhomogeneous GGE in Eq. (9). In Subsec. 4.1 we state for the reader's convenience the result and the equations necessary for the numerical computation of G(s, x, t). In Subsec. 4.2 we report the derivation of the results. In Subsec. 4.3 the results for the first three cumulants of the time-integrated current, according to Eq. (33), are reported. In Subsec. 4.4 the non-interacting limit of our general expression is studied.

SCGF for inhomogeneous GGEs: statement of the result
For the inhomogeneous and non-stationary GGEs ρ 0 , as in Eq. (9), the evaluation of G(s, x, t) is much more difficult and has not been considered before. In this Section we provide the first expression for the SCGF for inhomogeneous and non-stationary states in the form of ρ 0 in (9) in the limit z → ∞ according to Eq. (32).
Consider the inhomogeneous state ρ 0 in Eq. (9). Then G(s, x, t) defined in Eq. (32) can be computed as where the average current is evaluated on the homogeneous GGE ϑ(x, τ; s ) in Eq. (2) and it is therefore given by Eq. (15). The fundamental difference with respect to the homogeneous case, using the same extension of the notation introduced for Eq. (42), is that the GGE ϑ(x, τ; s ) depends not only on the parameter s as a consequence of the biasing of the measure in Eq. (9) via the exponential of the time-integrated current, but also on the (scaled) spacetime coordinates (x, τ) of the fluid cell, because, for every deformation parameter s , the state is inhomogeneous and non-stationary.
To be accurate, the bias of the measure by the exponential of the time integral of the current as in Eq. (32), depends not only on s, but also on the parameters x, t characterising the (scaled) space-time position of the integration interval (see, for instance, Fig. 1). Therefore, an average at the fluid cell (x , t ) in the deformed state should be denoted as 〈. . . 〉 ϑ(x ,t ;x,t,s) . For lightness of notation, we keep implicit the x, t dependence of the bias itself in the fluid-cell average notation; these can be considered as fixed parameters throughout. The dependence on s is important, as for instance this is integrated over in Eq. (43).
As in the homogeneous case, the s dependence is described by a flow equation for an sdependent state, which in the inhomogeneous case is however significantly more complex than Eq. (39). In terms of the fluid-cell Lagrange parameters described by β n (x , t ; s), it takes the following form: where the propagator Γ (x,τ)→(x ,t ) has been defined in Eq. (28) with τ as initial time according to the discussion after Eq. (29). In Eq. (44) we are further using the integral-operator notation introduced in Eq. (27). The one-particle form factor V j i * of the current j i * is given in Eq. (24). It is convenient to re-write this as a flow equation for the pseudo-energy . This directly follows upon differentiating the left and the right hand side of Eq. (3) with respect to s where in the first line we have used the chain rule and the relation ϑ( ) = dF ( )/d . In the second line we have again used the chain rule and the identity h dr n (λ) = ∂ /∂ β n , which follows upon differentiating with respect to β n Eq.
where we have used the completeness relation in Eq. (37). Using the decomposition of the propagator Γ in Eq. (28), the last equation can be written as where we have defined the set of times τ (x , t , λ, x) = {τ : U(x , t , λ, τ) = x}, with the characterstic curve defined in Eqs. (18) and (19). Notice that for the homogeneous GGE ϑ, the indirect propagator ∆ vanishes, as detailed after Eq. (83) in Appendix A. The characteristic curve, from Eq. (19), in this case is simply given by U(x , t , λ, τ) = x − v eff (λ)(t − τ), the set τ is composed by one element only and Eq. (47) reduces to Eq. (41). In the same limit, the time-integral in Eq. (43) trivializes, as the current average value is independent on time, and Eq. (42) is re-obtained. One therefore sees that Eqs. (44) and (47) generalize the results of Refs. [71,72], recalled in Subsec. 3.2, by including non-stationarity and inhomogeneous situations, when motion occurs at the Euler scale of hydrodynamics. Equations (43), (44) and (47) are the main results of this paper.
Equations (43) and (47) involve the propagator Γ (x,τ)→(x ,t ) which, as we have seen in Sec. 2.3, describes the motion of the quasi-particles between (x, τ) and (x , t ) in the inhomogeneous fluid background (under the homogeneous evolution Hamiltonian of the model considered). It is consequently clear that the Euler-scale expression for G(s, x, t) is expected to apply in the same limit where the expression for the two-point correlator of Subsec. 2.3 does. In particular, based on the findings of Ref. [37], where the formulas for two-point functions have been tested again numerical simulations of the microscopic hard-rod dynamics, we expect the expression for G(s, x, t) in Eqs. (43) and (47) to apply for smooth initial inhomogeneities, large enough z in Eq. (9), and long times. The actual verification of this statement by comparing the predictions coming from Eqs. (43) and (47) against simulations of the hardrod gas will be addressed in a future publication, together with the evaluation of G(s, x, t) for some specific quantum models, e.g., the Lieb-Liniger and the sinh-Gordon. In this work, we will solely present the general theory based on Eqs. (43) and (47) describing the large deviation theory of ballistically transported conserved quantities.
We conclude this Subsection by noting that, if the state ϑ 0 is invariant under simultaneous rescaling of space and time (x, t) → (a x, at), and therefore β n (a x , at ; 0) = β n (x , t ; 0) with a > 0 an arbitrary positive constant (e.g., the partitioning protocol initial state in Eqs. (11) and (12)), then the expression in Eq. (43) becomes a function of the scaling variable ξ = x/t. As a matter of fact, using the property of the propagator Γ , valid if β n (a x , at ; 0) = β n (x , t ; 0) for every n, x and t , Γ (x,τ)→(x ,t ) = a Γ (ax,aτ)→(ax ,at ) , which directly follows from Eqs. (28) and (81), the following scaling property is obtained for the Lagrange multipliers β n (x , t ; s) for arbitrary values of s by using Eq. (48) in Eq. (44) β n (a x , at ; s) = β n (x , t ; s) .
Equation (49) in turn implies that the same scaling property holds also for the pseudo-energy (x , t , λ; s) and any dressed quantity h dr i (x , t , λ; s): Exploiting the last identity, the expression for G(s, x, t) can be eventually written in the form In the next Subsection, we provide the derivation of the main results Eqs. (43), (44) and (47).

SCGF for inhomogeneous GGEs: derivation of the result
We start by defining the s-tilted ensemble as that obtained by biasing the measure of the inhomogeneous GGE in Eq. (9) with the time-integrated current in Eq. (30). This procedure is analogous to the one introduced in Refs. [52][53][54]93] for one-dimensional critical systems described by conformal field theory, and later extended in Refs. [62,73] to non-interacting systems. We will further comment in Subsec. 4.4 about the relation between the present approach and how it simplifies to the case of non-interacting systems.
Averages over the s-tilted ensemble will be denoted as 〈.
where ∆q i * is given in Eq. (30). From the definition (52) one also has with the integrand on the r.h.s. denoting the two-point connected correlation function over the s-tilted ensemble with the notation of Eq. (52). It is crucial to note that, although the s-tilting involves an integral over time, this does not affect the dynamics. The bias is to be understood as a modification of the measure, that is, of the distribution of states at time 0, or on any chosen time slice. According to (52), the measure is modified by a weight which is evaluated by evolving the observables in time (or equivalenty, evolving the distribution of states in time), and evaluating the exponential of the time-integrated current. The dynamics is still given by the original, homogeneous and time-independent Hamiltonian of the model under consideration, and thus, at the Euler scale, time slices are related to each other by the original GHD equation (17) of the model. This is important, as we can then use, below, the results for correlations in inhomogeneous and dynamical states reviewed in Subsec. 2.3, which, as emphasised, are based on the assumption of an homogeneous and time-independent dynamics.
The definition of the s-tilted measure in Eq. (53) is widely known in the context of large deviations in classical stochastic systems, see, e.g., Refs. [94,95], and in open quantum systems. In the latter framework it is named "s ensemble", see, e.g., Refs. [96,97]. There, the approach consists in relating the bias in s, as an exponential of the time-integrated current, to a modification of the Lindbladian ruling the evolution of the system. This operation can be done via the so-called "quantum Doob" transformation, as shown in Refs. [98,99], and it amounts to a biasing of the probability measures which makes rare events typical. The approach we pursue here is complementary to the quantum Doob transform, in the sense that here we relate the insertion of the exponential of the time-integrated current to a change of the (inhomogeneous) GGE measure. Fundamentally, the latter modification still produces a (inhomogeneous) GGE, whose s-dependence can be determined exactly in terms of the flow equation in Eqs. (44) and (47).
We now show that the biasing procedure in Eqs. (52) and (53) defines a flow in the manifold of the inhomogeneous GGEs. Since, by construction, the dynamics is unchanged, we may think of the manifold of inhomogenenous GGEs as described by a space-dependent GGE on any given time slice, {β n (x , t ; s)} x ,n . Let us choose the time slice t = 0. One can adapt the argument developed in Ref. [72] for the homogeneous case, summarized in Subsec. 3.2. Namely, one defines the Lie derivative L at the point 〈...〉 Eul ϑ 0 in the manifold of inhomogeneous GGEs as Using Eq. (26) with τ as initial time for the propagator Γ (x,τ)→(x ,t ) , and Eqs. (23) and (24), the Lie derivative can be written as Equation (55) shows that the Lie derivative lies on the tangent space to the manifold of inhomogeneous GGEs identified by the Lagrange multipliers {β n (x , 0)} x ,n . Comparing Eqs. (54) and (55) with (53), and remembering the definition of the Euler-scaling limit of two-point correlation functions in Eq. (21), one realizes that the s-tilted measure defines a flow directed along the Lie derivative as z → ∞. Accordingly, from Eq. (53), infinitesimal s modifications lie on the tangent space to that manifold, and given that at s = 0 the state 〈. . . 〉 (0) inh,z is an inhomogeneous GGE and lies on that manifold, then it remains so even after the s-tilting. As a consequence, one has for the Euler-scaling limit z → ∞ of Eq. (52) that where in the second equality we have used the local relaxation assumption (8) thereby expressing the average of the local operator O over the local homogeneous GGE in Eq. (2) at the fluid cell (x , t ; s). It is also worth to remark that in the second equality of (56) we have extended the notation introduced after (21) for the mode-occupation function ϑ 0 (s), corresponding to the state ρ 0 in Eq. (9) as z → ∞, by exposing the additional dependence on s due to the biasing of the measure. In order to fix the additional dependence on s due to the tilting in Eq. (52) one further needs to specify the flow equation. To do this we start from (53) together with Eq. (56) by choosing the local observable O(x , t ) as some conserved density q n (x , t ) of the model The homogeneous GGE at every fluid cell (x , t ) is specified by the Lagrange multipliers {β n (x , t ; s)} x ,n or, equivalently, by the average conserved densities 〈q n (x , t )〉 Eul ϑ 0 (s) = 〈q n (0, 0)〉 ϑ(x ,t ;s) in Eq. (14). Equation (57) can therefore be considered as the "equation of motion" of the state coordinates 〈q n (0, 0)〉 ϑ(x ,t ;s) for their trajectory, parametrized by s, in the manifold of inhomogeneous GGEs. For convenience, we express these coordinates on an arbitrary time slice t .
As a matter of fact, Eq. (57) relates the tangent vector of the trajectory, the l.h.s., to a function of the coordinate itself, the r.h.s.. Since the expression of the correlator on the r.h.s. is known, from Eq. (29) in Subsec. 2.3, the equation of motion is well defined and it fixes the flow of the inhomogeneous GGE in Eq. (9) as a function of s. Indeed we now show that Eq. (57) is equivalent to the flow equation (44). Exploiting the chain rule one has where in the last step we used vector-matrix notation and the symmetry of the correlation matrix C ϑ(x ,t ) , defined in Eq. (25). Notice that in Eq. (58) we have further exploited the fluctuation dissipation relation in Eq. (20) where the second equality comes from the fact that for equal-time connected correlation func- . The latter equality simply expresses that at the Euler scale fluid cells on the same time slice separated in space are uncorrelated, as shown in Ref. [30]. The inverse matrix C −1 ϑ(x ,t ;s) can be defined by means of the functions h j dr in Eq. (37) as Once the flow equation is established, proving Eq. (43) for the SCGF is simple. By applying a derivative with respect to s to G(s, x, t) in Eq. (32), with ρ 0 as given by Eq. (9), one obtains where in the first equality we have used the definition in Eq. (52) for the average over the stilted ensemble, while Eq. (56) has been exploited in the second and third equality. Integrating in s, and exploiting the fact that G(s = 0, x, t) = 0, from the definition in Eq. (32) of the SCGF, one can see that Eq. (61) reduces to Eq. (43).
Regarding the numerical evaluation of Eqs. (43) and (47), we remark that they are expressed in terms of quantities available from the TBA for s = 0, i.e., in the unbiased case, such as v eff (x, δ, λ; 0), h dr i * (x, δ, λ; 0) and V j i * (x, τ; 0). These functions can be efficiently evolved in time using the solution of the GHD equation in Eq. (17) via the characteristic curves in Eqs. (18) and (19), as discussed in Ref. [29]. The most difficult element to numerically compute, which, instead, is not directly provided by the TBA, is the indirect propagator ∆ (x,τ)→(x ,t ) . For the calculation of the latter, however, an efficient numerical scheme has been developed in Ref. [37] within the iFluid open-source code [100]. Accordingly, the only step which has not yet been pursued in the literature, so far, is the numerical solution of the flow equation in Eq. (47). This passage can be achieved by writing Eq. (47) as an integral equation (by formally integrating in s the right and the left hand side). The latter can be in turn solved by iteration starting from the knowledge of the right hand side of Eq. (47) for s = 0, where all the TBA functions and the propagator ∆ (x,τ)→(x ,t ) are known as just explained. The resulting expression for the pseudo-energy (x , t , λ; s) as a function of s can then be plugged into Eq. (43), together with the expression in Eq. (15) for the average current 〈 j i (0, 0)〉 ϑ(x,τ;s ) , to eventually compute G (s, x, t). We plan to carry on this numerical analysis in a future work.
In concluding this Subsection, we emphasize that the present derivation of G(s, x, t) is strongly based on the result of Ref. [30], recalled in Subsec. 2.3, for two-point functions in inhomogeneous and dynamical GGEs. Accordingly, Equations (43), (44) and (47) are restricted to integrable models, differently from the results in Eqs. (39) and (42) for homogeneous GGEs. The calculation of the scaled cumulant generating for inhomogeneous and dynamical states in generic, not necessarily integrable, models is an unexplored and challenging problem which goes beyond the analysis of the present manuscript.

Analysis of the cumulants
We report here the first three cumulants, which can be obtained from the series expansion in Eqs. (32) and (33) of the general expression given by Eqs. (43) and (47). The first cumulants provide important information about the shape of the probability distribution of J i * = ∆q i * /t and they are the easiest to access experimentally. Moreover, the cumulants can be computed numerically, e.g., with Monte-Carlo simulations in classical systems such as the hard-rod gas, and therefore they can be used to test the predictions of the Euler-scale formula in Eqs. (43) and (47) with simulations of the microscopic dynamics. Higher cumulants can be in principle derived as well, even if the derivation becomes combinatorially more cumbersome as the order increases.
The first cumulant is trivial from Eq. (43) and it is just the time integral of the GHD expression of the current in Eq. (15) For higher-order cumulants we need the useful relation, with X dr (x , t , λ; s) a generic dressed quantity (we drop for brevity all the independent variables), which directly follows upon differentiating with respect to s Eq. (6). ∂ /∂ s is specified by the flow equation (46) and (47). The second cumulant is related to the Gaussianity of the distribution close to the mean value and, from Eqs. (43) and (63), its expression reads which, using (46) turns out to be that is in agreement with the expression for the connected current-current two-point correlation function one obtains from Eq. (29), therefore providing an important consistency check of Eqs. (43) and (46). The expression of the third cumulant c 3 is related to the asymmetry or skewness of the distribution. Its expression was not known before and it is given by We notice that both Eq. (64) for c 2 and Eq. (66) for c 3 are written in terms of ∂ /∂ s, which is given by the flow equation (46) and (47). Once the latter is numerically solved, also the cumulants can be therefore evaluated numerically. We note that in the cases where the initial state ϑ 0 is invariant under the rescaling (x, t) → (a x, at), with a > 0 an arbitrary constant, and therefore Eqs. become scaling functions of ξ once they are rescaled by t. We have numerically checked, by simulating the hard-rod gas dynamics from an initial step inhomogeneity in the inverse temperature β 2 (x, 0) as in Eqs. (11) and (12) and in Fig. 1, that the first three cumulants c k (ξ), with k ≤ 3, of the particle flow (h 0 (λ) = 1) are functions of ξ. In this case, we have observed a linear growth as a function of t of the connected correlation functions, 〈[∆q 0 (ξt, t)] k 〉 c inh,z with ξ fixed and k ≤ 3, which implies that the cumulants are finite (the same scaling behavior is expected also for higher cumulants with k > 3) and the large deviation principle applies, as commented after Eq. (34). In the homogeneous and stationary limit of the results in Eqs. (43), (44) and (47), the same scaling behavior as a function of t for the cumulants was already observed in Ref. [71]. The fineteness of the cumulants in integrable models implies that there is no divergence in the derivatives of G(s, x, t) w.r.t. s and therefore no dynamical phase transition in the time-integrated current J i * statistics, see, e.g., Refs. [101][102][103]. The understanding of the precise reason behind the absence of divergences in the s-derivatives of G(s, x, t) in integrable systems is still under investigation and it will not be addressed in this manuscript.

The non-interacting limit
In this Subsection we show how to specify the general result in Eqs. (43) and (44) to noninteracting systems. In particular, we show that the flow equation for the Lagrange multipliers simply reduces to a shift as a function of s, as already noted in the homogeneous case in Refs. [71,72]. This property is characteristic of non-interacting systems and can be ascribed to the fact that the Euler equations in Eq. (13) are in this case linear.
The fundamental simplification happening in non-interacting systems is that no dressing is present since the scattering phase Θ(λ, µ), or equivalently T (λ, µ) in Eq. (3), vanishes identically. The flow equation in Eq. (44) therefore simplifies drastically. Specifically, the propagator Γ (x,τ)→(x ,t ) (λ, µ) in Eq. (28) for free models is uniquely given by the homogeneous contribution (see also Appendix A) where v eff (x, t, λ) = v g (λ) = dE(λ)/dλ is the group velocity since there is no dressing. We observe that for the calculation of G(s, x, t) in Eq. (43), where the mean current 〈 j i * (0, 0)〉 ϑ(x,τ;s ) is evaluated on fluid cells at the fixed space point x and times τ ∈ (0, t) (see Fig. 1), one has to specialize the propagator to equal space points x = x, i.e., Γ (x,τ)→(x,t ) (λ, µ). Note that in interacting systems the indirect propagator ∆ (x,τ)→(x ,t ) (λ, µ) vanishes as a consequence of the fact that the quasi-particles trajectories do not depend on the state (see also Appendix A).
with the last step following from (36) without the dressing. For the same reason A j i is independent of s and Eq. (69) can be trivially integrated In order to further simplify Eq. (69) we briefly recall some basic identities about the Lagrange mulipliers {β n }. As already commented after Eq. (57), the GGE state at every fluid cell (x , t ) can be equivalently described in terms of the Lagrange parameters {β n (x , t )} x ,n or with the densities 〈q n (0, 0)〉 ϑ(x ,t ) . As a matter of fact, it simple to see from Eq. (13), using the definition of the matrix C in Eq. (25), that the multipliers satisfy an hydrodynamic equation similar to the one the densities fulfill, as shown in details in Refs. [9,10], Since A j i in the absence of dressing is independent of space and time, Eq. (71) is linear. Accordingly, one can introduce the normal mode function N (x , t , λ) via a simple linear combination of the {β n (x , t )} x ,n , i.e., where h n are the bare single-particle eigenvalues, see Eq.
Note that Eq. (73) generalizes to the interacting case upon replacing v g (k) → v eff (x , t , λ). As a matter of fact, one can notice that Eq. (73) has the very same structure as Eq. (17). Indeed, the function ϑ(x , t , λ) is the normal mode function associated to the conserved densities 〈q n (0, 0)〉 ϑ(x ,t ,λ) , as already commented after Eq. (17) Let us now specialize the discussion to the partitioning protocol initial inhomogeneous state in Eqs. (11) and (12) with an initial thermal inhomogeneity in β 2 (x , 0) and all the other Lagrange parameters initially set to zero (see, for example, Fig. 1). For non-interacting models this protocol has been addressed, e.g., in Refs. [62,63,65,73,[104][105][106][107][108][109]. In this case, the solution of Eq. (73) is readily found to be with ξ = x /t , while N l,r (λ) are fixed by the initial conditions for the left and right half of the system. Then, from Eq. (72) computed at time t = 0, one has where we identified Note that i * = 2 in Eq. (74) and h 2 (λ) = E(λ) since we are now considering energy transport. Inserting Eqs. (75) and (77) into Eq. (74) and diving both sides by h 2 (λ) one eventually has where we defined β(ξ , λ) = N (ξ , λ)/h 2 (λ). We observe that for the partitioning protocol initial state β(ξ , λ) depends on x and t through the scaling variable ξ = x /t , and therefore the Lagrange multipliers are invariant upon simultaneous rescaling of space and time β n (a x , at ) = β n (x , t ) for every n, x and t . Accordingly, Eq. (49) is satisfied, as one can explicitly see from Eq. (78), and G(s, x, t) = G(s, ξ), with ξ = x/t according to the discussion before Eq. (51). Upon using Eq. (78) into Eq. (51) one obtains with for the filling function ϑ similarly to Eq. (75), and β r,l (s) = β r,l − s sgn(v g (λ)). For noninteracting systems, accordingly, the biasing of the measure necessary to compute the SCGF can be simply obtained by performing a linear shift of the inverse temperatures β l,r characterizing the partitioning protocol initial condition in Eq. (76). The latter statement in the case of the homogeneous stationary state ensuing from the partitioning protocol initial state, obtained upon setting ξ = 0 in Eqs. (79) and (80), was already known as "extended fluctuation relation", first proposed in Ref. [93]. The extended fluctuation relation has been shown to be valid in conformal field theory [53,54] and in free-particle models [62,63,93]. In Ref. [72], in particular, the validity of the extended fluctuation relation has been shown to be a characteristic feature of non-interacting systems, where the statistical properties of the time-integrated current are directly determined by the fluctuations of the initial state, induced by the shift of the Lagrange parameters β l,r , since the quasi-particle trajectories do not depend on the surrounding state where the quasi-particles propagate. Within this perspective, the result in Eqs. (79) and (80) constitute a generalization as a function of ξ to inhomogeneous and dynamical GGEs in Eq. (9) of the extended fluctuation relation, in agreement with the results of Ref. [73] (cf., Eqs. (84) and (85) therein). In Ref. [73], indeed, Eqs. (79) and (80) have been first derived using stationary phase methods for the partitioning protocol initial state with a thermal inhomogeneity in β 2 (x, 0), as in Eq. (76). The derivation of Eq. (74) presented in this Subsection is, instead, based on Euler-scale hydrodynamics and applies more generally to any initial state ρ 0 having the form in Eq. (9). Only in Eqs. (75)-(80) the result has been specialized to the partitioning protocol state with a step profile of β 2 (x, 0). More importantly, the Euler-scale hydrodynamics allows for the exact treatment of interactions, as we have seen in Secs. 4.1 and 4.2. In the case of interacting systems, however, quasi-particles trajectories, due to the nonlinearity of the GHD equation in Eq. (17), are modified by the dynamical inhomogeneous background within which they move. Accordingly, fluctuations do not depend directly on the initial-state fluctuations, but also on the whole motion of the quasi-particles between time 0 and t through the inhomogeneous fluid background, as expressed by the time integral in Eqs. (46) and (47).

Conclusion
In this manuscript we have addressed the problem of computing the Euler-scaling limit of the scaled cumulant generating function G(s, x, t), defined in Eq. (32) and depicted in Fig. 1, of the time-integrated currents associated to ballistically transported conserved densities in integrable models. The SCGF is of paramount importance since it gives access via the largedeviation principle in Eq. (31) to the full probability distribution p(J i * ) of the rare-atypical fluctuation of the rescaled time-integrated current J i * = ∆q i * (x, t)/t far from its mean value 〈J i * 〉. Moreover, its series expansion provides all the cumulants of the time-integrated current according to Eq. (33). We have extended the analysis of Refs. [71,72], which we have reviewed in Subsec. 3.2, for the calculation of the SCGF in homogeneous and stationary GGEs in Eq. (2) to the more complex and interesting case of inhomogeneous and dynamical GGE states of the form in Eq. (9), where we have defined the Euler-scaling limit of G(s, x, t) in Eq. (32). Our findings are detailed in Sec. 4. Specifically, in Subsec. 4.1 we have stated our main results given by Eqs. (43), (46) and (47), whose proof is provided in Subsec. 4.2. Equation (43) expresses the SCGF as an integral of the mean value of the current, given by the generalized hydrodynamics formula in Eq. (15). The crucial point is that the mean value is computed over an s-tilted inhomogeneous measure, see Eqs. (52) and (53), according the exponential of the time-integrated current. All quantities which depend on the inhomogeneous fluid state acquire, as a consequence, an additional dependence on the parameter s, coupled to the time integral of the current. This s-dependence can be determined by solving the flow equation, cf., Eqs. (46) and (47), which describes exactly how the initial measure is affected by the tilting procedure in terms of a flow, parametrized by s, in the manifold of inhomogeneous GGEs. The flow equation is based on the knowledge of the Euler-scale two-point correlation functions, whose expressions have been first derived in Ref. [30] and later numerically tested in Ref. [37] for the hard-rod gas. In Subsec. 4.3 we have provided in Eqs. (62), (64)-(66) the expressions of the first three cumulants of the time-integrated current. We have further numerically checked for the hard-rod fluid, in an initial inhomogeneous state with an inverse temperature β 2 (x, 0) profile as in Eqs. (11) and (12) and in Fig. 1, that the first three cumulants c k of the particle current, with k ≤ 3, are finite (higher cumulants are similarly expected to be finite). The numerical check of the finiteness of the cumulants is methodologically fundamental since it confirms, a posteriori, that the time integrals in Eqs. (62), (64)-(66) for the cumulants are convergent, that the series expansion in Eqs. (32) and (33) is well defined and that therefore the large deviation principle in Eq. (31) applies. In Subsec. 4.4 we have analyzed the non-interacting limit of the general formulas in Eqs. (43) and (44). In this case, due to the linearity of the hydrodynamic equations, the flow equation drastically simplifies and, in the case of the partitioning protocol initial state in Eq. (11), it reduces to a shift of the Lagrange parameters characterizing the initial state, see Eq. (78). We have further checked the the results obtained in this way reproduce those obtained in Ref. [73] via stationary-phase methods for non-interacting systems. In the latter case, the fluctuations of the time-integrated current are directly determined by fluctuations of the initial state since quasi-particles propagate along straight line characteristics with a velocity which is independent on the state. Due to the inhomogeneity and the interactions, on the contrary, quasi-particles move along curved trajectories whose shape depends on the surrounding dynamical fluid state. Consequently, the fluctuations of the time-integrated current depend on the whole trajectories of the quasiparticles, as expressed by Eqs. (46) and (47).
As a future perspective, we plan to numerically evaluate the general expressions in Eqs. (43), (46) and (47) for specific classical and quantum interacting integrable models, such as the Lieb-Liniger, the hard rods, and the sinh Gordon. The comparison between the expressions in Eqs. (62), (64)- (66) and the numerical simulations of the hard-rod model will be also carried out. From the analytical side it would be very interesting to extend our results in Sec. 4 beyond the Euler-scaling limit of infinite variation lengths of the inhomogeneities. This can be achieved by accounting for diffusive corrections to the hydrodynamic equation in Eq. (17), as done in Refs. [43,[86][87][88][89]. In the specific case of classical models, such as the hard-rod gas where diffusion has been investigated in detail in Refs. [11,14,75], this should somehow connect to the macroscopic fluctuation theory, see Refs. [110][111][112], which is also based on diffusive hydrodynamics. We surely plan to carry out this study in the future.
Funding information G.P. thanks the A. Della Riccia Foundation (Florence, Italy)-INFN-for financial support and King's College (London, United Kingdom) for hospitality in the period when this work has been carried on.

A Indirect propagator for two-point correlation functions
In this Appendix we report for completeness the necessary formulas for the evaluation of the indirect propagator ∆ ( y,0)→(x,t) defined in Subsec. 2.3, which enters in the calculation of G(s, x, t) in Eqs. (43) and (47).
The indirect propagator, as shown in Ref. [30], satisfies the integral equation with V O defined in Eq. (24) and the characteristic curve in Eqs. (18) and (19). The function D is usually named "effective acceleration". The latter has been first defined in Ref. [26] to account for weakly inhomogeneous terms in the Hamiltonian ruling the time evolution of the system, e.g., due to the presence of a confining potential. In the case of Eq. (81), in a complementary way, it expresses the inhomogeneity of the initial state in Eq. (9), while the dynamics is assumed to be dictated by an homogeneous and time-independent Hamiltonian (as commented at the beginning of Subsec. 2.3). The expression of the effective acceleration is given by [26] D(x, λ) = ∂ x ϑ(x, 0, λ) The function W ( y,0)→(x,t) is dubbed source term and it is given by where the star-dressing operation is defined as h * dr (λ) = h dr (λ) − h(λ), Θ(x) is the Heaviside step function, the set λ (x, t, y, 0) has been defined after Eq. (29), and x 0 , defined in Eq. (19), has to be fixed such that ϑ(x, t , λ) = ϑ(x, 0, λ) for x < x 0 and t ∈ [0, t], see Refs. [29,30]. In principle one should think that x 0 = −∞, although for numerical calculations, see Ref. [37], one fixes x 0 to a value sufficiently far from the point y and any other inhomogeneity. Note that the indirect propagator ∆ ( y,0)→(x,t) vanishes identically in the two limiting cases of an homogeneous initial state, i.e., ϑ(x, 0, λ) = ϑ(λ) and therefore D = 0 identically in Eq. (82), and for non-interacting systems, as already commented in Subsec. 4.4 in Eq. (68). In the case of non-interacting models the vanishing of ∆ ( y,0)→(x,t) in Eq. (81) follows from the fact that h * dr (λ) = h dr (λ) − h(λ) = h(λ) − h(λ) = 0 for any arbitrary function h(λ) since T = 0 and there is no dressing. In the two cases of homogeneous initial states and non-interacting systems correlations are then solely given by the direct term in Eqs. (28) and (29).