Generalized-Hydrodynamic approach to Inhomogeneous Quenches: Correlations, Entanglement and Quantum Effects

We give a pedagogical introduction to the Generalized Hydrodynamic approach to inhomogeneous quenches in integrable many-body quantum systems. We review recent applications of the theory, focusing in particular on two classes of problems: bipartitioning protocols and trap quenches, which represent two prototypical examples of broken translational symmetry in either the system initial state or post-quench Hamiltonian. We report on exact results that have been obtained for generic time-dependent correlation functions and entanglement evolution, and discuss in detail the range of applicability of the theory. Finally, we present some open questions and suggest perspectives on possible future directions.

We give a pedagogical introduction to the Generalized Hydrodynamic approach to inhomogeneous quenches in integrable many-body quantum systems. We review recent applications of the theory, focusing in particular on two classes of problems: bipartitioning protocols and trap quenches, which represent two prototypical examples of broken translational symmetry in either the system initial state or post-quench Hamiltonian. We report on exact results that have been obtained for generic time-dependent correlation functions and entanglement evolution, and discuss in detail the range of applicability of the theory. Finally, we present some open questions and suggest perspectives on possible future directions. GHD approach is valid at the hydrodynamic scale, but it can be extended to include genuine "quantum" effects, such as the spreading of entanglement and quantum correlations. Furthermore, differently from Luttinger-liquid hydrodynamics, it is not restricted to low energies. Since its introduction, GHD has experienced a rapid development, already partially covered in a review [40] and a set of lecture notes [41], proving to be a very versatile and powerful tool with many applications (most of which are discussed in other contributions to this volume), and being eventually experimentally verified [42,43] (see in particular the review by Bouchoule and Dubail in this volume). The aim of this manuscript, which sets it apart from the other aforementioned surveys, is to provide a self-contained pedagogical presentation of the theory, with a special focus on the study of inhomogeneous quantum quenches. In particular, we will consider two prototypical quantum quenches where translational symmetry is broken by either the initial state ("bipartitioning protocols") or the post-quench Hamiltonian ("trap quenches"), and review recent results in these two important types of problems. We point out that particular types of these quenches are also treatable by other means (see, e.g., Refs. [44][45][46]). Moreover, there exist inhomogeneous quantum quenches after which observables exhibit behaviours qualitatively similar to the ones considered here but that cannot be addressed using generalised hydrodynamics [47][48][49][50][51]. Both these aspects will not be covered in this review.

CONTENTS
The rest of this article is organised as follows. We begin in Sec. II by introducing the general formalism of GHD. There we present a pedagogical discussion focusing on a simple non-interacting model, which we use to highlight the main conceptual and formal points. We first review the case of a homogeneous quench (Sec. II B) and show how one can build on its solution to treat the simplest case of inhomogeneous quench dynamics: the bipartitioning protocol, i.e., the sudden junction of different homogeneous states (Sec. II C). More general inhomogeneous quenches are discussed in Sec. II D, while in Sec. II E we discuss the changes occurring in the presence of interactions. We continue with Sec. III, which contains a review of several works where GHD has been applied to the study of bipartitioning protocols, focusing in particular on the main physical implications of the results. Next, Sec. IV is devoted to the dynamics of quantum entanglement, while Sec. V focuses on recent developments to capture quantum fluctuations around a inhomogeneous background (namely, to compute generic time-dependent correlation functions for a particular class of states). Sec. VI presents a discussion about genuinely quantum effects in inhomogeneous quenches and how GHD can be extended to account for them. Finally, Sec. VII contains our conclusions and a brief discussion of some of the open questions.

II. GHD APPROACH TO INHOMOGENEOUS QUENCHES
In this section we provide a self-contained introduction to the main concepts of GHD, seen as a method to characterise the late-time non-equilibrium dynamics of inhomogeneous integrable systems. Since the latter relies heavily upon the formalism developed to treat the homogeneous case, we begin by considering homogeneous quenches and briefly review some of the aspects that are directly relevant for the inhomogeneous generalisation (for a comprehensive review we refer the reader to Ref. [52]). Before doing that, however, we introduce a simple system of non-interacting fermions which we will use as a paradigm. In Secs. II B-II D we use this system to exemplify the main features of GHD. The appropriate modifications of the treatment to account for interacting integrable interactions are discussed in Sec. II E.

A. A simple integrable model
Let us consider a system of spinless fermions on the lattice described by the following pairing Hamiltonian, with coupling J and chemical potential h Here we set the lattice spacing a to one, we assumed the number of sites L to be even, we denoted by c † x , c x a set of fermionic creation and annihilation operators fulfilling the canonical anti-commutation relations and we assumed periodic boundary conditions (c † L+x = c † x ). Note that in the above equations we set = 1: we will abide by this convention throughout the entire review except for Sec. VI, where the dependence will be restored for pedagogical reasons.
The Hamiltonian (1) is mapped onto a particular spin-1/2 chain (the celebrated Transverse Field Ising Model) through a Jordan-Wigner transformation [53]. The mapping from fermions to spins, however, is non-local and forces one to distinguish between two sectors that differ in the boundary conditions. Specifically, periodic boundary conditions for the fermions are mapped into periodic or anti-periodic boundary conditions for the spins and vice versa. This requires a slightly more refined analysis but does not change the main results that we are going to present, which are independent of the boundary conditions. To avoid this inessential complication we stick to the fermionic form (1).
Since the Hamiltonian (1) is a quadratic form of fermionic operators and in addition is translational invariant, it is directly diagonalised by a combination of Fourier and Bogoliubov transformations. Namely, it can be written as Here the sum runs over a discrete, or "quantized", set of rational multiples of π and we introduced the "dispersion relation" and the "Bogoliubov fermions" Here is commonly referred to as the Bogoliubov angle. Note that b † k , b p fulfil the canonical commutation relations (2) with x and y replaced by k and p. We also stress that the quantization conditions for the quasi-momenta appearing in the sum (3), i.e., e iLkj = 1 (8) do not couple different momenta. From the representation (3) we see that the Hamiltonian is non-interacting: it is written as an independent sum of mode operators describing the occupation of a single Bogoliubov mode (a mode k has energy ε(k)). This means that the eigenstates of the Hamiltonian are readily written as where |Ω is the vacuum state such that b k |Ω = 0 for all k.
Let us now proceed to illustrate three key concepts that will play an important role in the upcoming discussion.

Conserved charges
The Hamiltonian (1) commutes with the following set of operators Q n = k q n (k) n k − 1 2 n = 0, 1, 2, . . . , L − 1, where q n (k) are called "bare charges" or "single-particle eigenvalues" and are defined as q 2n (k) = 2Jε(k) cos(nk), q 2n+1 (k) = 2J sin((n + 1)k) , n = 0, 1, 2, . . . . (12) One can directly see that Q n commute with the Hamiltonian because they are linear combinations of the mode operators (note in particular that Q 0 = H). These conserved operators, or "charges", have the following three key properties: (i) They can all be written as sums of densities that are local in space, i.e., act non-trivially only on a finite number of neighbouring sites. Specifically, we have [33,54] Q n = L x=1 q n,x + const n = 0, 1, 2, . . . , L − 1, with [55,56] (14) Conserved operators with this property are called local conservation laws (or local charges).
(ii) Their number is extensive, i.e., it is proportional to the volume L.
(iii) They are linearly independent.
These are general features of integrable models. In fact, possessing an extensive set of independent local conserved charges can be considered the defining feature of an integrable system, although, in some cases, the notion of independence can be non-obvious, see, e.g., Ref. [57].
Let us now consider the commutator between q n,x and the Hamiltonian. Since q n,x is local and H has a local density, the commutator is also local. Moreover its sum over the entire chain gives [H, Q n ] = 0. Therefore, there must exist a local operator j n,x such that [H, q n,x ] = i(j n,x − j n,x−1 ), (16) where we included i in order for j n,x to be Hermitian. The local operator j n,x is called current operator associated with the charge Q n and is completely specified by (16) and the condition Ω|j n,x |Ω = 0. Note that (16) is nothing but a continuity equation and can be brought in the standard form by adopting the Heisenberg picture ∂ t q n,x (t) = i[H, q n,x (t)] = j n,x−1 (t) − j n,x (t) , where we used the standard definition of Heisenberg operators q n,x (t) = e iHt q n,x e −iHt . As evident from (16), the current does not simply depend on the charge, but on the very definition of its density. What is probably less intuitive is that even the total current J n = x j n,x , which is the sum over the entire chain of the current operator, depends on the definition of the charge density. In particular, our definitions of charge densities (14), (15) result in the following total currents [58,59] We stress that the above considerations apply to all integrable models regardless of interactions. In particular, since there is always an extensive number of local conservation laws, we have an extensive number of operatorial continuity equations (17). The main difference encountered in the interacting case is that one cannot find explicit expressions like (14)- (15) and (18)- (19) for generic charge densities and currents.

Thermodynamic description of expectation values
Another key property of integrable models is that the expectation values in their eigenstates are expressed in terms of the "momentum distribution" of the corresponding set of quasiparticles. To understand what this means in our simple example let us look at the following expectation value where the sum in the first term of the r.h.s. runs over the momenta of the eigenstate |{k j } N j=1 . Considering now the thermodynamic limit lim th : at fixed we find dk e −ik cos 2 (Θ h (k)/2)ρ(k) + e ik sin 2 (Θ h (k)/2))(1 − ρ(k)) , where ρ(k) is the momentum distribution of the quasiparticles, also known as "root density" in the literature of integrable models, which appears as the weight of the integrals replacing the sums in the thermodynamic limit Note that some authors define the root density as the limit of a sequence ρ(k) = lim th 1 L(k j+1 − k j ) k ∈ (k j , k j+1 ) .
From the definition of ρ(k), up to O(L 0 ), we have 2πLρ(k)∆k = # of momenta in [k, k + ∆k] present in the state |{k j } N j=1 .
Eq. (22) proves that, in the thermodynamic limit, the expectation value of c † x+ c x in the eigenstate |{k j } N j=1 depends only on the state's root density: no other property of the eigenstate needs to be retained. A completely analogous reasoning can be applied to the expectation value of all other quadratic combinations of the fermionic operators in |{k j } N j=1 (provided that their distance is finite in the limit). Moreover, using that the eigenstates fulfil Wick's Theorem for the fermions c † x , c x , this statement is extended to all local operators. Importantly, the correspondence between eigenstates and root densities is not one-to-one. In fact, exponentially many (in the volume) eigenstates of the Hamiltonian correspond to the same root density. This can be intuitively understood by noting that small changes to the distribution {k j } N j=1 do not change the root density. More precisely by counting all such ineffective changes one finds that the number of eigenstates corresponding to a given [60], where we introduced the so-called Yang-Yang entropy and is the total density, namely the density of "vacancies" (or slots) that can or not be occupied by the momenta [28,29].

Root density from the charges
The last key property that we want exemplify here is the so-called "string-charge duality" [37]. Namely, the fact that there is a one-to-one correspondence between the root density introduced in the previous section and the expectation values of all charge densities in the thermodynamic limit.
One direction is straightforward: since the charge densities are local operators, the thermodynamic limit of their expectation values can be expressed in terms of the root density. In particular, we have Note that (28) has a very simple kinetic theoretical interpretation in terms of the quasiparticles. To find the density of the n-th conserved charge one considers the contribution of a particle with momentum k (q n (k)) times the number of particles of momentum in [k, k + dk] divided by L (ρ(k)dk) and sums over all allowed values of k.
The key point here is that (28) can be inverted. Namely, one can use it to find ρ(k) in terms of { q n ρ } n=0,1,... . This is a straightforward consequence of the fact that {q n (k)} n=0,1,... is a complete (but not orthogonal) set of functions in L 2 ([−π, π]). Explicitly we have the following Fourier series expression for the root density Note that the string-charge duality implies that the root density can be written as the expectation of an operator written as an (infinite) linear combination of charge densities. In the non-interacting case one can re-sum the series and obtain with and g µ (x) is a smooth approximation of the periodic Dirac delta function, i.e.

B. Stationary states after homogeneous quenches
Having introduced our toy model (1), we can now move to consider homogeneous quench problems. Specifically, let us focus on the following protocol (i) Prepare the system in its ground state |Ψ 0 = |GS(h 0 ) for a given value h 0 of the chemical potential; (ii) At time t = 0 (suddenly) change the chemical potential to h = h 0 ; This means that for t > 0 the system is no longer in equilibrium and undergoes non-trivial evolution. A crucial point for the following discussion is that both |Ψ 0 and H(h) (and hence also |Ψ(t) ) are translationally invariant.
A natural question is whether, after the quench, the system can eventually go back to an equilibrium state, i.e., whether it can relax. It is easy to understand that this cannot happen for the system as a whole. Indeed, relaxation is naturally associated with loss of information while the evolution of the system is purely unitary and retains all the information (probabilities are conserved). Nevertheless, relaxation can still be observed considering finite subsystems in the thermodynamic limit because the system can act as an effective bath on its own finite parts.
In order to make this intuition quantitative, we probe the physics of local subsystems by looking at the expectation values of local observables. More precisely, given a local operator O x , we consider where we denoted by |n the eigenstates of the Hamiltonian, with associated energy eigenvalues E n . If the limit (33) exists for any local observable O x , then we say that the system relaxes locally [52]. In this case one can find a stationary stateρ s of H (i.e. [H,ρ s ] = 0) such that for every local operator O x . In fact, condition (34) does not specify a state uniquely and can be fulfilled by many different stationary states of H. However, all states fulfilling (34) are equivalent when reduced to finite subsystems in the thermodynamic limit. Therefore they can be regarded as different representations of the same stationary state, very much like the different ensembles of standard statistical mechanics [52].
For the quench problem under examination local relaxation can be explicitly proven. Indeed, one can find a simple integral representation for lim th Ψ 0 |O 0 (t)|Ψ 0 and evaluate the limit of infinite times [61,62]. This can be done every time that both the Hamiltonian before (H(h 0 )) and after (H(h)) the quench are quadratic in the same variables. Conversely, this is generically impossible in the presence of interactions. Nevertheless, assuming local relaxation, one can generically find a representation of the stationary state without solving the dynamics. This can be done in the following two steps H1 Assume that (34) holds for some stationary stateρ s . H2 Fixρ s by imposing the conservation of the expectation value of all charge densities. Indeed, since |Ψ 0 is translational invariant, taking the expectation value of (17) we have which implies Let us now show that H1 and H2 are indeed sufficient to fix the stationary state reached after the quench (i)-(ii). We begin by following Refs. [34,35] and representingρ s microcanonically. Namely we writeρ s = |Ψ s Ψ s | where |Ψ s is a judiciously-chosen eigenstate of the Hamiltonian. The expectation values in the thermodynamic limit are then fully characterised by a root density ρ s (k) and from (36) we have π −π dk q n (k)ρ s (k) = lim th Ψ 0 |q n,0 |Ψ 0 , n = 0, 1, 2, . . . , which is readily inverted using (29). To find an explicit solution we then only have to evaluate the "initial data" lim th Ψ 0 |q n,0 |Ψ 0 . In our case this is easily done by expressing |Ψ 0 in terms of the Bogoliubov fermions of the Hamiltonian H(h) (see, e.g., Ref. [52]) and explicitly evaluating the expectation values of (14) and (15). The result reads where (cf. (7)) In fact, in our case we do not need to re-sum (29). Noting that (38) takes the same form as the l.h.s. of (37) and recalling that {q n (k)} is a complete set we immediately conclude One can explicitly verify that ρ s (k) agrees with the result found by explicit solution of the dynamics [61,62]. We point out that, even if we demonstrated it in a particular example, the above procedure to obtain the stationary state from the assumptions H1 and H2 relies on the decoupled form of the charges (11), which are written as sums over independent modes. Therefore, it can be directly applied to all non-interacting systems (its generalisation to interacting integrable systems is described in Sec. II E). From a physical point of view, the main message of the above discussion is that, after a given transient (which depends on the details of the initial states, and the Hamiltonian parameters), finite subsystems approach a stationary state that can be determined without solving the dynamics. Note this stationary state is generically non-thermal: see, e.g., the comparison in Fig. 1 between (40) and the root density of a thermal state with the same energy density. Importantly, translational symmetry (of post-quench Hamiltonian and initial state) implies that the stationary state is the same for all local subsystems independent of their spatial position. In the next section we will see how this description has to be modified when translational invariance is broken explicitly.

C. Stationarity along rays after bipartitioning protocols
Let us now look at a different quench problem. We consider the case in which the system is initially separated in two parts, "left" and "right", each prepared in the ground state of the Hamiltonian for different values of the chemical potential, say h L and h R = h L . Then, at time t = 0, the two parts are joined together and left to evolve unitarily according to the Hamiltonian H(h), with h = h R = h L (see Fig. 2). More precisely, we consider the following initial state where |GS L (h) and |GS R (h) are respectively the ground states of and q (h) 0,x is the energy density operator with field h (cf. (14)). Inhomogeneous quench problems of this kind, i.e. where the initial state is composed by the junction of two different homogeneous pieces, have been extensively studied in the literature. Of particular interest have been the sudden junction of two half chains prepared at different temperatures [38,39, and at different chemical potentials (or fillings) [38,44,50,, but more general initial states have also been considered [38,[111][112][113][114][115][116][117], particularly in relation to studies on the entanglement growth (see Section IV for references). Here we will refer to inhomogeneous quenches of this kind as "bipartitioning protocols".
The time-evolving Hamiltonian (1) is still invariant under translations but, crucially, the initial state is not. This means that the expectation values lim th Ψ(t)|O x |Ψ(t) maintain a dependence on x and we cannot directly use the strategy of the previous section to find their large-time limit. As we will see, however, the above strategy can be successfully modified at the expense of "diagonalising" the continuity equations for a system with infinitely many conservation laws.
Let us first consider point H1 . Since the expectation values retain a non-trivial dependence on the position we cannot assume (34). At the same time it is still reasonable to expect finite subsystems to reach local equilibrium Pictorial representation of the bipartitioning protocol. After the sudden junction of two homogeneous halves a growing region around the junction is affected by the inhomogeneity. This region is contained in a light cone spreading from the junction at the maximal speed attainable in the system, i.e., the speed of the fastest quasiparticles.
at large enough times. To conciliate these observations, one might think to impose a condition like (34) with an x-dependent stationary state. This idea, however, is too naive: an inhomogeneity in the stationary state inevitably produces dynamics, which makes the infinite-time limit ill defined. In fact, for the specific bipartite geometry of the initial state, true stationarity can still emerge along light cones or"rays", namely for observables moving away from the junction at fixed speed. This can be established by analytic calculations in simple cases (see, e.g., [104,[118][119][120]) and by direct numerical calculations. Formally, we have where we introduced the scaling limit lim sc,ζ : This ballistic scaling can be intuitively understood by recalling that in integrable models the dynamics is interpreted in terms of moving quasiparticles. In this interpretation the ray dependence comes naturally by noting that observables on a given ray ζ receive a blend of quasiparticles from the two edges that is fixed by ζ. The ray-dependent stationary stateρ s (ζ), known as locally quasi-stationary state (LQSS), has been introduced in [118]. The modification of point H2 is more direct. Instead of imposing the conservation of the expectation value of all charge densities, we require them to fulfil the continuity equation which is just the expectation value of (17). Putting all together, we arrive at the following strategy to predict the late-time properties of the system: B1 Assume that, in the scaling limit lim sc,ζ , every local subsystem is asymptotically described byρ s (ζ). Namely, assume that (43) holds for every local observable O x .
To show that this strategy is able to determine the stationary state we proceed as in the homogeneous case. We consider the scaling limit, plug (43) in (45), and represent the stationary state microcanonically obtaining −ζ∂ ζ q n,0 ρ ζ + ∂ ζ j n,0 ρ ζ = 0, where · ρ denotes the thermodynamic limit of the expectation values in an eigenstate with root density ρ(k), while ρ ζ (k) is the root density associated with a given ray ζ.
To proceed, we need to express the expectation values in (46) in terms of the root densities. The expression for the charge densities is reported in (28) while that for the currents reads [58] j n,0 ρ = dk q n (k)ε (k)ρ(k). Note that, as for the charge densities (28), also the expectation values of the currents can be interpreted in a kinetic theory fashion. Indeed, viewing the group velocity ε (k) as the classical velocity of the quasiparticles with momentum k, we see that (47) is the expression for the flux of charge Q n generated by the motion of the quasiparticles.
Putting all together and using the completeness of the bare charges {q n (k)} we then find The boundary condition for this equation can be found by noting that, since there is a finite speed for the propagation of signals [121], observables infinitely far from the junction relax as if the system were homogeneous. Namely whereρ s,L andρ s,R are respectively the stationary states reached after the homogeneous quenches h L → h and h R → h. Using the result (40), Eqs. (48), (49), and (50) are solved by where Θ(x) is the step function. Once again, one can directly verify that ρ ζ (k) agrees with the result found via explicit solution of the dynamics. Eq. (51) encodes complete information about the local properties of the system at large times after the quench. In particular, due to the fact that root densities fully specify the expectation values of local observables, it allows one to access their profiles throughout the whole light cone ζ ∈ (ζ min , ζ max ), where ζ min , ζ max respectively correspond to the minimum and maximum velocities of the quasiparticles. We note that, despite the discontinuous step function in (51), since the position of the step changes smoothly with k, the profiles are continuous in ζ, see, e.g., the example reported in Fig. 3.

D. Quasistationary states after general inhomogeneous quenches
Let us now look at a more general inhomogeneous quench. We consider the case in which the system is prepared in the ground state of a Hamiltonian of the form (1), but with a chemical potential h x depending non-trivially on the position. At t = 0 the chemical potential is then changed to a homogeneous value h x = h for all x and the system is left to evolve unitarily.
In this case there is no scaling limit in which the expectation value becomes exactly stationary. Intuitively, however, it is natural to expect that a form of quasi-stationarity emerges asymptotically in time. Namely, one expects that (43) can be turned into a statement of the form where denotes the leading contribution for large times, i.e., much larger than the time scale τ th of local relaxation. In writing Eq. (52) one assumes that at large enough times, and at the leading order in time, the state becomes locally stationary and homogeneous. Therefore, it can be replaced by a space-time dependent stationary stateρ s (x, t). In order for this assumption to be consistent, the stateρ s (x, t) must be slowly varying, i.e., there must exist a volume element × τ of the space-time around (x, t) such that where a and τ th are respectively the lattice spacing of the chain (1) (which we restored for the sake of clarity) and the local relaxation time, while L and L τ are the length and time scales for the variation ofρ s (x, t), cf. Fig. 4. The condition on length scales in (53) is known as local density approximation (LDA) [6]. Note that, for systems on the continuum, the lattice spacing a is replaced by the averaged interparticle distance d. Moreover, we remark that the above condition allows for the substitution (52) only for expectation values of local observables, i.e., observables with support of the order of d.
The assumption (52) directly leads to an equation for the root density. Indeed, repeating the steps of the previous section, with (43) replaced by (52), we find where ρ x,t (k) is the root density representing microcanonically the stateρ s (x, t). In the second step we used that the equation is non-trivial only when ∂ t ρ x,t (k) and ∂ x ρ x,t (k) are of the same order in time and therefore higher derivatives go beyond the accuracy of (52). We see that the final result (54) is a simple non-collisional Bolzmann equation for ρ x,t (k), interpreted as the distribution function of non-interacting classical particles moving with velocity ε (k). For a given initial condition imposed at time t = t 0 (large enough to be in the asymptotic regime) Eq. (54) determines the root density for all rescaled times t ≥ t 0 giving a quantitative description of all the local properties of the system at large times after the quench.
The discussion above is heuristic (for example, it assumes a finite relaxation time scale, τ th , while this quantity typically diverges for integrable systems). However, it can be made more precise introducing an appropriate scaling limit. A convenient way to proceed is to introduce a length scale Λ for the variation of h x , namely with g(x) smooth function, and the rescaled variables where v M = max k ε (k) is the maximal velocity of the quasiparticles. In this language, the analogue of (43) is obtained by taking the weak inhomogeneity limit Λ → ∞ (see, e.g., Appendix B of Ref. [122]) and leads to the following equation for the root density in rescaled variables where, once again, ρ ξ,τ (k) is the root density representing microcanonically the stateρ s (ξ, τ ). This more rigorous point of view becomes necessary to properly account for subleading corrections to (54), which originate from higher orders in the gradient expansion, see Refs. [122,123]. For a more thorough discussion of this point we refer the reader to Sec. VI.

E. Generalization to interacting integrable models
In the previous sections, we have explained the main ideas underlying GHD in non-interacting theories. Here we discuss how the conceptual structure carries over to the interacting integrable case. In particular, here we considered the so-called Bethe ansatz integrable models [28]. A key feature of these systems is that their spectral properties can always be understood in terms of stable quasiparticles, which display many analogies with the excitations of non-interacting theories. In particular, stable quasiparticles are parametrised by quasi-momenta, or rapidities, λ j . However, due to the interactions, the latter cannot be quantized independently from one another: for an N -quasiparticle state, the simple relations (8) are typically replaced by [28,29] e iLk(λr) = F r (λ 1 . . . , λ N ) , where F r depend on the specific model considered, while k(λ r ) is the physical momentum associated with rapidity λ r . The quantization conditions (61) are customarily called "Bethe equations". Although the structure of eigenstates in the presence of interactions becomes significantly more complicated, expectation values of conserved quantities are expressed as simple sums over {λ j }. For instance, the momentum and energy associated with a given eigenstate are where the single-particle momentum k(λ) and energy ε(λ) are model-dependent functions.
In general λ j need not to be real, and can take arbitrary values in the complex plane. In addition, there might also be different species of quasiparticles connected to different physical degrees of freedom. For instance, in a Bethe-ansatz integrable system of particles with spin (like the Hubbard model [124] or the Yang-Gaudin model [28]) one has two distinct species of quasiparticles parametrised by disjoint sets {λ j } and {µ j }. Roughly speaking one is connected with spin and the other with charge degrees of freedom. For simplicity, in this section we will restrict to the case of a single species.
For large volumes, the rapidities {λ j } arrange themselves in regular patterns in the complex plane which are obtained combining a number of "basic" configurations where the rapidities stay at a fixed distance from one another. The latter can be specified by a real rapidity and are interpreted as bound states formed by the elementary quasiparticles. Moreover, in the thermodynamic limit the values that these real rapidities can take become dense and one can describe a solution of (61) in terms of sets of quasi-momenta distributions ρ n (λ). Here, n takes discrete positive integer values, and ρ n (λ) is interpreted as the distribution of the quasi-momenta for a bound-state of n quasiparticles.
In summary, in all Bethe-ansatz integrable models eigenstates can be described using root densities as discussed in Sec. II A 2 but, differently from the non-interacting case, for each eigenstate there is now a set of functions {ρ n (λ)} (rather than a single one) with λ ∈ [−Λ, Λ] and n = 1, . . . , N b . The maximal values Λ and N b that λ and n can take depend on the details of the model and they both can be infinite.
Importantly, due to the non-trivial quantization conditions (61), the "available" values of λ that could be occupied are not distributed uniformly. Accordingly, differently from the non-interacting case (cf. (27)), the distribution ρ t,n (λ) of vacancies becomes non-trivial. The precise relation between ρ n (λ) and ρ t,n (λ) is found from (61) (see e.g. [28,29]) and reads as where k n (λ) is the momentum of an n-particle bound-state, while T n,m (λ − µ) encodes all information about the interactions (it is proportional to the logarithm of the scattering phase shift [28,29]). Note that Eq. (63) does not specify uniquely the functions ρ n (λ), and stationary states must be determined by an independent equation, which is typically expressed in terms of the function For example, in the case of thermal stationary states such additional equation takes the form where ε n (λ) is the energy of a bound-state of n quasiparticles. In analogy with the non-interacting case, the root densities completely specify the thermodynamic properties of the system. For instance, given a local charge Q n , the corresponding expectation value on the state described by {ρ n (λ)} is simply This formula has once again an intuitive kinetic theoretical interpretation, generalizing Eq. (28) to the interacting case. In fact, a generalization exists also for the expectation value of local currents, which reads where the velocities {v n (λ)} are "dressed" by the interactions as described by the following integral equation where {ρ t,n (λ)} are given in (63). Differently from Eq. (66), which follows immediately from the definition of the root densities, formula (67) is highly non-trivial. It was first conjectured in Ref. [38,39] but its rigorous proof has been accomplished only very recently [125] (see, however, Refs. [126][127][128] for relevant partial results). For more detail on this aspect we refer the reader to the contributions by Cubero, Yoshimura, and Spohn; and Borsi, Pristyák, and Pozsgay to this special issue.

Homogeneous quenches in interacting integrable models
Despite the conceptual framework being completely analogous to the non-interacting case, the study of quantum quenches in interacting integrable systems is significantly more complicated on the technical level. In fact, explicit results are typically restricted to simple families of initial states. In essence, this is due to two main complications.
(i) Strictly local conservation laws are in general not enough to uniquely specify the root densities, and one also needs to consider quasi-local conserved operators [129,130]. These charges are again expressed as sums of densities (cf. (13)) but the densities do not have finite support and exhibit exponentially decaying tails. It is customary to denote the combined set of local and quasi-local charges by {Q n,s } s=1/2,1,3/2,...;n=0,1,... where Q n,1/2 are the usual local charges. The expectation value of any charge Q n,s on a stationary state described by {ρ m (λ)} can again be written in the form (66) for some appropriate functions {q n,s,m (λ)}. Importantly, the set of all local and quasi-local charges is complete, in the sense that (ii) The analogue of (38), i.e. an explicit expression for initial-state expectation values, is in general not available. Currently, this can be found only for low entangled states [131].
When these two complications can be overcome, i.e. a complete family of quasi-local conservation laws is known and the initial-state expectation values can be computed exactly, one can follow the "string-charge duality" logic of Sec. II B and arrive at a full description of the steady state in terms of the root densities ρ n (λ). Note that in the interacting integrable case there is typically no exact solution for the full dynamics to compare with and one should test the assumption (34) against numerical results. All the cases considered so far confirmed the validity of (34), see e.g. Refs. [129,[131][132][133][134]. Finally, we should stress that the string-charge duality is not the only available approach to determine the postquench stationary state. Indeed, two complementary methods are given by the so-called "Quench Action" [35,36] and "Quantum-Transfer Matrix" approaches [135][136][137][138]. All these methods yield the same results but, depending on the specific model and initial state, some of them might be difficult (or even impossible) to implement. For further details on the latter methods we refer the reader to the relevant literature, see e.g. [36,37,135].

Inhomogeneous quenches in interacting integrable models
As for the case of homogeneous quenches, the GHD logic outlined in Secs. II C-II D can be directly extended to the interacting integrable case. Specifically, using Eqs. (66) and (67) for the expectation values of quasi-local charge densities and related currents one finds the following continuity equation for the root densities {ρ n,x,t (λ)} of the locally quasi-stationary state In particular, in the case of bipartitioning protocols we can explicitly take the scaling limit lim sc,ζ and obtain These equations differ from their non-interacting counterparts (54) and (48) because of the presence of space-time dependent velocities. This is a direct consequence of the interactions (cf. (68)) -in fact, it is the only interaction effect in (70) and (71) -and has a very intuitive explanation. The motion of a given quasiparticle is perturbed by the scatterings with the others: this results in a change in its averaged velocity. Naturally, the change depends on {ρ n,x,t (λ)}, the set of densities of the different species of quasiparticles at the space-time point (x, t), yielding space-time dependent velocities. Assuming that, for any λ and n, the equation ζ = v n,ζ (λ) has a unique solution, Eq. (71) can be immediately solved in terms of η n (λ) as follows [38] where Θ(x) is the step function, while the "left" and "right" functions η n,L (λ) and η n,R (λ) are those characterising the state at infinite distance from the junction on the right and on the left hand side, respectively. Note that (72) is still an implicit solution because v n,ζ (λ) depends on η n,ζ (λ). The standard procedure to treat it is by iteration, see Sec. III. Interestingly, the interaction-related complications outlined above do not arise when deriving Eqs. (70) and (71). Indeed, in the derivation one only needs to assume that a complete set of quasi-local charges exists, without ever needing their explicit form. The aforementioned problems, however, emerge when trying to find appropriate boundary conditions. For example, let us consider Eq. (71). Repeating the arguments of Sec. II C we have that a unique solution is found by imposing the boundary conditions (49) and (50). These are nothing but homogeneous quench problems and, as discussed before, in the interacting case they can be solved only for special classes of states. In particular, the simplest cases to treat are bipartitioning protocols where one joins two different stationary states with known root densities (in this case the homogeneous quench problems associated with boundary conditions are trivial). This includes the highly studied cases of two half chains prepared at different temperatures or at different filling. In the case of Eq. (70) the problem of finding initial condition proved itself to be even harder. Up to now it has been solved only when the system is initialised in a slowly varying stationary state, see, e.g., Refs. [139][140][141][142][143][144]. In this case ρ n,x,0 (λ) is determined in terms of the equilibrium root densities within a local density approximation.

III. LOCAL PHYSICS OF INHOMOGENEOUS QUENCHES
In this section, we provide a survey of some recent applications of the GHD approach to the study of inhomogeneous quenches in interacting integrable lattice systems. We concentrate the discussion to the case of bipartitioning protocols and only focus on the scaling limit (44) where observables become functions of the ray ζ. Moreover -although we will also mention results in other models -we will predominantly focus on the paradigmatic case of the XXZ Heisenberg chain (cf. (73)). We present the main physical results obtained in this regime and their qualitative interpretation, outlining connections with findings by alternative methods. Relevant results have also been obtained in other settings, models and limits but they are not covered in detail here. We refer the reader to the other contributions to this special issue (in particular see those by Bastianello, De Luca, and Vasseur; Bulchandani, Gopalakrishnan, and Ilievski; De Nardis, Doyon, Medenjak, and Panfil; Bouchoule and Dubail) for other examples of applications of GHD to the dynamics of local observables after inhomogeneous quenches.
We begin by exemplifying the procedure discussed in Sec. II E 2 for the case of a bipartitioning protocol in the paradigmatic case of the XXZ Heisenberg chain. The latter describes a system of spins on the lattice that interact as described by the following Hamiltonian Here {σ α,x } α=1,...,3 act as Pauli matrices on the local space C 2 at site x and like the identity elsewhere, while h is a magnetic field. The Hamiltonian (73) is related to a chain of spinless fermions with a quartic interaction term via a Jordan-Wigner transformation (we again neglect issues arising from the boundary conditions), with ∆ = 0 corresponding to the non-interacting point. The quasiparticle content of the model depends on the anisotropy parameter ∆, with a particularly simple structure observed for ∆ > 1. In that case, one has an infinite number of bound states, N b = ∞, and λ ∈ [−Λ, Λ] = [−π/2, π/2]. The driving terms and kernels of Eqs (63), (65), and (68) are given by ε n (λ) = −π sinh(η)a n (λ) + 2hn, where we set η = arcosh ∆. Let us consider a bipartitioning protocol, i.e., we initialise the system in the state and evolve it with Hamiltonian (73). The GHD prescription to compute the profiles of local charges and currents can be summarised in the following steps: I. Determine the left/right thermal stationary states (and hence the corresponding η L,R (λ)) by solving the homogeneous quench problem. For example in the case where the left and right halves of the system are initialised in two different thermal statesρ one can obtain η L,R (λ) from Eq. (65). n,ζ (λ) using Eqs (63) and (68), and obtains a new ansatz η (1) n,ζ (λ) using the right-hand side of (72). One proceeds in this way to obtain subsequent approximations η (k) n,ζ (λ), until convergence is reached. III. For each ray ζ, use the knowledge of η n,ζ (λ) to compute ρ n,ζ (λ) and v n,ζ (λ) from Eqs (63) and (68), and finally obtain the values of the charges and currents from Eqs (66) and (67). As discussed in Sec. II A 2, the knowledge of ρ n,ζ (λ) allows for the computation of all local properties of the system beyond the density of conserved charges and their currents. In practice, however, explicit formulae expressing local observables in terms of root densities are scarce. Important examples have been found for simple few-point operators in the Heisenberg spin chain [145,146], in the Lieb-Liniger model [147][148][149][150][151][152], and in the sinh-Gordon field theory [153][154][155].
These steps are very simple to implement numerically and with straightforward modifications they can be applied to any integrable model treatable with the formalism of Sec. II E. This has been explicitly demonstrated in multiple studies of bipartitioning protocols in concrete models [38, 39, 63, 83-85, 87, 89, 105, 106, 156-161], see also Ref. [162] for a versatile, open-source numerical framework for solving typical equations appearing within GHD. Furthermore, while the above prescription typically allows one to access the values of the profiles numerically, there exist cases where fully analytic solutions can be obtained, as we discuss in the following.

A. NESS
Arguably, the most interesting aspect of bipartitioning protocols in integrable systems is that they allow for the realisation of non-equilibrium steady states (NESS)s -i.e. steady states supporting non-trivial currents -in the context of isolated quantum lattice systems (i.e. without resorting to external driving). This contrasts with what happens in generic (non-integrable) lattice systems where the only local conservation law is the Hamiltonian. In the latter case assumption (43), together with some physical requirements on the form of H (for example invariance under space inversion or time reversal), allows one to prove that the current vanishes at late times (see, e.g., Sec. IX B of Ref. [40]). The latter fact is in agreement with expectations coming from the Fourier law [163,164] -which predicts a current proportional to the temperature gradient -and the currently available numerical evidence [69,77,165]. The emergence of non-trivial NESSs after bipartitioning protocols in integrable systems was first pointed out in the non-interacting case [97] and then in conformal field theories [67]. Proving that the NESS survives (integrable) interactions has been the first stark success of GHD [38,39]. In particular, in the context of lattice systems, the emergence of a NESS was first demonstrated for the XXZ Heisenberg model in Refs. [38,105], where GHD allowed for a detailed study of the dependence of charges and currents on the interaction parameter.
In the language of the previous section, the NESS is the LQSS associated with the ray ζ = 0 (cf. (43)). Namely, it is the state that captures the late-time properties of any finite region at infinite times after the quench. This state has been extensively investigated in non-interacting systems [65, 66, 68, 70-73, 75, 97, 100, 104, 114, 119, 166-170] and conformal field theories [67,[171][172][173][174][175][176][177][178] (see also the dedicated reviews [175] and [179]). An important result of these studies is the determination of all higher cumulants of the NESS currents, which give access to the full counting statistics of the charges transferred through the junction (see also [180,181]).
Within GHD, thermodynamic properties of the NESS can be studied directly from the general theory introduced in the previous section. For a given bipartitioning protocol, the value of the currents are generically found to be nonmonotonic functions of the interaction parameter, see the example reported in Fig. 5. Note that it is not uncommon to see NESS currents growing with the interaction strength. Another interesting property of the NESS currents in integrable systems is that they cannot generically be written as sums of functions involving properties of a single lead only, i.e.
This is in contrast to what happens in free systems and CFTs, and can be viewed as a transparent signature of the interaction. The simplification discussed above happens only in some special cases, see e.g. the low-temperature regime for thermal reservoirs discussed in Sec. III C. While the GHD equations can be typically solved only numerically, closed analytic expressions for the profiles may be obtained in special cases. In the Heisenberg chain this happens trivially for ∆ = 0, where the system becomes non-interacting. A more interesting example was found in Ref. [98], which considered a quench from a "domain-wall" state, i.e., a bipartitioning protocol where the left and right halves of the system are initialised in completely polarized states, in opposite z-directions. While, in this case, transport at the hydrodynamic scale is trivial for ∆ > 1, all local observables display non-trivial ballistic profiles in the regime ∆ < 1, and fully analytic expressions may be obtained. Arguably, the most interesting result of Ref. [98] is that the NESS has a nowhere differentiable dependence on the strength of interactions. In particular, the magnetisation density and current profiles exhibit jumps in correspondence to values of the anisotropy ∆ for which arcos(∆)/π is a rational number: this is in agreement with the nowhere differentiable spin Drude weight computed in the linear response regime [182][183][184][185] (for further details see [40] and the contribution by De Nardis, Doyon, Medenjak, and Panfil to this special issue). The analytic solution of Ref. [98] also allowed the authors to analyze the behaviour around the edge of the magnetisation profile, ruling out the presence of a Tracy-Widom scaling, typical of non-interacting behaviour (see Sec. VI A).

B. Phenomenology of the profiles
Although the structure of the GHD equations is completely general, the details of the underlying microscopic model encoded in in the parameters N b and Λ and functions k n (λ) and T n,m (λ) (cf. Sec. II E) give rise to a manifold phenomenology, with significant qualitative differences in distinct integrable systems. Examples of bipartitioning protocols have been studied, and often checked against independent numerical methods, in spin chains and lattices [38,63,83,89,105,106,[156][157][158]160], quantum gases [39,85,161], hard-rod systems [186], classical [187][188][189] and quantum [39,87] field theories.
In general, from the structure of the GHD equations, it is easy to see that different bound-states of quasiparticles give rise to non-analyticities {ζ ± n } n inside of the light cone and at its boundaries. These can be understood in terms of the quasiparticles' motion, as the non-analytic points ζ ± n correspond to the maximum and minimum velocities of the n-quasiparticle bound-states. The precise nature of such non-analyticities depends on the initial state and the model considered.
Once again, this could be already appreciated from GHD studies in the Heisenberg chain [38,98,105]. In the gapless regime of the model, ∆ < 1, one has a finite number of bound states for rational values of arcos(∆)/π, so that a finite number of non-analyticities appear [38]. On the contrary, the number of bound-states is always infinite for ∆ ≥ 1, giving rise to a series of non-analytic points, which accumulate inside the light cone [105]. Typically, ζ ± n are easily visible for small values of n, see Figs. 5-6 for concrete examples. As n → ∞ the velocities converge to a λ-independent value, i.e.
Accordingly, also the sequence ζ + n converges to a ray ζ ∞ which corresponds to the velocity of the heaviest quasiparticle.
Interestingly, it was shown in Ref. [105] that, depending on the initial state, the profiles of magnetisation and charges that are odd under spin-flip may exhibit abrupt jumps at the ray ζ ∞ . This peculiar behaviour is ultimately related to the structure of the Bethe Ansatz for ∆ > 1 and can be heuristically explained by saying that information about the overall sign of the magnetisation is carried by the heaviest quasiparticles [105]. An abrupt jump in O ζ signals that the expectation value of O x varies on length scales shorter than t (i.e. proportional to t a with a < 1), implying that the transport of O x is sub-ballistic. This is in agreement with the numerical findings of Refs. [91,190], which identified diffusive spin-transport (a = 1/2) for ∆ > 1 and superdiffusive one (a = 2/3) for ∆ = 1. In particular, the former type of transport has later been described in GHD by introducing appropriate subleading corrections to (70) [159,191], while the latter is still subject to intensive research in relation to the observed Kardar-Parisi-Zhang scaling of the profiles [192][193][194][195][196][197].
Non-analyticities also appear in the case of multiple quasiparticle species, which corresponds to integrable models with internal degrees of freedom describing for example particles with spin. Due to their physical relevance, inhomogeneous quenches in these systems have been widely investigated, resulting in applications of GHD to the Hubbard model [63,160,198,199], spinful Fermi and Bose gases [85,161], sine-Gordon model [87], non-linear sigma model [194], and a special point of the two-component Bariev model [200]. In these cases, non-analytic points correspond to either different species of quasiparticles or bound-states of quasiparticles of the same species. In fact, it is natural to wonder whether the presence of different species could be directly inferred from the profiles, or, in other words, if bipartitioning protocols could detect separation effects. For instance, in the case of spinful fermions one could ask whether there exist some local observable whose light-cone profile only shows the effect of one of the two species of excitations.
It turns out that, for bipartition protocols at finite energy densities, non-analyticities of all quasiparticle species are typically present in the profiles of arbitrary observables, so that no strict separation happens [85]. However, some separation effects become manifest in special cases. One of them has been pointed out recently in the study of the Hubbard model, where an interesting phenomenon called clogging emerges for some fine-tuned initial states [63,160]. In essence, clogging consists in the fact that a vanishing charge current coexists with nonzero energy currents (or vice versa), within a finite region of the light cone. Its existence has been proven analytically in Ref. [63] in the case where half of the system is initially at half-filling and at infinite temperature, and it has been numerically observed in the high-temperature regime. In addition, different initial configurations resulting in clogging were studied in Ref. [160], where it was confirmed that it could also take place in the NESS. In the next section, we will study a different case where some separation effects become visible, namely the low-temperature regime.

C. The low-energy limits
The low-temperature regime of the GHD equations turns out to be particularly interesting and has been subject to several investigations over the past few years [83,85,89,161,194,198].
In particular, after a bipartitioning protocol from two thermal states at low (but different) temperatures, the GHD equations become analytically solvable, yielding qualitatively different results depending on whether or not the spectrum of the post-quench Hamiltonian has a gap. In the gapped phase, variations of the profiles are found to be exponentially small in the temperatures and are described by non-trivial functions of ζ [83]. In the gapless  [85]. The initial state is obtained by joining together two thermal states at different (low) temperatures, and the same values of chemical potential and non-vanishing magnetic field.
regime, instead, the leading order contributions for the profiles are polynomial functions of the temperature, which turn out to be universal. In fact, in this limit GHD allows one to recover the predictions of conformal field theory (CFT) [67,174,175]. In order to illustrate this, let us consider for concreteness the Hamiltonian (73) in the gapless regime, which is realized for values (∆, h) in a strip h min (∆) ≤ |h| ≤ h max (∆) (see Fig. 8), and initialise the system by joining together two thermal states with inverse temperatures β L , β R . As β L , β R → ∞ only low energy modes around the Fermi point are expected to contribute and the system is expected to be described by a CFT. This implies that the only relevant quasiparticles are those moving at the Fermi velocity and the profiles of local observables take the form of a "three-step staircase". In particular the NESS values (corresponding to the central step) can be computed exactly in the conformal limit, yielding the following results for the energy density and current [67,174,175] where c and v are respectively the central charge and the speed of light in the CFT (c = 1 for the critical Heisenberg chain, while v is a non-trivial function of ∆ and h). Note that the NESS current (and also the energy density) shows 0. the "non-interacting" structure discussed in Sec. III A: it is the sum of two terms, each one depending solely on the initial state of one lead. Predictions (83) can be recovered from a low-temperature expansion of the GHD equations, which also allows one to access higher-order corrections in β −1 L,R [83]. In fact, GHD reveals modifications to the CFT picture for local observables different from the energy density and current. In particular, at the edges of the light cone of generic observables there appears a region of width T L/R = β −1 L,R where the leading contribution is linear, rather than quadratic, in the temperature. An example is given in Fig. 7, where we report the magnetisation profile for an inhomogeneous quench in the critical Heisenberg chain.
Surprisingly, it was found that such a broadening of light cone could also be described by a universal function, which can be derived via a non-linear Luttinger-Liquid approach [89]. Within this paradigm, one approximates the Hamiltonian spectrum near the Fermi points, using fermionic quasiparticles with dispersion relation where r = ±1, while v and m * are phenomenological parameters respectively associated with the velocity and mass of the quasiparticles. From this approach, one finds that, in a neighbourhood of the light cone quantified by , the leading contribution for the profiles of a given observable O is while d is a constant which depends on the observable of interest. Once again, Eq. (85) can be exactly recovered based on low-temperature expansions of the GHD equations [83]. Note that, also including this correction, the NESS currents are still of non-interacting type.
The technical steps involved with the low-temperature analysis of the GHD equations are largely independent of the details of the specific model considered. However, as we have mentioned, in this limit qualitative differences emerge in models with more than one quasiparticle species. This is best exemplified in the Yang-Gaudin model of spinful fermions [29], where the study of bipartitioning protocols at small temperature reveals spin-charge separation effects [85] (similar features are also observed in Fermi-Bose mixtures [161]). In this case, the profiles of local observables display a five-step form, with two distinct light cones propagating from the junction, see Fig. 7 for an example. Observing profiles with this structure, one can argue that they are produced by two decoupled nonlinear Luttinger Liquids, rather than a single one. It is important to stress, however, that for external magnetic field h = 0 local observables couple the two theories: this is due to the fact that the decoupled Luttinger Liquids do not describe individual spin and charge excitations, but a combination of the two [124,201]. For h = 0, spin and charge completely decouple at the level of the Luttinger Liquid description [201]. However, by setting h = 0 in the postquench Hamiltonian, and constructing the initial state by joining together two Gibbs states at different temperatures, it is not possible to create a magnetization imbalance. Therefore, in the ensuing dynamics the magnetization remains frozen (and equal to zero), without any light cone. In conclusion, within the bipartitioning protocol, it is not possible to observe a genuine separation of spin and charge in the form of two distinct light cones. On the other hand, it was recently shown that GHD is capable to predict such a separation for more general inhomogeneous initial states, where the gas is initially confined in a trap potential [202]. In this case, for vanishing post-quench magnetic field and low temperatures, an initial spin-charge imbalance lead to the formation of two separate light cones for spin and charge, whose real-time dynamics can be quantitatively captured by the GHD equations [202].
Finally, we mention that low-temperature limits of the GHD equations have also been investigated in integrable quantum field theories [87,194], where they allowed to clarify the connection between GHD and the semiclassical approach developed by Sachdev and collaborators [203][204][205][206]. Specifically, Ref. [87] considered bipartitioning protocols in the sine-Gordon model recovering the predictions of Refs. [207,208], based on the semiclassical approach, as a lowenergy limit of the GHD equations. In this limit the transport of the topological charge was found to be sub-ballistic. Away from the low-energy limit, however, the numerical solution of the GHD equations showed that transport is always ballistic, in conflict with the semi-classical predictions [87].

IV. QUANTUM ENTANGLEMENT GENERATED BY INHOMOGENEOUS QUENCHES
The dynamics of quantum correlations are generically very hard to describe exactly, both in homogeneous and inhomogeneous settings. This is essentially due to the fact that they go beyond the purely hydrodynamic description that arises at large times. For example, although GHD gives us the exact asymptotic values of one-point functions after quenches from bipartite initial states, connected equal-time correlation functions between points located at different rays are subleading in the scaling limit (52). There are, however, some exceptions to this empirical fact where non-trivial correlations can actually be accessed. For instance, dynamical correlation functions along ballistic light-cones are indeed accessible within GHD [209,210] (see also the contribution by Doyon, De Nardis, Medenjak and Panfil in the present volume). Moreover, as discussed in Sec. V, one can recover an effective description for time-dependent quantum correlations for quenches starting from a particular class of initial states.
In this section we consider another of such remarkable examples. In particular we show how, retaining some genuine quantum correlations generated at the time of the quench, one can describe the linear growth of several entanglement measures. The key for this to happen is the presence of well-defined quasiparticles protected by integrability. After a quantum quench, EPR (Einstein-Podolsky-Rosen) correlations are created between quasiparticles with opposite quasimomenta. The balistic propagation of these quasiparticles transports these correlations leading to the growth of entanglement. An important remark is that this "quasiparticle picture" for the entanglement spreading applies to quenches, both homogeneous and inhomogeneous, in which the steady state is described by a statistical ensemble with finite density of thermodynamic entropy, i.e., entropy per volume. For stationary states with zero entropy density entanglement-related quantities exhibit a sublinear growth as function of time which is not captured by this approach (see Sec. V).
More specifically, to quantify the entanglement we use the Rényi entropies [211][212][213][214], defined as S (n) where n is a real parameter while ρ A (t) is the density matrix of the system at time t reduced to a subsystem A, i.e.
withĀ the region complementary to A and |Ψ(t) the evolved state of the system. The Rényi entropies characterise the spectrum of ρ A (t), sometimes called entanglement spectrum [214], encoding information on how entanglement is shared between A andĀ. In particular, in the limit n → 1 one recovers the von Neumann entropy Beside their theoretical interest, the quantities (86) are also experimentally relevant. Indeed, in the last few years it has become possible to address the dynamics of entanglement-related quantities with cold-atom experiments [215][216][217][218], and Noisy Intermediate Scale Quantum (NISQ) computers [219].
In the remaining part of this section we show that, combining GHD with a simple quasiparticle picture, one can describe exactly the asymptotic dynamics of the von Neumann entropy in a particular scaling limit. The structure of the section is as follows. In Sec. IV A we introduce the quasiparticle picture. In Sec. IV B this is applied to describe the entanglement spreading after homogeneous quenches, both for interacting and non-interacting systems. In Sec. IV C we discuss the entanglement dynamics after an inhomogeneous quench in free-fermion systems. Finally, in Sec. IV D we consider inhomogeneous quenches in interacting integrable systems.
A. Quasiparticle picture: a semiclassical description of entanglement spreading The quasiparticle picture was originally proposed in the context of CFT [220] to describe entanglement dynamics after global quenches from homogeneous initial states. In essence, the idea is that a homogeneous quench produces an extensive number of quasiparticle excitations, which are responsible for propagating the entanglement throughout the system. The quasiparticles that are produced far apart are assumed not to be entangled with one another, i.e. they do not contribute to the coherent quantum correlations. Only quasiparticles that are produced at the same point in space are entangled. As the entangled quasiparticles propagate through the system, quantum correlations spread accordingly. A further simplifying assumption is that only pairs of entangled quasiparticles are produced, the two members of the pair being emitted with velocities of opposite sign. Finally, entangled quasiparticles travel as free particles, i.e., they do not interact.
Let us now consider the case in which all the quasiparticles have velocity with the same magnitude v. This corresponds to models with perfect linear dispersion like CFTs. Within the quasiparticle picture the von Neumann entropy between a region A and its complement at a generic time t is proportional to the number of entangled pairs that are shared between them. Let us consider a finite interval A of length embedded in an infinite chain. At short time t ≤ /(2v) the number of shared entangled pairs is proportional to the horizontal width of the shaded areas in Fig. 10 (a) at that time, which is 4vt. On the other hand, for t > /(2v) the number of entangled pairs is proportional to . In conclusion, one has that Here s is the "entanglement content" of each pair of entangled quasiparticles (it is the contribution of a single pair times the density of pairs). Note that, due to translational invariance, the only piece of information about A we need to know is the length of the interval, while the position of A within the chain is not important. For this reason we used the notation S (t) in (89). The interpretation of (89) is straightforward, and it is shown in Fig. 9. For 2vt < , S (t) grows linearly. All the pairs that originated in the region of the t = 0 axis shaded by the two light cones are shared between A andĀ. At any time 2vt > the number of shared entangled pairs saturates to an extensive value (S (t) ∝ ). Eq. (89), at this level, has to be regarded as a phenomenological description of the entanglement dynamics. In the next sections we will show, however, that with minimal modifications the quasiparticle picture can be made quantitatively accurate in specific integrable systems.

B. Entanglement dynamics after homogeneous quenches in integrable systems
Let us now discuss how to apply the quasiparticle picture to microscopic integrable systems. We first focus on homogeneous quantum quenches in non-interacting models. To promote Eq. (89) to a quantitative prediction, these have to be fixed from the microscopic data of the integrable model under consideration. As we will see, to do that one only needs "thermodynamic information" about the system. Let us begin by considering non-interacting systems. In this case it is quite natural to associate the entangling quasiparticles with the free modes that diagonalise the Hamiltonian. Unlike (89), such single-particle modes have a nontrivial dispersion, i.e. a mode k has energy ε(k) that depends on k (see, e.g., Sec. II A). Their velocities are then given by the group velocities of these modes which is the same quantity appearing in (47). Note that for non-interacting systems v(k) depends only on the model's dispersion, and not on the pre-quench state.
Adopting the quasiparticle picture, we restrict to the situation where only locally generated pairs are entangled, and we take them to have opposite momentum [221]. This means that correlated pairs are specified by a single momentum k. Assuming an even dispersion relation ε(k) we then find the following generalisation of (89) Now, to completely fix (91) we just have to determine the entanglement content of the pair k, which we denote by s(k) to stress that it now depends of the quasimomentum. This can be fixed by imposing that the integrand of (91) relaxes to the density of thermodynamic entropy of the GGE that describes the steady state after the quench [222]. For free fermionic systems the latter is nothing but the integrand of the Yang-Yang entropy in Eq. (26), namely we have where ρ(k) is the root density of the GGE. Crucially, for quenches in free fermionic systems the validity of (91)-(92) can be proved [223]. Namely, by solving exactly the dynamics one can show that (91) gives the leading order contribution for large t and , and, in particular, it becomes exact in the space-time scaling limit t, → ∞, with /t fixed.
As anticipated before, Eq. (91) contains only thermodynamic information about the system and the quench. Indeed, the group velocity of the quasiparticles is a property of the system's dispersion relation, and the entropy of the GGE is a thermodynamic quantity. For systems with a maximal velocity v M (for example those fulfilling the hypothesis of the Lieb Robinson bound [121]), Eq. (91) predicts a linear growth at short times 2v M t < , followed by a saturation to the volume law scaling S (t) ∝ at asymptotically long times. In contrast with (89) the saturation is not abrupt because the quasiparticles have a dispersion. Note that, while in the original quasiparticle picture only pairs of entangled quasiparticles are assumed to contribute to the entanglement (and such an assumption is verified in several quantum quenches in non-interacting models [223] and for some initial states in integrable interacting systems [135,222], see also [224,225]), this is not true in general. Specifically, for some initial states multiparticle entanglement, for instance involving triplets of quasiparticles, can be present. Some of these situations have been investigated in free-fermion systems, and the quasiparticle picture extended accordingly [226][227][228]. Let us now discuss quenches from homogeneous initial states in interacting systems. In particular, let us consider the case of Bethe Ansatz solvable models with a single species of quasiparticles (see Sec. II E). In this case the quasiparticle prediction (91) is modified as follows [222] where the sum runs over the quasiparticle bound states, labeled by n, while λ denotes the rapidities. Apart from the sum over n, Eq. (93) has the same structure as (91). Also the nature of the entangled quasiparticles is the same as compared to free systems: the quasiparticles in (93) are constructed as the low-lying (particle-hole) excitations around the stationary thermodynamic macrostate, i.e., the GGE, that describes the steady state after the quench. However, the properties of the quasiparticles depend on the specific initial state, unlike the free-fermion case. For instance, the group velocities of the quasiparticles can be calculated by solving the system of integral equations in (68), and now depend on the pre-quench state (because the GGE depends on the initial state). Finally, in analogy with non-interacting systems, s n (λ) in (93) is set to be equal to s YY [ρ n (λ)], the density of the Yang-Yang entropy of the GGE described by {ρ n (λ)}. The latter is given by s YY [ρ n (λ)] ≡ ρ t,n (λ) ln ρ t,n (λ) − ρ n (λ) ln(ρ n (λ)) − (ρ t,n (λ) − ρ n (λ)) ln(ρ t,n (λ) − ρ n (λ)), where, in contrast to the free case (cf. Eq. (63)), ρ t,n (λ) is a nontrivial function of λ. We remark that Formula (93) can also be generalized to the case of integrable models with multiple species of particles [229,230]. Note that, differently from Eq. (91), Eq. (93) is a conjecture. Up to now it has been verified numerically in several quenches in the XXZ Heisenberg model [158,222] and in nested spin chains [230], as well as analytically for some solvable quenches in the quantum cellular automaton "Rule 54" [231][232][233].
Note that, while the quasiparticle picture has been extended to study the Rényi entropies in the steady state [234][235][236], also for quenches from inhomogeneous initial states [237], the full dynamics of the Rényi entropies does not appear to be captured by a formula similar to (93). We also remark that the quasiparticle picture can be easily adapted to describe the dynamics of the mutual information between two intervals [156,158,229,238]. Moreover, it is possible to show that the so-called logarithmic negativity, which is a genuine entanglement measure between two disjoint intervals, becomes the same as the mutual information constructed from the Rényi entropy with Rényi index n = 1/2, cf. (86). This has been checked in free-fermion and free-boson models in Ref. [239]. Finally, we mention that, very recently, the quasiparticle picture has been applied to describe entanglement dynamics in simple free-fermion systems subject to gain and loss dissipation [240].

C. Entanglement dynamics after inhomogeneous quenches: non-interacting systems
Let us now see how the quasiparticle picture can be applied to inhomogeneous setups. We begin considering the case of non-interacting systems. Specifically, we focus on a bipartitioning protocol in the toy model introduced in Section II A. The initial state is obtained by joining two semi-infinite chains (left and right) prepared in two homogeneous initial states. The setup is depicted in Fig. 10. The core of the original quasiparticle picture, namely the idea that entanglement propagates via quasiparticles initially produced in pairs in the bulk of the two chains, remains the same. In contrast to the homogeneous case, however, pairs produced in the left and right chains have different entanglement content. On the other hand, since the system is evolved with a free homogeneous Hamiltonian the velocity of the quasiparticles does not depend on the region where they are produced. Putting all together, one arrives at the following quasiparticle prediction for the entanglement of an interval [x 1 , Note that, because of the absence of translational invariance, we have to keep track of both x 1 and x 2 , defining the boundary of the subsystem under consideration. Accordingly, our notation for the entanglement entropy has changed to S [x1,x2] (t). The integral over k is over the Brillouin zone [−π, π], v(k) is the quasiparticles' group velocity, and s x (k) is the entanglement content of the quasiparticles, which is to be determined.
Before doing that, let us comment on the physical interpretation of (95) (see Fig. 10 (b)). Let us focus on the first term in (95). This describes the entropy contribution of the quasiparticles (members of entangled pairs that are shared with the complement of A) with v(k) < 0 that at time t are within A. They were originated around the edge of the interval at x 2 , and at a generic time t are in the region [max(x 2 + 2v(k)t, x 1 )]. The second term takes into account the quasiparticles with v(k) > 0 that were originated near the edge on the right of the interval at x 1 and are members of an entangled pair shared via that edge. At the generic time t the allowed position of the quasiparticle is in the interval [x 1 , min(x 1 + 2v(k)t, x 2 )].
We now discuss the entanglement content of the quasiparticles. Within the quasiparticle picture there is no creation of entanglement during the evolution, but entanglement is simply created soon after the quench, and then "transported" by the quasiparticles. Indeed, because of the argument x−vt of s x−v(k)t in (95) the entanglement contribution is traced back to the sites where the entangled pairs were produced. As it turns out [241] , the entanglement content s x of the quasiparticles is again fixed by the Yang-Yang entropy transported from the two initial chains Here ρ L/R are the root densities corresponding to the GGEs that describe the bulk of the two chains, i.e., for x/t → ±∞, and s YY is the associated Yang-Yang entropy density (cf. Eq. (92)). The entanglement content of the quasiparticles can be rewritten in terms of the root density ρ x/t (k) associated with the GGE that describes the system at time t and distance x from the origin (cf. (51)). To do so one uses the fact that the Yang-Yang entropy density satisfies the same GHD equation (48) to write where in the last step we used (51). Plugging this in (95), one obtains Importantly, despite the inhomogeneous setup, Eq. (98) still predicts a linear growth at short times, followed by a saturation to a volume-law entropy. The inhomogeneous initial condition is reflected in a non-monotonic behaviour at intermediate times [241]. For later comparison with the interacting case, we specify the above results to the entanglement production rate between two initial semi-infinite chains. Namely we fix x 1 = 0, x 2 → ∞ in (98) and neglect the contribution of the right boundary (which corresponds to choosing open boundary conditions). This yields At this point it is straightforward to rewrite Eq. (99) in terms of quantities depending only on the GGE that describes the interface between the two chains, namely Several remarks are in order. First, the entanglement production rate (i.e., the rate at which the entanglement entropy grows) depends only on the physics at the interface between the two chains, which is described by ρ x/t=0 . This is expected because within the quasiparticle picture the entanglement growth reflects quasiparticles crossing the interface. Second, the quantity on r.h.s. of Eq. (100) is nothing but the rate at which the two chains exchange thermodynamic entropy. Implying that the latter coincides with the entanglement production rate. We anticipate that while the first property remains true for interacting integrable systems, the second one does not.

D. Entanglement dynamics after inhomogeneous quenches: interacting systems
Let us now discuss the entanglement dynamics after bipartitioning protocols in interacting integrable systems. Once again, the main idea of the quasiparticle picture is the same as for free models: entangled pairs of quasiparticles are produced in the bulk of the two chains, with entanglement content given by the density of thermodynamic Yang-Yang entropy in the GGE for x/t → ±∞. Moreover, following the reasoning employed in the homogeneous interacting case (cf. Sec. IV B), here we assume that the quasiparticles are nothing but low-lying excitations on the (quasi) stationary state, which in this case is the LQSSρ s (ζ). This allows us to identify the velocity of the excitations with v n,x/t (λ), obtained by solving (63) and (68) with ρ n (λ) replaced by the root density ρ n,x/t (λ) of the LQSS (cf. Sec. II E 2).
The main complication introduced by this identification is that the trajectories of the quasiparticles are not straight lines. This means that, given a pair of entangled quasiparticles with rapidities λ and −λ forming an entangled pair, it is a nontrivial task to trace back their trajectories and identify the subsystem in which they were produced. This step is, however, necessary to assign to the entangled pair the correct entanglement content. Still, as we are going to show, it is possible to provide compact analytical results for the entanglement dynamics. Moreover, these take a rather revealing form for the entanglement production rate between the two chains.
Let us start by defining as X n,λ (x, t) the position at time t of a quasiparticle created at position x with rapidity λ and bound-state index n. The trajectory of the quasiparticle is determined by the following classical equation of motion [157] d dt X n,λ (x, t) = v n,X n,λ (x,t)/t (λ) .
Since v n,x/t (λ) can be determined for any value of x, t by solving (68), the trajectory of a generic quasiparticle can be obtained numerically from (101) by knowing its initial position. Alternatively, on can use the the so-called "flea gas" method [242] to simulate numerically the motion of the quasiparticles and of the entangled pairs [156]. From the analytical point of view, the integration of (101) seems a daunting task. Yet, in some simple situations it is possible to obtain explicit expressions for the dynamics of the entanglement entropy [157]. Here we do not provide the full derivation of the results but we illustrate the main steps of the reasoning and discuss the physical implications. Let us begin by noting that Eq. (101) can be integrated as with x the initial position of the quasiparticle and t 0 a transient time sufficiently large to ensure the validity of the hydrodynamic description [158]. In order to assign the correct entanglement content to the quasiparticles one has to trace back the quasiparticles' trajectory expressing their final position X n,λ (x, t) in terms of their initial one x. First, from (102) one can derive [157] x where ζ n,λ is the solution of ζ = v n,ζ (λ), while v min = min n min λ v n,−∞ (λ), and v max = max n max λ v n,∞ (λ). Eq. (103) gives the initial position x of a generic quasiparticle in terms of its final one X n,λ (x, t) at time t. The method used above is similar to the so-called method of characteristics to solve first order partial differential equations [243]. The crucial step to derive the quasiparticle picture is to establish a relation between the trajectory of the two quasiparticles forming an entangled pair, i.e., having rapidity λ and −λ. To this end, one defines a function J n,λ (ζ) such that J n,λ (x/t)t is the position at time t of the quasiparticle labelled by (n, λ) that started at the same point in space with the quasiparticle labelled by (n, −λ), which at time t is at position x. We do not provide the explicit expression of J n,λ (x/t), because is rather unwieldy, and we refer the interested reader to the original reference [157].
In terms of J n,λ (x/t) one can express the dynamics of the von Neumann entropy of a given subsystem A = [x 1 , x 2 ] as follows [157] S [x1,x2] (t) = t n dλ {sgn(J n,−λ (ζ 1 ) − ζ 1 )sgn(ζ 2 − J n,−λ (ζ 1 ))(ζ 1 − v n,ζ1 (λ))s n,ζ1 (λ) −sgn(J n,−λ (ζ 2 ) − ζ 1 )sgn(ζ 2 − J n,−λ (ζ 2 ))(ζ 2 − v n,ζ2 (λ))s n,ζ2 (λ)} , where we set ζ i = x i /t. The entropies s n,ζ 1/2 (λ) are obtained by using (63), (72) and the definition of the Yang-Yang entropy density (94). Eq. (104) is the analog of (98) for interacting integrable systems. The first term in (104) takes into account contributions to the entanglement entropy due to quasiparticles at the boundary at x 1 , whereas the second one accounts for the boundary at x 2 . Eq (104) becomes much simpler for the entanglement production rate between two semi-infinite chains. By fixing x 1 = 0, x 2 → ∞ in (104) and neglecting the contribution of the right boundary one finds where we assumed sgn(v n,±∞ (λ)) = sgn(λ), which is the typical situation in bipartitioning quenches [157]. Eq. (105) predicts a linear growth of the entanglement entropy, and depends only on the macrostate with x/t → 0, similar to free fermions (cf. Eq. (100)). An interesting result is that in the presence of interactions the growth rate of the entanglement entropy is different from the exchange rate of thermodynamic entropy, given by S exch = t n dλ|v n,0 (λ)|s n,0 (λ). As noted in the previous subsection this does not happen in the non-interacting case. The origin of this difference is depicted in Fig. 11, and it can be traced back to the dressing of the quasiparticles velocities due to the interactions (see (68)). In particular, Fig. 11 shows the possible trajectories of the quasiparticles forming entangled pairs, which are straight lines outside the light cone, while they are curved inside of it. Three different situations are possible: In the first one (see Fig. 11 (a)) the quasiparticles forming the entangled pair are in different subsystems at t → ∞. Apart from the dressing of the velocities, this is the same situation observed in quenches from inhomogeneous initial states in free-fermion models, and in quenches from homogeneous initial states in interacting integrable ones. However, as shown in Figs. 11 (b) and (d), interactions can be so strong that the quasiparticles forming an entangled pair are in the same subsystem at t → ∞. In particular, in the example depicted in Fig. 11 (b) both members of the entangled pair generated on the left chain are on the right one at t → ∞, while the pair in Fig. 11 (d) is confined within the left chain where it originated. The latter configuration does not contribute to the entanglement. This happens because only the pairs that are shared between A andĀ contribute to their mutual entanglement, and the pairs in Fig. 11 (d) are never shared. On the other hand the configuration in Figs. 11 (b) does contribute to (105). Specifically, an entangled pair as in Fig. 11 (b) contributes to the entanglement until both members of the pair are in the right subsystem. This leads to the factor sign(λ) in (105). An important remark is that the situation in Fig. 11 (b) occurs for "slow" quasiparticles, for which the dressing of the velocities is stronger [157]. In particular, this is true for bound states. Indeed, typically, the larger the bound-state size n is, the smaller is its group velocity. An interesting fact is that for some inhomogeneous quench protocols by tuning the interaction strength, and hence changing the quasiparticles velocities, it is possible to completely suppress the contribution of the bound states to the entanglement dynamics and to transport in general. In this case the scenario in Fig. 11 (d) holds for bound states with arbitrary λ. Moreover, there is a "critical" value of interactions strength below which bound-state transport is permitted. Interestingly, the critical interaction depends on the bound-state size n, and it decreases on increasing n. This means that by lowering the interaction strength it is possible to progressively allow for the transport of bound states with larger and larger n. These effects have been investigated in Ref. [237] (see also Ref. [158]).
Let us finally discuss some numerical checks of (105). We consider once again the XXZ Heisenberg spin chain, and focus on the initial state obtained by joining together the Néel and the tilted ferromagnetic states, which are defined as where we denoted by |↑ , |↓ the two local basis states. We focus on the entanglement entropy between two semiinfinite chains prepared in (107) and (108). In Fig. 12 we report exact numerical data for the dynamics of the entanglement entropy between the two chains obtained by using the time-dependent Density Matrix Renormalization Group (tDMRG) [244]. The figure shows S [x1,x2] (t) plotted versus time. The different lines correspond to different values of the anisotropy ∆ of the chain (cf. (73)) and different tilting angles θ (cf. (108)). Clearly, the numerical data exhibit a nontrivial transient dynamics. Note in particular the presence of oscillating corrections. Still, at moderately long times t ≈ 10 a clear linear behaviour sets in. In all the cases, this appears to be accurately captured by Eq. (105), which is reported as red straight lines in Fig. 12.

V. QUANTUM FLUCTUATIONS AROUND INHOMOGENEOUS BACKGROUNDS
While in the previous section we have assumed a quasiparticle picture to gain insight into the entanglement dynamics, in this section we consider another approach to complements GHD and access exact leading-order quantum correlations out of equilibrium. This approach, first introduced in Ref. [143], is based on the combination of GHD and a recently-developed inhomogeneous Luttinger-Liquid description [177,245,246] (see also Ref. [167]), and allows one to compute exactly generic time-dependent correlation functions when evolving from low-energy states.
The initial idea was that, at low energy, inhomogeneous systems can still be treated using CFT methods at the price of working in a curved background. While the first works were restricted to equilibrium situations [177], the method was later extended to time-dependent problems [246]. Meanwhile, it was realized that, giving up conformal invariance, the same framework allows us to treat much more general situations both in [245] and out-of-equilibrium [143]. These ideas are better illustrated considering a specific model. In particular, following most of the relevant literature, here we focus on the case of the 1D Bose gas. This choice is also experimentally motivated as 1D Bose gases in external potentials are now routinely realised in cold-atom experiments, see e.g. [42,[247][248][249][250] and the contribution by Bouchoule and Dubail to this special issue.
We organize the discussion as follows. In Sec. V A we focus on trapped 1D Bose gases at equilibrium, while in Sec. V B we move to the out-of-equilibrium situation.

A. Luttinger-Liquid treatment of trapped 1D Bose gases at equilibrium
A 1D Bose gas with point-wise interactions -also known as Lieb-Lininger model [28,251] -in the presence of an external potential V (x, t), is described by the following Hamiltonian where µ is a chemical potential and Ψ † = Ψ † (x), Ψ = Ψ(x) are operators that create/annihilate a boson at position x, satisfying the canonical commutation relations There are several differences with respect to the model considered in the previous sections: besides its bosonic nature, this model is defined in the continuum and, for a generic value of g, it is interacting. Moreover, while in the absence of confining potential the model is Bethe-Ansatz solvable [29], V (x, t) breaks integrability.

Homogeneous Luttinger Liquid
In the homogeneous case (i.e. V (x, t) = 0) and at equilibrium, the low energy physics (i.e. the physics on length and time scales that are large compared to microscopic ones) of the system (109) is captured by the Luttinger-liquid (LL) theory [201] or, equivalenty, the CFT of a free compactified boson φ [252]. For later convenience, let us briefly review how such a bosonic field, sometime referred to as "height field" [253], emerges from the microscopic degrees of freedom. In general, any local operator O(x) in the microscopic model can be written as a sum of local operators in the CFT with A O,φj non-universal constants, and {φ j } j CFT operators, which can be ordered from the most to the least relevant one, according to their scaling dimension. If we denote the microscopic density of the physical bosons bŷ ρ(x) = Ψ(x) † Ψ(x), then its leading order in the CFT operators is given bŷ with ρ 0 the average (constant) density. The boson correlations are fixed (in imaginary time) by the following action Here z = x + ivτ ,z = x − ivτ , and τ is the imaginary time. The parameters of the model are K, known as Luttinger parameter and related to the interactions, and v, the sound velocity [201,254].

Inhomogeneous Conformal Field Theory in the Tonks-Girardeau limit
Let us now turn to the case of non-zero external potential. We begin by considering the limit of infinite repulsion, g → +∞, also known as Tonks-Girardeau (TG) limit. This limit describes hard-core bosons or, equivalently (after a Jordan-Wigner trasformation) free fermions. Specifically, by introducing the fermionic operators the Hamiltonian (109) becomes quadratic Note that bosonic and fermionic densities are the same under the mapping (114), i.e.
Following the historical development (cf. Ref. [177]) we begin by considering a gas at equilibrium at zero-temperature in a time-independent potential V (x). We also assume that the inhomogeneity changes "slowly enough", namely we are in the regime of Eq. (53): at this stage, since we are at equilibrium, it is enough to assume it in the "space" direction, which is nothing but LDA [6].
In this regime, the system can be considered homogeneous in a small patch (i.e. a small spatial region) and one can use results from the homogeneous theory. Matching these patches together, however, imposes constraints on the global low energy theory. In particular, Ref. [177] showed that the only consistent field theory defined globally on the entire domain such that its propagator has the required local behaviour everywhere is a CFT of a Dirac fermion or, equivalently, a free boson, in a curved space-time. The associated action reads where g ab is the metric field (a, b ∈ {z,z}) and we used that at the free Fermi point K = 1. Importantly, the latter fact is not modified by the presence of the trapping potential. As in the homogeneous case, the bosonic field φ is related to the density fluctuations, but now around an inhomogeneous average. Namely, at the leading order in the scaling dimensions, the density is now written aŝ where ρ 0 (x) is fixed by the LDA. In our case, where the inhomogeneity comes from a generic external potential V (x), it is given by ρ 0 (x) = µ − V (x). The curved metric in (117) is fixed by such (inhomogeneous) classical limit of the density. In particular, there exists a set of coordinates (known as isothermal coordinates) for which the metric takes the simple form with (z ,z ) given explicitly by In such "stretched" coordinates the system appears again uniform, except for the (overall) conformal factor π 2 ρ 0 (x) 2 in Eq. (119) (which can be eliminated by a Weyl transformation). Namely, the action is still of the form (117), with g ab replaced by the flat (euclidean) metric. The main difference with respect to the homogeneous case, is that now the sound velocity characterizing the CFT is not constant, but depends explicitly on the position, i.e. v(x) = πρ 0 (x). This feature leads to the emergence of curved light cones, as pointed out in [178]. In principle this approach gives access to all the correlation functions. Specifically, in the case of the Bose gas it leads to explicit predictions for the entanglement entropy (both in the ground state [177] and in the low-lying excited states [255]) and the one-particle density matrix Ψ † (x)Ψ(0) [256] for a generic form of the external potential.
Other systems which can be described in terms of a CFT in a curved space are those where the Hamiltonian is deformed via an envelope function f (x), i.e.
Examples studied in the literature include sine-square and Moebius deformations [257], and the so-called rainbow model [258]. In the latter case, the underlying conformal invariance has been exploited to get analytic predictions and explicit results for entanglement-related quantities (i.e. different entanglement measures, entanglement Hamiltonian, and contour [259]).

Inhomogeneous Luttinger Liquid for generic interaction strength
Let us move away from the TG limit, and look at generic values of the interaction strength, g (cf. (109)). Also in this case we expect that the low energy physics is described by a quadratic action. This is ultimately related to the universality of the Luttinger-Liquid description and can be understood using the following argument (cf. Ref. [245]).
First we assume that the local degrees of freedom of the system of interest are captured by a single bosonic field φ associated to a local Hamiltonian. Next, we also assume that physical observables, and hence the action, are invariant under constant shifts φ → φ + const. This automatically leads to an action of the form where L is the Lagrangian density. At this point, assuming that the action is minimised by a unique classical configuration φ cl , we expand at second order around the minimum and define with metric tensor written in the coordinates (x, τ ), and K(x) ≡ [4π √ det ∇ 2 L] −1 is interpreted as a space-dependent Luttinger parameter. All higher order terms in the expansion of the same action have scaling dimension larger than 2, and therefore are irrelevant in two dimensions in an RG sense. We refer to this model as inhomogeneous Luttinger Liquid (or iLL in short): literally, the inhomogeneous (and in general also time-dependent) generalization of the standard Luttinger-Liquid model [254,260]. Note, however, that also other names appeared in literature (e.g., it was dubbed inhomogeneous Gaussian free field in [245]).
The action in (117) is a special example of iLL which is also conformally invariant. This is essentially due to the fact that only the sound velocity of the Luttinger Liquid becomes position dependent, v → v(x), while the Luttinger parameter remains constant. The latter, however, turned out to be a special property of the infinite interaction limit (free fermionic point): for finite interaction also the Luttinger parameter generically acquires a spatial dependence, K → K(x), breaking conformal symmetry [245]. Nonetheless, since the resulting iLL theory is still quadratic, it can be used to numerically evaluate all correlation functions in terms of the two-point function, while the latter is found by solving an appropriate generalised Poisson equation [245]. All the inhomogeneous parameters are again obtained from the microscopic model relying on separation of scales, and using the exact Bethe Ansatz solution of the homogeneous case [261][262][263]. This approach has been used to compute the density profile (including density ripples) [245], the one-particle density matrix [245], and the entanglement entropy [264] in the ground state of the trapped Lieb-Liniger model (and its anyonic generalization [265]) for generic values of the interactions. Finally, we remark that there are special situations where also at finite values of the interactions K remains constant [178], and therefore the model retains conformal invariance. An example is given by the stationary state resulting from a bipartitioning protocol, where two XXZ chains prepared in two ferromagnetic states with opposite magnetizations are joined together. There CFT methods were exploited in Ref. [266] to analytically characterise the leading behaviour of the entanglement entropy. Moreover, in some cases, K has been observed to vary very slowly [267] and therefore it is still possible to rely on conformal symmetry to get analytic results.

B. Luttinger-Liquid treatment of quenches in 1D trapped Bose gases
The inhomogeneous and time-dependent LL point of view outlined above can be extended to non-equilibrium situations. Every time one has a well-defined (semi)classical limit (which at equilibrium is given by the LDA), one can include quantum fluctuations by exhibiting an action, whose minimum reproduces that classical limit, and whose quadratic expansion gives quantum fluctuating corrections. This action, in turn, defines a path integral, and, therefore, a standard "quantization" procedure which has been applied in many different context, from superfluidity [268] to Hall liquids [269]. The problem, then, becomes how to write explicitly an action satisfying the above requirements (whose existence and uniqueness is anyway not guaranteed). Below we are going to illustrate the solution to this problem, as found in [143].
Before proceeding, however, it is worth recalling the main logic behind this quantization. Indeed, our starting point is a fully quantum model (in our case the one described by the Hamiltonian (109)), where correlations are all contained in the time evolved quantum state so that no further quantization is needed at this level. However, in typical situations the full state cannot be constructed explicitly and its correlations are practically inaccessible. The idea is then to circumvent the problem by first finding the (semi)classical limit of the theory (the hydrodynamic solution in our case) and interpreting small fluctuations around such solution as classically propagating linear sound waves. The latter are eventually "re-quantised" by constructing the aforementioned action. Note, in particular, that this is not the action associated to the initial model. In fact, the relation between the two is not generically understood.
As before we will illustrate the method focusing on 1D Bose gases in traps, now subject to a quantum quench to induce the dynamics.

Quantum Hydrodynamics
We begin once again by considering the model in (109) in the TG limit, and focus on a harmonic trap quench. Namely, we prepare the system in the ground state of (115) with and suddenly change the frequency to a different value ω = ω 0 (note, however, that for harmonic traps the quench problem can be solved for generic smooth functions ω(t) [246]). Since in the TG limit the theory is quadratic and the initial state is Gaussian, all observables can be computed in terms of two-point functions. In particular, a very convenient quantity to consider is the so-called Wigner function [270] n x,t (k) = dy e iky Ψ † which fulfils the following evolution equation While this equation is exact for harmonically trapped non-relativistic particles, the same is not true when considering a relativistic system or free fermions on a lattice (see Sec. VI). In addition, from (126) it is clear that n x,t (k) contains the same information as the two-point function Ψ † F (x, t)Ψ F (y, t) , while it does not capture "off-diagonal" contributions (cf. Sec. VI) coming from Ψ F (x, t)Ψ F (y, t) (when it is non-zero) or higher point correlation functions (present for non-Gaussian initial states).
Here, however, we are not interested in the exact microscopic dynamics of n x,t (k). Instead, we interpret n x,t (k) as a coarse grained (or semiclassical) slowly varying distribution function in position and momentum, very much like ρ x,t (k) in Eq. (54). In fact, the evolution equation (127) is nothing but the GHD equation in the presence of a slowly varying harmonic potential [271]. We remark that such a coarse grained interpretation of the Wigner function has been extensively used to characterise the semiclassical regime of Fermi gases in 1D [272][273][274].
Our goal here is to simplify (127) by using the structure of the initial state. To this end we observe that our initial state (the ground state in the harmonic trap with ω(t) = ω 0 ), is a particular example of zero-entropy state, namely, it has zero entropy density (such expression was first introduced in [139]). Those are states whose Wigner function (in the semiclassical limit) reduces to a characteristic function, parametrised by a curve Γ t (the Fermi contour ) Locally, they look like split Fermi seas [275][276][277], labelled by a finite number of Fermi points {k a } a=1,··· ,2Q (see Fig. 13). Importantly, since entropy is conserved at the hydrodynamic level, these states remain zero-entropy states under GHD evolution. Note also, in connection to Sec. IV, that evolution from such states would give zero entanglement entropy within the quasi-particle picture (which captures linear growth only). However, they still display a sublinear growth of entanglement entropy that can be exactly accessed with the method reviewed in the present section [278]. When evolving from zero entropy states the dynamics of the Wigner function is encoded in that of the contour Γ t , or, equivalently, in that of the Fermi points. Therefore the evolution of n x,t (k) can be described in terms of the 2Q evolution equations of the Fermi points {k a } a=1,··· ,2Q , which take the form of Burgers equations [279] In particular, since in our quench both initial and trapping potentials are harmonic, the Wigner function is just an ellipse that rotates in time [278] (see Figure 13 (a)). This means that for any given position x, at any time t there are always no more than two Fermi points. As a result, the two associated Burgers equations can be re-expressed in terms of the local density ρ(x, t) (continuity equation) and the hydrodynamic velocity u(x, t) (Euler equation), thus recovering the conventional hydrodynamic equations Here P = π 2 ρ 3 /(3m) is the quantum pressure (which explicitly depends on , here set to 1) of the Tonks-Girardeau gas at zero temperature, the only "quantum" input at this stage. We stress that the reduction of GHD to standard hydrodynamics is not related to the free nature of the TG gas, and in fact holds for zero-entropy states of generic interacting integrable models, as long as the number of Fermi points is two. In general, however, during the evolution the number of Fermi points might increase, leading to shocks in the solution of the hydrodynamic equations (note that shocks do not occur, instead, in zero temperature GHD, as pointed out in [139]). An action reproducing the Equations (130) has been constructed in Ref. [246] (see also Ref. [279]). The result is again an iLL of the form (123). Interestingly, in the TG limit the action is still of the form (117). Namely, it is a CFT (only the sound velocity is inhomogeneous), with metric explicitly given in terms of the hydrodynamic parameters as follows This action was used to compute several correlation functions, including a closed-form expression for the n-particle density matrix [246]. The same technique has also been used to engineer non-equilibrium scenarios where entanglement entropy shows a non-standard (namely, non-linear) growth [280]. A similar approach has also been used in [281,282] to study transport and correlation functions starting from an inhomogeneous temperature profile. Finally, Ref. [283] proposed a generalisation of the approach to the case of a CFT with random sound velocity: this can be used to elucidate how purely ballistic waves in standard CFT acquire normal and anomalous diffusive contributions.

Generalized Quantum Hydrodynamics
Even in the TG limit, and for evolutions from zero-entropy states, it is possible to find cases where the simplified description (130) does not hold. Indeed, it is sufficient to look at situations where more than two Fermi points occur [139,278]. This is the case, for instance, of a quench from a double-to a single-well potential, see Figure 13 (a)

x f D a z 4 R K U u S K L R a F q S Q Y k 9 n f Z C A 0 Z y g n l l C m h b 2 V s B H V l K F N p 2 R D 8 J Z f X i W t y 6 r n V r 7 W q V + k 8 d R h B M 4 h X P w 4 A r q c A c N a A K D I T z D K 7 w 5 0 n l x p 2 P R W v B y W e O 4 Q + c z x / 5 t 4 2 R < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " l 4 T v r s q C 0 9 j s y e F + b V V m S l p u h O M = " > A A A B 6 n i c b V B N S 8 N A E J U r 1 q / q h 6 9 L B Z B E E o i B T 1 J w Y v H i v Y D 2 l A 2 2 0 2 7 d L M J u x O h h P 4 E L x 4 U 8 e o v 8 u a / c d v m o K 0 P B h 7 v z T A z L 0 i k M O i 6 0 5 h b X 1 j c 6 u 4 X d r Z d s / K B 8 e t U y c a s a b L J a x 7 g T U c C k U b 6 J A y T u J 5 j Q K J G 8 H 4 9 u Z 7 i 2 o h Y P e I k 4 X 5 E h 0 q E g l G 0 0 o P p X / T L F b f q z k F W i Z e T C u R o 9 M t f v U H M 0 o g r Z J I a 0 / X c B P 2 M a h R M 8 m m p l x q e U D a m Q 9 6 1 V N G I G z + b n z o l Z 1 Y Z k D D W t h S S u f p 7 I q O R M Z M o s J 0 R x Z F Z 9 m b i f 1 4 x f D a z 4 R K U u S K L R a F q S Q Y k 9 n f Z C A 0 Z y g n l l C m h b 2 V s B H V l K F N p 2 R D 8 J Z f X i W t y 6 r n V r 7 W q V + k 8 d R h B M 4 h X P w 4 A r q c A c N a A K D I T z D K 7 w 5 0 n l x p 2 P R W v B y W e O 4 Q + c z x / 5 t 4 2 R < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " l 4 T v r s q C 0 9 j s y e F + b V V m S l p u h O M = " > A A A B 6 n i c b V B N S 8 N A E J U r 1 q / q h 6 9 L B Z B E E o i B T 1 J w Y v H i v Y D 2 l A 2 2 0 2 7 d L M J u x O h h P 4 E L x 4 U 8 e o v 8 u a / c d v m o K 0 P B h 7 v z T A z L 0 i k M O i 6 0 5 h b X 1 j c 6 u 4 X d r Z d s / K B 8 e t U y c a s a b L J a x 7 g T U c C k U b 6 J A y T u J 5 j Q K J G 8 H 4 9 u Z 7 i 2 o h Y P e I k 4 X 5 E h 0 q E g l G 0 0 o P p X / T L F b f q z k F W i Z e T C u R o 9 M t f v U H M 0 o g r Z J I a 0 / X c B P 2 M a h R M 8 m m p l x q e U D a m Q 9 6 1 V N G I G z + b n z o l Z 1 Y Z k D D W t h S S u f p 7 I q O R M Z M o s J 0 R x Z F Z 9 m b i f 1 4 x f D a z 4 R K U u S K L R a F q S Q Y k 9 n f Z C A 0 Z y g n l l C m h b 2 V s B H V l K F N p 2 R D 8 J Z f X i W t y 6 r n V r 7 W q V + k 8 d R h B M 4 h X P w 4 A r q c A c N a A K D I T z D K 7 w 5 0 n l x p 2 P R W v B y W e O 4 Q + c z x / 5 t 4 2 R < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " l 4 T v r s q C 0 9 j s y e F + b V V m S l p u h O M = " > A A A B 6 n i c b V B N S 8 N A E J U r 1 q / q h 6 9 L B Z B E E o i B T 1 J w Y v H i v Y D 2 l A 2 2 0 2 7 d L M J u x O h h P 4 E L x 4 U 8 e o v 8 u a / c d v m o K 0 P B h 7 v z T A z L 0 i k M O i 6 0 5 h b X 1 j c 6 u 4 X d r Z d s / K B 8 e t U y c a s a b L
J a x 7 g T U c C k U b 6 J A y T u J 5 j Q K J G 8 H 4 9 u Z 3 3 7 i 2 o h Y P e I k 4 X 5 E h 0 q E g l G 0 0 o P p X / T L F b f q z k F W i Z e T C u R o 9 M t f v U H M 0 o g r Z J I a 0 / X c B P 2 M a h R M 8 m m p l x q e U D a m Q 9 6 1 V N G I G z + b n z o l Z 1 Y Z k D D W t h S S u f p 7 I q O R M Z M o s J 0 R x Z F Z 9 m b i f 1 4 3 x f D a z 4 R K U u S K L R a F q S Q Y k 9 n f Z C A 0 Z y g n l l C m h b 2 V s B H V l K F N p 2 R D 8 J Z f X i W t y 6 r n V r 3 7 W q V + k 8 d R h B M 4 h X P w 4 A r q c A c N a A K D I T z D K 7 w 5 0 n l x 3 p 2 P R W v B y W e O 4 Q + c z x / 5 t 4 2 R < / l a t e x i t > s < l a t e x i t s h a 1 _ b a s e 6 4 = " G c j P u N 6 9 P f n s T N V s K t Y A P W P a z 6 k = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B i y W R g p 6 k 4 M V j R f s B b S i b 7 a Z d u t m E 3 Y l Q Q n + C F w + K e P U X e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y R S G H T d b 6 e w t r 6 x u V X c L u 3 s 7 u 0 f l A + P W i Z O N e N N F s t Y d w J q u B S K N 1 G g 5 J 1 E c x o F k r e D 8 e 3 M b z 9 x b U S s H n G S c D + i Q y V C w S h a 6 c H 0 L / r l i l t 1 5 y C r x M t J B X I 0 + u W v 3 i B m a c Q V M k m N 6 X p u g n 5 G N Q o m + b T U S w 1 P K B v T I e 9 a q m j E j Z / N T 5 2 S M 6 s M S B h r W w r J X P 0 9 k d H I m E k U 2 M 6 I 4 s g s e z P x P 6 + b Y n j t Z 0 I l K X L F F o v C V B K M y e x v M h C a M 5 Q T S y j T w t 5 K 2 I h q y t C m U 7 I h e M s v r 5 L W Z d V z q 9 5 9 r V K / y e M o w g m c w j l 4 c A V 1 u I M G N I H B E J 7 h F d 4 c 6 b w 4 7 8 7 H o r X g 5 D P H 8 A f O 5 w / 8 v 4 2 T < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " G c j P u N 6 9 P f n s T N V s K t Y A P W P a z 6 k = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B i y W R g p 6 k 4 M V j R f s B b S i b 7 a Z d u t m E 3 Y l Q Q n + C F w + K e P U X e f P f u G 1 z 0 N Y H A 4 / 3 Z p i Z F y R S G H T d b 6 e w t r 6 x u V X c L u 3 s 7 u 0 f l A + P W i Z O N e N N F s t Y d w J q u B S K N 1 G g 5 J 1 E c x o F k r e D 8 e 3 M b z 9 x b U S s H n G S c D + i Q y V C w S h a 6 c H 0 L / r l i l t 1 5 y C r x M t J B X I 0 + u W v 3 i B m a c Q V M k m N 6 X p u g n 5 G N Q o m + b T U S w 1 P K B v T I e 9 a q m j E j Z / N T 5 2 S M 6 s M S B h r W w r J X P 0 9 k d H I m E k U 2 M 6 I 4 s g s e z P x P 6 + b Y n j t Z 0 I l K (b) (note that this setup is also relevant from the experimental point of view, see, e.g., Ref. [42,43,284]; for instance, it can model the celebrated quantum Newton cradle experiment [8]). In this case the hydrodynamic regime is still described by the Wigner evolution equation, Eq. (127), and, for zero-entropy states, it is still equivalent to a set of Burgers equations for the Fermi points, Eq. (129) (now four of them, i.e. Q = 2). However, the latter cannot be recast in terms of the two equations for ρ(x, t) and u(x, t) anymore.
To re-quantize the theory we then have to write an action whose equations of motion reproduce the Burgers equations (namely the GHD equations) and whose second order expansion describes (classical) fluctuations around the GHD solution. We refer to this framework as "Generalized Quantum Hydrodynamics" (note that the theory was originally dubbed "Quantum Generalized Hydrodynamics" [143]; the rearrangement of terms is here aimed at reducing some potential confusion in relation to the results reviewed in the next section). In this more general situation, the assumption for the local degrees of freedom to be captured by a single bosonic field φ becomes questionable. Indeed, from the point of view of bosonization [201], it is necessary to associate a chiral excitation (field) to each Fermi point. Therefore, if there are more than two components it is not possible to end up in a field theory of a single bosonic field.
The crucial observation is that, in phase-space, the problem can be recast in that of quantising fluctuations of incompressible regions (see [143] for details), which is well known in the literature on the quantum Hall effect [285][286][287][288]. Borrowing from those results, it is possible to exhibit an action with the desired property, which turns out to be the one of a chiral boson ϕ living on the Fermi contour (equivalently, around a given point x, it appears as the sum of 2Q chiral fields {ϕ a }). Parametrising this contour as (x(s), k(s)), and introducing a density operator which measures the excess number of occupied states around (x(s), k(s)), one finds indeed the following commutation relations Passing to the usual coordinates (x, t), and to the Hamiltonian formalism, the Hamiltonian H GqHD -where the subscript stands for Generalized Quantum Hydrodynamics -generating (129) can be taken as where A ab ≡ ∂ε a /∂k b is known as flux Jacobian (with ε a , k b energy and momentum, respectively, of a small excitation around a given Fermi point) and is diagonal in the TG limit. The program described above can be directly generalised to fully interacting integrable models. In fact, it was first carried out for a 1D Bose gas (109) with generic interaction g [143]. In this case the role of n x,t (k) is played by the "filling function" where ρ t,x,t (λ) and η x,t (λ) are defined in terms of ρ x,t (λ) in (63) and (64) respectively (we suppressed the index n assuming no bound states). To proceed one replaces (127) with the general GHD equation for ϑ x,t (λ) in the presence of external potentials [271], and considers the class of zero-entropy states of the interacting system [275][276][277], now labelled by a finite number of "Fermi rapidities" {λ a } a=1,··· ,2Q . The associated filling function ϑ x,t (λ) takes again the form of a characteristic function, and small fluctuations can be seen as deformations of the Fermi rapidities {λ a }. It has been shown in Ref. [143] that the momenta {k a } of quasiparticles satisfy an exact continuity equation (generalizing the Burger equation valid in the TG limit). Therefore they are the right quantum variables to be quantised also in the interacting case. The final result is exactly Eq. (134), where the flux Jacobian A ab can generically couple different chiral excitations. The Hamiltonian (134), together with the symplectic structure (133), defines the theory of Generalized Quantum Hydrodynamics, which has been applied in both the free [278] and in the interacting [143] case to get results for correlations functions and entanglement entropy. Note that it takes the form of a multi-component, spatially inhomogeneous, time-dependent Luttinger Liquid, but now, importantly, constructed on top of the exact profile (at the Euler scale) given by GHD, as derived from the full quantum initial model. Note, finally, that once again, (134) acquires conformal invariance in the TG limit.

VI. DYNAMICAL GENERATION OF QUANTUM CORRELATIONS
If we try to reintroduce in the GHD equation (70) describing time evolution in the low inhomogeneity limit, we realise that it does not appear explicitly: the equation is essentially classical. Therefore, that equation can only describe how quantum correlations are transported over time, rather than generated. We have seen in the previous section that, for a certain class of initial states, one can describe the propagation of initial-state correlations also in the presence of interactions by means of the theory of generalised quantum hydrodynamics. Instead, we discuss here how quantum correlations can be generated dynamically, even when they are not present in the initial state (for example, when the system is initialised in a product state). To do that we will assigning a precise, quantitative meaning to the space-time dependent root density that we have qualitatively introduced in the previous sections. Consequently, the forthcoming discussion is less elementary than the rest of the review.
We identify two mechanisms for the generation of quantum correlations: 1. There are quantum corrections to the asymptotic generalised hydrodynamic equation (70).
2. The state is not completely characterised by root densities.
Let us start with clarifying Point 1. Despite the asymptotic GHD equation being classical, the full quantum evolution is doubtlessly not classical. In fact, the dependence on in the von Neumann equation disappears practically only if the density matrix is stationary. This indicates that there must be quantum corrections to the GHD equation itself. Point 2 is clear already in the homogeneous limit, as it is equivalent to the statement that not every state is stationary (see Section II B). Trivial as it sounds, this point can nevertheless be a source of confusion: in relevant scaling limits, some observables are described by expressions that depend only on the root densities, a remarkable example being the growth of the entanglement entropy after a quantum quench from a product state -cf. Sec. IV. In this case the initial state is not stationary but the asymptotic expression in Eq. (93) depends on the state only through the root density. This paradox is resolved after understanding two points: (i) quantities that time evolve independently are coupled by the condition that the state is pure; (ii) the growth of entropy is essentially due to dephasing, which gets rid of the quantities that cannot be interpreted as root densities (in the non-interacting case, they are parametrised by the field Ψ defined below) (see, e.g., Refs [222,223]).
Points 1 and 2 are intimately connected and can be better understood by adopting the point of view described in the second part of Sec. II D, i.e., interpreting (70) as the leading order contribution in a low inhomogeneity limit. In this way the terms generating quantum corrections will be identified with the subleading contributions in the limit. Let us describe how this can be done in the case of free fermion systems, like the one described by the Hamiltonian H(h) in Eq. (1).
To access the subleading contributions in (70) one has to first understand what is the meaning of a space-time dependent root density beyond the leading order. To this end, there are two possibilities to envisage (a) Lifting the root density to a function of space and time is just an effective way to capture the asymptotic evolution.
(b) There is a way to define a space-time dependent root density so that it exactly describes the time evolution of a class of states.
In the first scenario it is usually sufficient to define the root density so as to approach the standard root density in the homogeneous limit. An example of such a line of attack is provided by Ref. [123] (and, to some extent, by Ref. [159] in an interacting system), which attempted to go beyond the low inhomogeneity limit by fixing the definition of root density a priori. The main drawback of a similar description is the impossibility to describe time evolution exclusively in terms of the root density. For example, if we enforce definition (31) for a given choice of functions g µ (k), the evolution of the root density will not be closed (see the discussion around (155)).
On the other hand, at least for non-interacting fermions or equivalent spin chain models, Approach (b) can be successfully applied and one can provide a nonperturbative definition of a special class of states described entirely by the root density [122]. Specifically, the root density can be associated with a Wigner function [270] of the Bogoliubov fermions diagonalising the system and the special states are such that no correlations other than those captured by the root density are present. The latter condition is necessary because the relation between root density and correlation matrix is not one-to-one -cf. Point (2). These states extend the notion of locally quasistationary states beyond the low-inhomogeneity limit (59).
In this section we will assume that the Hamiltonian is homogeneous, like that in Eq. (1). A similar approach can however be applied also to fermions described by inhomogeneous quadratic Hamiltonians [122], and, in particular, to free fermions in a trap [289][290][291][292].
Let us consider the 2-by-2 block of the correlation matrix associated to sites ( , n) where c n can be the spinless fermions of our example with Hamiltonian (1). For the LQSSs of the Hamiltonian (1), this can be parametrised as follows where ρ x,e/o (p) can be interpreted as the even and the odd part (w.r.t. p) of the root density, respectively. This matrix has indeed two important properties: first, in the homogeneous limit it describes a generic stationary state, and, second, it is closed under time evolution [122]. We point out that the space dependence in the root density (highlighted by the subscript m/2 in (137)) is the result of a particular convention, explained below, in assigning a position to a product of operators acting on different sites (e.g., c † c n ), which is however irrelevant in the low-inhomogeneity limit. As long as quadratic operators like c c n are concerned, arguably, the most intuitive definition of position is the average x = ( + n)/2. This convention has for example the advantage that the position changes in a simple way under chain inversion. There is, however, a possibly unexpected consequence: quadratic operators with odd range have integer positions, whereas those with even range lie on an effective lattice with half-integer positions. To avoid the complication of specifying which lattice the operator belongs to, we allow x to run over all positions, independently of whether they are physical or not. Given that, the root density ρ x (p) does not need to capture the expectation value of the operators that cannot have position x. And this holds true even for values of x that are neither integers nor half-integers, in which case ρ x (p) could in principle assume any value. Ref. [122] proposed to use these degrees of freedom to promote ρ x (p) to an entire function of x ∈ C. This simplifies the asymptotic expansion in the limit of low inhomogeneity and makes it more transparent the correspondence with the kinetic interpretation of the root density.
Keeping this in mind, we can invert (137) and find ,n∈Z e i(n− ) p +i(2y−n− )q tr Γ ,n e −i 1 2 Θ h (p− q)σ z I + σ y 2 e i 1 2 Θ h (p+ q)σ z , (138) where K(x) is any entire function that equals 0 at nonzero integers and 1 at x = 0 (in fact, one could also consider more generic interpolations, but we aim at simplicity). The choice of K(x) is arbitrary because the root density is defined so as to capture only the expectation values of the operators at either integer or half-integer positions. This degree of freedom becomes superfluous in the homogeneous limit, where we would like to drop any dependence on x. It is then convenient to choose K(x) satisfying two auxiliary properties: and z∈ 1 2 Z We call this (2N + 1)-th order generalised hydrodynamics. In particular, truncating the expansion at the third order (N = 1), setting back b tob, and integrating over q, we find the solution to third-order generalised hydrodynamics ρ (3) x,t (p) = dyAi[y]ρ (3) x−E h (p)t+ y where the physically irrelevant dependence onK is understood and we have also shown the first orders of the lowinhomogeneity expansion using that the dispersion relation is even. In particular, in the low-inhomogeneity limit e 2iE h (p)t/ Ψ x,t (p) (159) satisfies the Schrödinger equation of a free particle with effective mass 2/E h (p). As for the complete GHD equation, also the dynamical equation of the auxiliary field, (158), can be readily inverted Ψ (K) x,t (p) = z∈ 1 2 Z π −π dq 2π e 2iqz e − i (E h (p+ q)+E h (−p+ q))t Ψ (K) x−z,0 (p) .
The solution to (158) truncated at a given order can be obtained in the same way as the solution to n-th order GHD; in particular, expanding (158) up to the third order (which vanishes for even dispersion relations like E h (p)) we have In contrast to what happens for the root density (150), the leading correction to the auxiliary field is captured by a Gaussian kernel.
We are now in a position to understand that a different definition of root density,ρ x (p), mixing the fields ρ x (p) and Ψ x (p) and approaching ρ x (p) in the low-inhomogeneity limit, which can be formally expressed as would satisfy a different dynamical equation, irremediably coupled with independent degrees of freedom. In addition, even neglecting the contributions independent ofρ x (p), the dynamical equation would have a different lowinhomogeneity expansion, which is exactly the problem affecting this definition of inhomogeneous root density.

A. Behaviour at the light-cone edges
In this section we consider the physics around the edges of the light cone when the root density describing the initial state is piece-wise continuous. Note that this does not exactly correspond to the bipartitioning protocols considered so far, where the initial density matrix takes the formρ(0) =ρ L ⊗ρ R : correlations between left and right operators are still present around the discontinuities of the root density (and of the auxiliary field). Note however that a state with a piece-wise constant root density and the corresponding partitioned state are almost equivalent everywhere except in regions, localised around the junctions, of extent proportional to the correlation length. The bipartitioning protocol will be considered at the end of this section. We will follow the approach of Ref. [123] but we will define the root density as in the previous section. Specifically, using the sinc kernel, the correlation matrix is given by Eq. (137) with the root density ρ x,t (p) = 1 4π + lim where we called ± (p) = ρ R/L (p) − 1/(4π ) the shifted root density at time zero at the right and at the left of the discontinuity, respectively. Expectation values around a light-cone edge are completely characterised by the correlation matrix restricted to a moving subsystem that contains the edge. We consider the scaling limit in which the time is large and the displacement from the edge scales as t 1/3 . For pedagogical purposes, we start by assuming first-order generalised hydrodynamics (which, we remind the reader, corresponds to expanding the argument of the exponential in (163) at the first order in q). By carrying out the scaling limit, we find that the correlation matrix behaves as where Γ R is the correlation matrix associated with ρ R (p) (hence, at the right of the light-cone edge) and we defined the dimensionless displacement as function asymptotically behaves as O(t −2 ). This cancellation is convenient, as it allows us to use (173) even in quench settings, where the state evolves in time also in the bulk. Finally, let us clarify the relation between the results presented in this section and the bipartitioning protocolŝ ρ(0) =ρ L ⊗ρ R discussed in the previous sections. To this end, we assume that in the two states joined together connected correlations decay exponentially. Focusing on the part of (155) involving the root density, the correlation matrices of the two settings differ in a term of the form dyAi [y] m∈Z [−π,π] 2 d 2 q π e i( +n−m)q1 δρ where δρ x,0 (p) decays exponentially in |x|. We can express the latter property as a smoothness condition on its Fourier transform: whereρ ω ( q 2 ) is a smooth function of ω. By summing over m in (174) Since all functions are smooth, a change of variables to is sufficient to show that the expression is O(t − 2 3 ), and hence subleading with respect to (168). In conclusion, (168) and (173) describe the light-cone edge even in the standard bipartitioning protocols where there are no initial correlations between operators acting on opposite sides of the junction.

VII. CONCLUSIONS AND OUTLOOK
We have given a pedagogical account of GHD seen as a theory to describe the large-time behaviour of quantum many-body systems after inhomogeneous quenches. We have described the basic ideas of the theory -paralleling the now understood case of homogeneous quantum quenches -and its simplest predictions concerning the expectation values of local observables at infinite times. We have discussed certain extensions of the theory to keep track of the evolution of quantum correlations. Specifically, we have showed how GHD can be combined with the quasiparticle picture of Ref. [220] to describe the entanglement growth after inhomogeneous quenches, and how it can be combined with the inhomogeneous Luttinger-Liquid framework of Ref. [245] to account for the evolution of quantum correlations for a certain class of small inhomogeneous quenches. Finally we have discussed a systematic approach to determine higher order corrections to the GHD equations in non-interacting systems. In spite of the extremely rapid development of GHD (the original papers appeared less than five years ago and the theory is already the subject of a special issue) there are still many open questions. Here we discuss some of them, following their order of appearance in the review.
• Arguably one of the most pressing open questions concerning the development of the theory is the "initial condition" problem for Eq. (70). Namely, there is currently no standard procedure that, given a slowly varying inhomogeneous initial state |Ψ 0 , produces a suitable initial condition for Eq. (70), enabling a quantitative description of the late-time dynamics. As discussed in Sec. II E 2, for interacting systems this can be done only for a very limited class of initial states.
• As discussed in Sec. IV the quasiparticle picture for the entanglement spreading gives us the leading-order description in the space-time scaling limit. An interesting question is to understand subleading corrections. These are expected due to the fact that the quasiparticles follow classical trajectories only on average. In fact, due to the interactions, the motion of a given quasiparticle performs a random walk around the classical trajectory. This effect is responsible for the diffusive correction to the GHD equation [159,191,294] but its role in the entanglement dynamics is not yet understood. A related question is whether and how exotic transport properties of integrable systems at some special points, such as superdiffusion, are reflected in the entanglement dynamics.
Other open issues pertaining to the study of entanglement are (i) describe the full-time dynamics of the Rényi entropies based on a quasiparticle structure, (ii) clarify how integrability breaking affects the entanglement spreading [295], and (iii) incorporate the effects of dissipation in the quasiparticle picture [240,[296][297][298]. Finally, another interesting direction is to investigate the entanglement spreading during the out-of-equilibrium dynamics with a time-dependent Hamiltonian, for instance, by using the approach developed in Ref. [299].
• As discussed in Sec. V, the current construction of Generalized Quantum Hydrodynamics is only valid for zero-entropy states. It would be interesting to extend the theory to include other initial states, such as, e.g., states at finite (low) temperature. In analogy with the standard Luttinger-Liquid theory, we expect it to be feasible and to potentially reveal an interesting competition between thermal and quantum fluctuations. Moreover, as a theory of linear fluctuations, Generalized Quantum Hydrodynamics cannot give access to all quantum corrections. There are different ways to take them into account, corresponding to different kinds of (subleading) contributions. A natural next step, would be to investigate corrections coming from non-linear quantum fluctuations (namely, taking into account the "curvature" of the dispersion relation): this is expected to lead to a generalization of standard non-linear Luttinger-Liquid theory [300]. This is related to the problem of studying the effects of irrelevant and marginal perturbations of the LL theory, which are far from understood out of equilibrium (see, e.g., [301] and references therein). Finally another interesting direction pertains to the experimental validation of such extension of GHD in the context of cold-atom experiments, which would be a natural next step to pursue after the recent experimental confirmations of GHD itself [42,43].
• The non-perturbative formulation of generalised hydrodynamics described in Section VI strongly relies on the applicability of Wick's theorem. The possibility of generalising it to interacting integrable systems is an open question. In particular, in the presence of interactions (in the Hamiltonian), it is unclear whether a class of inhomogeneous states described by higher-order GHD can be at all defined. We also mention that the nonperturbative formulation of Section VI has not yet been extended to the non-interacting time evolution of states that do not satisfy the Wick's theorem. Another aspect that deserves more investigation is the study of connected correlation functions, which, as shown in Section VI, could depend strongly on contributions that are generally subleading when considering the expectation value of local operators.