A short introduction to Generalized Hydrodynamics

These are notes based on lectures given at the 2021 summer school on Fundamental Problems in Statistical Physics XV. Their purpose is to give a very brief introduction to Generalized Hydrodynamics, which provides a description of the large scale structure of the dynamics in quantum integrable models. The notes are not meant to be comprehensive or provide an overview of all relevant literature, but rather give an exposition of the key ideas for non-experts, using a simple fermionic tight-binding model as the main example.


I. INTRODUCTION
The notes are based on lectures given at the summer school on Fundamental Problems in Statistical Physics XV. Their purpose is to give a very brief introduction to the exciting recent developments in Generalized Hydrodynamics [1,2]. They are not meant to be comprehensive or provide an overview of all relevant literature. They are several recent reviews and lecture notes [3][4][5] which provide much more detailed discussions and I strongly encourage the interested reader to consult in particular the following 1. B. Doyon, Lecture notes on Generalized Hydrodynamics, SciPost Phys. Lect. Notes p. 18 (2020).

A. Quantum quenches and experiments
A quantum quench is a particular protocol for driving many-particle quantum system out of equilibrium. It is defined as follows.
1. The starting point is a many-particle system in a large, finite volume L with Hamiltonian H.

2.
The system is then prepared in an initial state |Ψ(0) , or an initial density matrixρ(0), that has non-zero overlaps with exponentially many (in system size) eigenstates of H. The initial state should have good clustering properties 1 and is often taken to be a lowly entangled state.
3. At later times the quantum state describing the system is then given by the solution of the time-dependent Schrödinger equation 4. The objective is to study expectation values of local operators O A in the thermodynamic limit lim L→∞ Ψ(t)|O A |Ψ(t) .
Here we define local operators as acting as the identity outside a finite, connected spatial region in the infinite volume limit. For a quantum spin chain operators of the form σ α1 j1 . . . σ α j where j k ∈ [a, b] with a, b fixed are local. As we will see later locality is a very important concept in non-equilibrium dynamics. I note that this notion of locality is very strong and a weaker notion of quasi-locality [6] sufficies for our purposes 2 Often the Hamiltonian depends on a parameter h such a magnetic field or interaction strength, and a popular way of defining a quantum quench is then to take |Ψ(0) as the ground state of H(h 0 ), and consider time evolution under the Hamiltonian H(h 1 ) with h 1 = h 0 . This corresponds to an instantaneous "quench" of h at time t = 0 from h 0 to h 1 .
The theoretical quantum quench protocol introduced above is inspired by cold atom experiments. In order to make the connection more concrete I now present a brief cartoon of experiments on ultra-cold bosonic Rb atoms carried out e.g. in Jörg Schmiedmayer's group in Vienna [7][8][9][10][11]. The Hamiltonian describing the atoms is to a good approximation where V (r j , t) is a confining potential that is varied in a time-dependent way. The potential is separable in the sense that V (r j , t) = 1 2 mω 2 x 2 j + V ⊥ (y j , z j , t) and V ⊥ can be tuned in such a way that the transverse degrees of freedom are essentially projected to the ground state(s) of the single-particle Hamiltonian By choosing the transverse confining potential to be very tight and having a single minimum the system at time t = 0 can be prepared in a low temperature thermal state of the one dimensional Hamiltonian where the interaction strength c is proportional to g and overlap integrals of the eigenfunctions of the confining potential, see e.g. Ref. [11]. In the experiments V ⊥ (y j , z j , t) is then changed in a time-dependent fashion so that one ends up with a tight double-well potential in the transverse direction, cf. [11] for a theoretical description. Neglecting the higher transverse modes (as they have very high energies) then leads to a Hamiltonian that in second quantization takes the form Here a = 1, 2 label the two wells, cf. The splitting of a Bose-Einstein condensate, as realized by a radial deformation of an initially harmonic potential into a double well [46]. The two gases in the final picture are completely decoupled, with no more overlap between the respective wave functions. Animations of the full dynamics are available online [33].
wherein 1 (0) = 0, 1 (T ) = 1. The control parameter mimics the situation in experiments, where the double well potential is controlled by changing the RF field amplitude through an RF current in a wire. For = 0 we recover the static harmonic potential, whereas = 1 corresponds to a fully separated double well with no wave function overlap between the two halves of the system. Since the Rabi-frequency is strictly positive in experiments we employ the same saturation function as in the previous example (cf. Fig. 4).
As the trapping potential is significantly changed during the splitting the atoms are radially displaced from their equilibrium position in the harmonic trap. Consequently, strong dipole and breathing oscillations are usually observed in experiments. This poses a strong limitation to the use of such systems as interferometers [56]. The minimization of such excitations is therefore one of the main motivations for our optimization.

Numerical simulations: single-parameter control
In a first step we again consider the case where the Rabi-frequency is increased linearly (see Fig. 7a). This procedure is identical to the one that is typically used in experiments [49,53]. At the final time t = T the infidelity has only decreased slightly as can be seen from Fig. 7b. Moreover, the infidelity shows the expected strong oscillations for t > T. A snapshot of the density at time t ⇤ = 22.5 ms is illustrated in Figs. 7c-e, revealing that there is large discrepancy between the computed state and the desired state d .
Next, we consider the result of the optimal control algorithm. We find that, irrespective of the specific choice of 0 , the algorithm always converges to approximately the same minimizer of the cost functional. The corresponding time-evolution of the Rabi-frequency is shown in Fig. 7f. We observe that the Rabi-frequency remains zero for the first few milliseconds. In fact, only about three milliseconds of the optimization time T are used for the transformation of the external potential. This behavior persists even if we increase the optimization time T , with the Rabi-frequency vanishing for an even longer initial period of time. The precise timescale depends on  wherein 1 (0) = 0, 1 (T ) = 1. The control parameter mimics the situation in experiments, where the double well potential is controlled by changing the RF field amplitude through an RF current in a wire. For = 0 we recover the static harmonic potential, whereas = 1 corresponds to a fully separated double well with no wave function overlap between the two halves of the system. Since the Rabi-frequency is strictly positive in experiments we employ the same saturation function as in the previous example (cf. Fig. 4).
As the trapping potential is significantly changed during the splitting the atoms are radially displaced from their equilibrium position in the harmonic trap. Consequently, strong dipole and breathing oscillations are usually observed in experiments. This poses a strong limitation to the use of such systems as interferometers [56]. The minimization of such excitations is therefore one of the main motivations for our optimization.

Numerical simulations: single-parameter control
We illustrate the splitting procedure for N = 2000 atoms and T = 6 ms.
In a first step we again consider the case where the Rabi-frequency is increased linearly (see Fig. 7a). This procedure is identical to the one that is typically used in experiments [49,53]. At the final time t = T the infidelity has only decreased slightly as can be seen from Fig. 7b. Moreover, the infidelity shows the expected strong oscillations for t > T. A snapshot of the density at time t ⇤ = 22.5 ms is illustrated in Figs. 7c-e, revealing that there is large discrepancy between the computed state and the desired state d .
Next, we consider the result of the optimal control algorithm. We find that, irrespective of the specific choice of 0 , the algorithm always converges to approximately the same minimizer of the cost functional. The corresponding time-evolution of the Rabi-frequency is shown in Fig. 7f. We observe that the Rabi-frequency remains zero for the first few milliseconds. In fact, only about three milliseconds of the optimization time T are used for the transformation of the external potential. This behavior persists even if we increase the optimization time T , with the Rabi-frequency vanishing for an even longer initial period of time. The precise timescale depends on the parameters of the trap, as the optimization algorithm tries to find a compromise between longitudinal and radial directions. some initial density matrixρ(t 0 )) that is not an eigenstate of H(t 0 ): in this sense the situation is analogous to a quantum quench. The system is now left to evolve in time governed by the Hamiltonian (6). At a time t 1 the confining potential is switched off and the two clouds of atoms start to expand freely in three dimensions. Eventually they overlap and at a time t 2 the density of atoms is measured. The expansion can be easily modelled as the atoms effectively do not interact. One therefore can integrate the Heisenberg equations of motion for the measured observable (the particle density) backwards and relate it to an operator in the split one dimensional Bose gas at time t 1 . One finds that the measured density is [12] ρ tof (x, r, t 2 ) ≈ 2 a,b=1 where r denote the transverse directions and g ab (t, r, x 1 , x 2 ) are known functions. Repeating the experiment many times then provides access to e.g. the expectation value (7) in the split, one-dimensional Bose gas after a period of non-equilibrium evolution.

B. Homogeneous quantum quenches and local relaxation
We first focus on the simpler case of homogeneous systems where both the post-quench Hamiltonian and the initial state are translationally invariant [13] 3 . The first question we ask is whether after a quantum quench a many-particle system somehow relaxes, i.e. whether if we wait long enough the quantum mechanical probability distributions describing the outcomes of measurements become time independent. This is equivalent to the question whether the double limit lim t→∞ lim L→∞ Ψ(t)|O|Ψ(t) (8) exists for all Hermitian operators O. We note that the order of limits is crucial here. It is easy to see that this limit cannot exist for all observables. Indeed, let |n be the eigenstates of the Hamiltonian describing the time evolution of our system and E n the corresponding energies. Then Now we can choose "observables" that never relax, e.g.
Indeed, we have which shows that the expectation value of this particular observable is a periodic function of time. However, the operator O is typically highly non-local in space. This suggests that we should restrict our attention to local measurements and concomitantly local operators O A . For these our double limit generally exists, i.e.
The physical picture underlying this fact is as follows. As O A is a local operator it acts like the identity outside some finite spatial region A. In the infinite volume limit the complement of A simply acts like a bath on A and eventually leads to relaxation. One can reformulate this observation in terms of density matrices as follows. The density matrix of the entire systemρ(t) = |Ψ(t) Ψ(t)| in our case is a pure state and hence can never become time independent. On the other hand, the reduced density matrixρ describing the region A on which our observable acts (Ā is the complement of A) is a mixed state and hence can become time-independent at late enough times in the thermodynamic limit. A natural question to ask at this point whether it is possible to describe the late-time limits of the expectation values of local operators in terms of a statistical ensemble. In other words, is it possible to find a time-independent density matrixρ SS such that for any local operator O A acting non-trivially only on a finite subsystem A We note that in analogy to equilibrium statistical mechanics (where we can choose from micro-canonical, canonical and grand canonical ensembles)ρ SS is not unique. In order to determineρ SS we employ the following ergodicity principle: under time evolution the reduced density matrix of a finite subsystem will retain the minimal possible amount of local information on the initial state. As we are dealing with an isolated quantum system with a Hamiltonian that has a local density H j one piece of local information that is always retained is the energy density Ψ(0)|H j |Ψ(0) .

Local Conservation Laws
A local conservation law is a Hermitian operator I (n) that commutes with the Hamiltonian of our system and has a density I (n) j that is a local operator as defined above, i.e.
We will be particularly interested in the situation where we have many local conservation laws that are mutually compatible, i.e.
We stress that the conservation laws we have in mind here are extensive. The existence of a local conservation law has important consequences for the steady state density matrixρ SS in translationally invariant cases. By (15) we have Translational invariance then implies that Combining (18) with (17) we conclude that where in the last step we have used that I (n) j are local operators. This tells us thatρ SS retains information about the expectation values of all local conservation laws in the initial state.

Thermalization
As we are dealing with an isolated quantum system energy is always conserved This is the minimal amount of information on the initial state |Ψ(0) that gets retained under the dynamics. If there are no conserved quantities other than energy the system thermalizes at late times after a quantum quench. The steady state density matrix is then given by a finite temperature (equilibrium) ensemble constructed as follows. We define a Gibbs density matrixρ and fix the effective temperature β −1 eff by requiring it to correspond to the energy density established by the choice of initial state Under this choice we haveρ We could have equally well chosen a micro-canonical description where |n are energy eigenstates with energy E n and N is a normalization factor that ensures that Tr(ρ MC ) = 1.
Finally we note that averaging over a micro-canonical shell is not required as long we use a typical energy eigenstate to define our micro-canonical ensemble, which then takes the simple form Drawing an energy eigenstate at random out of our micro-canonical window provides us with a typical state with a probability that is exponentially close (in L) to one.

Non-equilibrium steady states and Generalized Gibbs Ensembles
If we have additional conservation laws with local densities I (n) j the system cannot thermalize because which tells us that the system retains more information on the initial state than just its energy density. What should the ensemble describing the steady state then be? The answer to this question is provided by Rigol et al [14] (see also Jaynes [15]): we should maximize the entropy under the constraints (26). This leads to a generalized Gibbs ensemblê where the Lagrange multipliers λ n are fixed by We note that solving this system of equations is a difficult task in general. Alternatively we may employ a generalized micro-canonical ensemble [16,17]:ρ where |Φ L be a simultaneous eigenstate of all local conservation laws I (n) such that

Approach to the steady state
In integrable models how fast the expectation value of a given operator approaches its steady state value depends on its locality properties relative to the elementary excitations in the model [18]. Expectation values of operators that are local relative to the elementary excitations approach their stationary values in a power-law fashion [13]. The spatial "size" of the operator under consideration 4 together with the maximal propagation velocity of elementary excitations v max provides a time scale t = /2v max after which relaxation to the steady state begins [19]. The rule of thumb is that the more local an operator is, the faster its expectation value relaxes towards its steady-state value.

C. Summary
Integrable models with short-range interactions prepared in homogeneous initial states that have good clustering properties relax locally to non-thermal stationary states, which are completely specified by the expectation values of the (quasi)local conservation laws. In terms of the time-evolving density matrixρ(t) of the system this statement reads whereĀ is the complement of an arbitrary, finite subsystem A. If the following we will be interested in the situation where our initial density matrix is inhomogeneous. The idea there is that we still have local relaxation, but in a space and time dependent way, i.e.

II. BRIEF INTRODUCTION TO QUANTUM INTEGRABLE MODELS I: FREE THEORIES
The simplest integrable models are free theories and arguably the simplest free theory is the tight-binding chain where c j and c † j are fermionic creation and annihilation operators on site j and Imposing periodic boundary conditions quantizes the momenta The momentum space operators obey canonical anticommutation relations of the form The 2 L energy eigenstates are conveniently expressed in the momentum representation A. Local conservation laws and associated currents It follows from the representation (33) that all mode occupation numbers are conserved They are in one-to-one correspondence with mutually commuting, extensive conservation laws with local densities by We The mode occupation operators can be expressed as linear combinations of the conservation laws aŝ The Heisenberg equations of motion for the densities of the conservation laws take the form of continuity equations B. Macro-states In the thermodynamic limit physical properties are described in terms of macro-states defined as follows. Consider L, N very large with D = N/L fixed and let 0 ≤ ρ(k) ≤ 1 be a given function. The family of energy eigenstates is called a macro-state. We have We see that the extensive parts of the eigenvalues only depends on the density ρ(k), i.e. is the same for all micro-states. Conversely, if we fix the expectation values of all local conservation laws in some quantum state in the thermodynamic limit then relation (41) fixes a unique macro-state (46)

Counting micro-states
There are clearly many ways of fulfilling (43) for a given density ρ(k). The total number of possible k j values in the interval [k, k + ∆k] is ∆n vac = [L∆k/2π], where [x] denotes the integer part of x. Of these "vacancies" ∆n p = [ρ(k)L∆k/2π] are occupied. The number of ways of distributing ∆n p particles among ∆n vac vacancies is This first of all shows that in general exponentially many (in L) micro-states correspond to each macro-state. Using that the entropy our macro-state is given by S = ln(# of micro-states), the contribution arising from reordering momenta in the interval [k, k + ∆k] is where we have used Stirling's formula in the second step. As we are interested in very large system sizes we therefore have where we have defined the density of holes as ρ h (k) = 1 − ρ(k).

Typical vs atypical states
It is clear from the above construction that any positive function ρ(k) gives rise to a macro-state and these generically have finite entropy densities in the thermodynamic limit. Importantly these macro-states are in general not thermal. The thermal states are obtained by extremising the free energy per site This gives Thermal states are by construction the maximal entropy states, i.e. the most likely states, at a given energy density.
In free theories other macro-states exist at the same energy density, but their entropies are smaller than the one of the thermal state. This means in particular that if we select a micro-state with a given energy density at random, this state will be thermal with a probability that is exponentially close (in system size) to one. So typical states at a given energy density are thermal, but there are exponentially many "atypical" states as well, which differ from the thermal state by the values of the higher conservation laws and hence have different local properties. This situation generalizes to interacting integrable models [20], where such atypical states can have very interesting properties [21].

C. Expectation values of local operators in the thermodynamic limit
Let |k 1 , . . . , k N be a micro state corresponding to a macro state with density ρ(k). Then expectation values of fermion bilinears depend only on the macro-state up to finite-size corrections By Wick's theorem this fact gets lifted to any multi-point correlation function involving a fixed, finite number of fermion operators. This in turn means that expectation values of any finite number of fermion operators calculated in different micro states corresponding to the same macro-state differ only by finite-size corrections that vanish in the thermodynamic limit. Given a local operator O j we introduce the following notations for later convenience D. "Excitations" over (finite entropy density) macro-states Let us consider a micro state |k 1 , . . . , k N corresponding to the macro-state with density ρ(k). We can consider "excitations" over the micro state by composing elementary particle and hole excitations: The energy and momentum of this excitation relative to those of |k 1 , . . . , k N are where the dispersion (p) is given by (34).
The energy and momentum of this excitation relative to those of |k 1 , . . . , k N are In contrast to excitations over the ground state the energy of these "excitations" can be either positive or negative.
In the thermodynamic limit we obtain families of excited states parametrized by their particle and hole momenta. Importantly they propagate with a group velocity that is simply the derivative of the "bare" dispersion relation v(p) = (p) .
This is a special feature of free theories and is not the case for interacting integrable models.

III. BRIEF INTRODUCTION TO QUANTUM INTEGRABLE MODELS II: INTERACTING THEORIES
We now want to generalize the results discussed above for free theories to interacting integrable models, i.e.
2. Construct simultaneous eigenstates of H and Q (n) .
3. Construct macro-states in the thermodynamic limit.
4. Work out stable excitations over these macro-states and determine their group velocities.
The example we will work with is the δ-function Bose gas [22], also known as the Lieb-Liniger model [23] Here Φ(x) is a complex Bose field with canonical commutation relations [Φ(x), Φ † (y)] = δ(x − y). In first quantization the Hamiltonian for N bosons reads It is customary to work in notations where = 1 = 2m and we do this from here on. The Lieb-Liniger model is integrable and local conservation laws can be constructed using the quantum inverse scattering method [22], and the first few read [24] The Heisenberg equations of motion for the densities of these conserved quantities take the form of continuity equations The first two currents are A. Simultaneous eigenstates of Q (n) Simultaneous eigenstates of the Hamiltonian and the conservation laws are constructed as follows. Working in the position representation the time-independent Schrödinger equation for N -particle states reads When the separation between any two particles is large the solution is simply a superposition of plane waves. The structure of the interaction term is such that the plane-wave solutions can be consistently matched whenever two particles meet such that the exact eigenfunctions take the form of the celebrated Bethe ansatz where the normalization is If the model were not integrable we would have a superposition of plane waves with a set single-particle momenta {k 1 , . . . , k N } in the sector x 1 < x 2 < · · · < x N , but with a different set {p 1 , . . . , p N } in the sector x 2 < x 1 < x 3 < · · · < x N , such that the total energy and the total momentum are the same The presence of the higher conservation laws Q (n) ensures that the set of single-particle momenta does not change between different sectors. The states |χ λ corresponding to the wave functions (66) are in fact simultaneous eigenstates of all conserved charges As we need a description of the system in a large, finite volume we now impose periodic boundary conditions As usual these lead to the quantization of the single-particle momenta, but unlike in the free case the resulting quantization conditions are very non-trivial This set of equations is referred to as Bethe equations. Importantly, as a result of the interactions the "single-particle momenta" k j are state dependent and correlated with one another via the quantization conditions (70). Taking the logarithm of these equations we obtain where I j are integers (half-odd integers) for N odd (even). The key point is that each set {I j } of distinct (half-odd) integers is in one-to-one correspondence with a solution {λ j } of the Bethe equations, which in turn provides the wave-function χ λ (x 1 , . . . , n N ) of a simultaneous eigenstate of the Hamiltonian and the conserved charges. This set of states is complete.

B. Macro-states
We now want to construct macro-states along the same lines as for free theories. The main complication is that the quantization conditions (70), (71) are non-trivial and state-dependent. The way to deal with this complication is to work with the (half-odd) integers I j , which are independent from one another. In analogy 6 with (43) we may define a density for z j = I j /L through A positive function χ(z) specifies a macro-state, and corresponding micro-states can the obtained by choosing sets {I j } distributed according to χ(z). In practice we require a formulation in terms of the "particle" distribution function of the roots λ j of (70) defined by This can be related to χ(z) by turning the sum over roots in (71) into an integral over the particle density ρ(λ) in the thermodynamic limit In the thermodynamic limit we have The function z(λ) is called counting function. It is easy to see that z(λ) is a strictly monotonically increasing function.
We are now in a position to relate χ(z) to ρ(λ) by equating the numbers of roots and integers in corresponding intervals The function ϑ(λ) is known as occupation function. This gives It is customary to define a hole density ρ h (λ) by This establishes that rather than specifying a macro-state by the corresponding function χ(z), we can specify it through its root density ρ(λ).

Expectation values of conserved charges in the thermodynamic limit
Let |λ 1 , . . . , λ N be a micro state corresponding to a macro-state with root density ρ(λ) 7 . By translational invariance the densities Q (n) (x) of the conserved charges fulfil In the thermodynamic limit these expectation values do not depend on the choice of micro-state. For later convenience we introduce the notation For local operators these expectation values are independent of the sequence of micro-states chosen in the limiting procedure and depend only on the root density ρ(λ) that specifies the macro-state [22,25,26].

C. Stable excitations over macro-states
We now turn to the construction of excitations over macro-states. Our discussion follows the textbook [22]. We start by specifying a macro-state through its root density ρ(λ). We then determine the counting function from (75), invert it, and determine χ(z) from (76) Finally we use χ(z) to generate a distribution of I j /L that corresponds to our given macro-state 8 . We can then construct families of "excited states" as we now explain for the particular example of a "particle-hole" excitation. The latter is obtained by changing one of the (half-odd) integers I j as shown in Fig. 3 gives rise to a two-parameter family of eigenstates. The Bethe equations for the eigenstates characterized by {I j } and {Ĩ j } are respectively The two equations involving I p and I h simply fix the positions of the corresponding momenta λ p and λ h . Turning sums into integrals and using that both micro states correspond to the same macro-state, i.e. are described by the same root distribution ρ(λ), we have Focusing on the subset of N − 2 equations for which I j =Ĩ j and taking differences we have Next we use that λ j − λ j = O(L −1 ) to Taylor-expand, which gives Defining a shift function by we can turn (85) into an integral equation in the thermodynamic limit We can write this as where Using the shift function we can work out the difference in eigenvalues of any of the conserved charges (68) In particular n = 1, 2 correspond to the excitation momentum and energy respectively. We stress that the latter can be positive or negative. We see that ∆Q (n) can be expressed in the form Crucially, we observe that the contributions of the particle and the hole excitation are additive. It also should be clear that the above construction straightforwardly generalizes to excitations involving multiple particles and holes. Specifying n = 1, 2 in (91) tells us that the energy and momentum of a particle with rapidity λ p are given by q 2 (λ p ) and q 1 (λ p ) respectively. The corresponding group velocity is then v ρ (λ p ) = q 2 (λ p ) q 1 (λ p ) . (92) Here we have introduced the index v ρ to indicate that, in contrast to the case of non-interacting models, the group velocity now depends on the macro-state under consideration. Inspection of (92) and (91) shows that the group velocity in fact depends on the macro-state only through the occupation function ϑ(λ) (76).

IV. HYDRODYNAMIC DESCRIPTION OF NON-EQUILIBRIUM DYNAMICS
In order to see how the non-equilibrium dynamics of many-particle quantum systems relates to classical hydrodynamics we follow Ref. [3].
Classical hydrodynamics in 1+1 dimensions describes the dynamics of fluids on intermediate time and length scales. It combines the continuity equation expressing mass conservation with the Euler equation Here ρ(x, t) is the density of the fluid, v(x, t) the velocity field and the pressure P is assumed to be a function of the density only. Classical hydrodynamics can be reformulated in terms of continuity equations for conserved quantities by introducing the momentum p(x, t) = v(x, t)ρ(x, t) and its associated current j(x, t) = P + v 2 (x, t)ρ(x, t) In these lectures we want to discuss how to obtain similar descriptions for many-particle quantum systems. The general idea is as follows. Letρ(0) be the initial density operator of a many-particle system with (local) time-independent Hamiltonian H. Then the time evolution operator is and the Schrödinger equation implies that the density matrix at time t iŝ The quantities of interest are e.g. expectation values of local operators where we have defined the Heisenberg picture operator As the time evolution (96) The Heisenberg equations of motion for the density of the conservation laws take the form of continuity equations where J (n) (x) are currents associated with the conserved quantities Q (n) . Let us now consider the expectation values in some quantum mechanical state The basic idea is that we now choose to look at our system at length and times scales that are large compared to some mesoscopic "Euler" scale. We will then assume that our system has relaxed locally to a state that by our "minimal information principle" is fully characterized by the expectation values of the conserved quantities, i.e. the q n (x, t). This implies in particular that and hence Assuming that the matrix A is diagonalizable and changing variables (locally) from q n (x, t) to ϑ n (x, t) defined by R n,m = ∂qn ∂ϑm we arrive at evolution equations for the "hydrodynamic normal modes" The difficulty in this description is hidden in the fact that the normal model velocities v n (x, t) depend on the conservation laws in an unknown way. As we will see, it is however possible to work them out exactly in the case of integrable models. This then allows one to determine ϑ n (x, t), which by our assumption of local relaxation/minimal information completely fixes the density operator and thus expectation values of local operators! A. "Derivation" of GHD in free theories using a local density approximation In a homogeneous macro-state with density ρ(k) we have (using translational invariance) in the notations of (53) Here we have used (39) and (42). Now imagine that we prepare our system in an equilibrium density matrixρ that is not homogeneous but varies very slowly in space. Locally this density matrix will "look like" a macro-state with a position-dependent density ρ x,0 (k), x = ja 0 with a 0 the lattice spacing, and in particular Tr ρQ Tr ρJ Let us now turn to the situation at sufficiently late times after our quantum quench. We expect that our system will have relaxed locally, and as long as the system varies very slowly in space we again can employ a description in terms of an appropriately chosen macro-state, which however now will depend on x and t Tr ρQ Using the equations of motion (42) for the charge densities in (108) we obtain an evolution equation for ρ x,t (k) As the (n,α) (k) form a basis of periodic functions on the interval [−π, π] the term in brackets must vanish, i.e.
These are indeed the GHD equations for free fermions. The argument above is based on the principle of local relaxation. We stress that there is no finite time scale over which integrable models relax locally, instead (108) become exact only in the so-called Euler scaling limit x, t → ∞, x/t fixed. At finite times there are power-law corrections to (108) and at sufficiently short times (108) generically do not hold.

B. Dynamics for weakly inhomogeneous initial states in free theories
The above "derivation" was rather hand waiving, and we now want to do a more proper job. The solution of the Heisenberg equation of motion for the fermion annihilation operator in momentum space is Using this we can calculate the two-point functions for an arbitrary initial density matrixρ(0) Assuming thatρ(0) is Gaussian all higher point correlators can then simply be obtained using Wick's theorem. This is the full description of the dynamics. In order to analyze weakly inhomogeneous initial states we introduce a Wigner function (for L → ∞) by W (z, p; t) = Tr ρ(0)Ŵ (z, p; t) , The Wigner function is defined for half-integer values of z and has the following properties: [W1.] The single-particle Green's function can be recovered from the Wigner function using [W2.] For homogeneous statesρ(0) the Wigner function reduces to the density of states in momentum space [W3.] For statesρ(t) that are slowly varying in space and have short-ranged correlations the Wigner function can be interpreted as a position-dependent density of states in momentum space. To see this we note that the Wigner function can be expressed as (116) Let A be a region centered around z such that its linear size |A| is large compared to the correlation length ξ(z) of the stateρ(t), but small compared to the scale ζ(z) on whichρ(t) becomes inhomogeneous in space. Let's further define an integer m(z) such that ξ(z) m(z) ζ(z). All length scales depend on the position z. Then we have As within |A| our state is by assumption approximately homogeneous and has a finite correlation length this corresponds to a position-dependent density of states in momentum space, as advertised.
[W4.] The Heisenberg equations of motion for the Wigner function are linear Given some initial values W (z, p; 0) the solution to (118) is obtain by Fourier methods

Euler scaling limit
Let us now consider the Euler scaling limit and assume that exists. Then we have from (118) This has the same form as the evolution equation (110), and by virtue of [W3] ρ ξ,τ (p) can indeed be thought of as a "ray-dependent" mode occupation function. We note that the solutions to (122) have simple expressions in terms of the "initial conditions" at τ = 0 In order for the Euler scaling limit to be non-trivial (in the sense that one obtains ray-dependent results) the initial momentum space Green's function Tr[ρ(0)c † (p)c(q)] must exhibit some singularities as (112) can otherwise be determined by a simple stationary-phase approximation with saddle points only along the real axis.

Expectation values of local operators in free theories
Let us consider expectation values of local operators along rays in spacetime. The basic building block is the single-particle Green's function G(j, , t) = Tr ρ(t)c † j c = dp 1 dp 2 (2π) 2 e −ip1j+ip2 Tr ρ(0)c † (p 1 , t)c(p 2 , t) .
Let us then go over to the scaling variables (120) and change integration variables to In the limit Λ → ∞ we have We then turn the sum over z into an integral which gives the final result This is precisely what our GHD "phenomenology" predicts: Interpreting ρ ξ,τ (p) as a mode occupation function along the ray characterized by ξ and τ we can construct a corresponding macro-state |ρ ξ,τ by a family of Fock states with momenta distributed according to ρ ξ,τ (p). GHD then predicts that lim Λ→∞ G(j, ; t) = ρ ξ,τ |c † j− c 0 |ρ ξ,τ = dp 2π e −ip(j− ) ρ ξ,τ (p) .
As shown below in the scaling limit the anomalous Green's function vanishes. This is expected [13] because the time evolution operator has a U(1) symmetry c j → c j e iϕ related to particle number conservation. The initial density matrix may break this symmetry, resulting in a non-vanishing anomalous Green's function F (j, ; 0). However, as the system relaxes locally the U(1) symmetry gets restored, which is why F (j, ; t) must vanish at late times. This means that the expectation value of an arbitrary local operators along a given ray can be calculated by Wick's theorem using the regular Green's function only. For example, along the ray t = Λτ ,

Anomalous Green's function
We can see that in the scaling limit the anomalous Green's function vanishes by considering We have • The anomalous Green's function can be recovered from g using • The Heisenberg equations of motion give • Going over to scaling variables this becomes which shows that W A (Λξ, p; Λτ ) vanishes in the scaling limit. Going beyond the scaling limit we see that the anomalous Green's function in position space in fact vanishes as a power law in t, i.e. slowly. This is important, as it imposes a restriction on the applicability of GHD: it will only work at times when the anomalous Green's function is small 9 .

Application of GHD away from the scaling regime
As we have seen above, in the Euler scaling limit GHD becomes exact. In practice one of course would like to apply it away from the scaling limit as well. To that end one writes the GHD equation in the original variables and then uses that the solution of this equation is of the form The "initial" value ρ x,0 (p) is then approximated by the actual initial conditions Finally, the single-particle Green's function is approximated as Using Wick's theorem this allows us to obtain approximate results for multi-particle Green's functions as well.
We will see below how well this works in an explicit example.

Corrections to GHD
For free fermions it is possible to obtain a systematic expansion around the GHD limit [27,28]. In order to do so we carry out a gradient expansion of (118) beyond the leading order, i.e.
Importantly, we also must account for the fact that away from the scaling limit the anomalous Wigner function does not vanish.

C. Example: Free expansion
Let us now consider the explicit example, where our initial density matrix correspond to a product state in position space where all sites j ≤ 0 are occupied We time-evolve this state with the tight-binding Hamiltonian (33). This problem is nice as it can be solved analytically [29,30]. Using that one finds where J n (z) are Bessel functions. Using the addition theorem for Bessel functions this can be simplified Similarly one has for j > 0 and n ≥ 1 The application of GHD proceeds as follows. We first determine the initial conditions for the Wigner function where θ H (x) is the Heaviside function. Following the prescription (138) we then use this to fix Here we have replaced the discrete lattice co-ordinate z by a continuous variable x. The rationale for doing this is that GHD is expected to work on a length scale that is large compared to the lattice spacing. The solution of the GHD equations is then Finally we use the relation between the Wigner function and the Green's function that holds in the scaling limit to obtain an approximate expression away from the scaling regime We similarly can determine the expectation values of the densities of the conserved charges (39), e.g. q (2n,0) (j + n, t) ≡ lim The result for the density profile in the Euler scaling limit is where v max = 2J is the maximal group velocity. In Fig. IV C we compare this asymptotic result to the density profile along rays j/t =fixed at finite times Jt = 10 and Jt = 100. We observe that the GHD result provides a rather good approximation already at times Jt ∼ 10. To get a more precise understanding how the exact result approaches GHD in the scaling limit we show the expectation values of q (2n,0) (10, t) for n = 1, 2 in Fig. IV C. Initially these expectation values are zero because there are no fermions at any of the sites j > 0. At a time v max t * n = 10 + n the propagating front arrives at site j = 10+n and the expectation value q (2n,0) (10, t) becomes sizeable. At later times the expectation values approach the asymptotic GHD result in a power-law fashion. In particular there is no characteristic time scale for local relaxation: the system relaxes locally in a power-law fashion. This is expected to be a generic feature for integrable models.

V. GHD IN INTERACTING INTEGRABLE MODELS
For interacting integrable models we proceed along the same lines as in the free fermion case: 1. Identify local conservation laws Q (n) = j Q (n) j .
4. Work out stable excitations over macro-states and determine their group velocities.
5. Apply a "local density approximation" to the continuity equations for the densities Q (n) j . We'll consider the implementation of this programme for the particular example of the Lieb-Liniger model. The first four steps are discussed in section (III). In a homogeneous macro-state with density ρ(k) we have shown in section (III B 1) that in the thermodynamic limit we have The expectation values of the current densities are more difficult to determine. A simple physically intuitive guess is where v ρ (λ) is the group velocity of the stable particle excitations with momentum p over the macro-state with density ρ(λ), cf. eqn (92) in section III C. The main difference to the case of free fermions is that this velocity now depends on the macro-state under consideration. It turns out that the guess (153) is indeed correct [31,32]. Now we proceed as in the case of free theories and first consider the case where we prepare our system in a density matrixρ(0) that locally corresponds to a homogeneous macro-state but varies slowly in space. Locally it will "look like" a macro-state with a position-dependent density ρ x,0 (k), so that Next we consider the situation at sufficiently late times after our quantum quench, where we assume that our system has relaxed locally and varies very slowly in space. Then we again can employ a description in terms of an appropriately chosen macro-state, which however now will depend on x and t Tr ρ(0) Q (n) (x, t) = Tr ρ(t) Q (n) (x) ≈ dλ 2π q (n) (λ)ρ x,t (λ) , The evolution equation for ρ x,t (k) then follows from combining (155) with the continuity equations (62) relating the expectation values of the charge and current densities, which gives dλ 2π q (n) (λ) ∂ t ρ x,t (λ) + ∂ x v ρ x,t (λ)ρ x,t (λ) = 0 , n = 0, 1, 2, . . .
As the q (n) (λ) form a basis of functions on the real line these imply the following form for the GHD equations in the Lieb-Liniger model We stress that the group velocities v ρ x,t (λ) depend on ρ x,t (λ) in a non-linear way. As we have observed before, the group velocities depend on the macro-state only through the occupation function ϑ x,t (λ) (76). Moreover, for a given root density we can determine the corresponding occupation function from (77), and vice versa. Using this relation between ρ x,t and ϑ x,t we can therefore obtain an equivalent formulation of GHD in terms of the latter A. Integrating the GHD equations In order to use GHD one requires an initial value for the occupation function ϑ x,0 (λ). In practice this imposes serious restrictions on the situations that can be handled. This is easy to understand: from a physical perspective the occupation function encodes the expectation values of the densities Q (n) (x) of the infinite number of conserved charges. These are known in thermal equilibrium, and using a local density approximation in presence of a potential V (x) can be determined in a broad class of inhomogeneous thermal equilibrium states. Certain other situations can be handled by the method introduced in Ref. 33 for homogeneous quantum quenches with lowly entangled initial states, but the problem of determining ϑ x,0 (λ) given a general initial state remains unsolved.
Given an initial value ϑ x,0 (λ) the GHD equation (157) needs to be integrated using that the group velocities are given in terms of ϑ x,t (λ) through (92). An open source Matlab code for carrying out this programme was introduced in Ref. 34, see https://integrablefluid.github.io/iFluidDocumentation/ Once ϑ x,t (λ) and ρ x,t (λ) have been obtained numerically they can be used to determine physical properties. The simplest ones are the expectation values of the charge and current densities (155). Extensions to certain other operators are also possible [35].

B. Extensions of GHD
GHD has been extended in many interesting directions, and we refer the interested reader to the reviews collected in the recent special issue Ref. 5. One such extension that is very important for applications of GHD to cold-atom experiments is to approximately take into account the effects of an external potential V (x) in the Lieb-Liniger model. This modifies the Hamiltonian to and importantly breaks integrability. Hence GHD no longer provides an exact description in the Euler scaling limit. However, as long as V (x) can be considered as a small perturbation, it makes sense to consider how it modifies the GHD equations. The answer was derived in Ref. 36 The obvious question is on which time and length scales this evolution equation provides an accurate description of the non-integrable dynamics induced by (159). The answer is not known, but on physical grounds one would expect that if the initial state of the system is only weakly inhomogeneous and the potential V (x) is weak as well, (160) may provide a good approximation at sufficiently short times. This is because eventually the fact that V (x) breaks integrability will be felt and a description based on continuity equations for conserved quantities must fail. On the other hand, as we have seen above, in general GHD is expected to work only at sufficiently late times unless the initial state is judiciously chosen [37]. It is important to keep these considerations in mind when applying (160).

VI. APPLICATIONS OF GHD TO COLD-ATOM EXPERIMENTS
GHD and its extensions has been applied to the description of cold-atom experiments that are described by the Lieb-Liniger model with confining potential (159), both in the weakly [38] and in the strongly [39] interacting regimes. A comprehensive review of these works is presented in the recent review [40]. In Fig. VI we show some of the results obtained in Ref. 38: a gas of very weakly interacting 87 Rb atoms is prepared in a thermal equilibrium state in a one-dimensional double-well potential, i.e. like in the example discussed in section I A the transverse confinement is very steep, but the potential along the y-direction in Fig. I A(a) is a double-well rather than harmonic. At time t = 0 the confining potential along the y-direction is changed to a shallow harmonic potential as in Fig. I A(a), and the two clouds of gases start expanding along the y-direction. The density is then measured as a function of time. The main result, and message of this section, is that the GHD prediction for this expansion is found to be in very good agreement with experiment.