Locally quasi-stationary states in noninteracting spin chains

Locally quasi-stationary states (LQSS) were introduced as inhomogeneous generalisations of stationary states in integrable systems. Roughly speaking, LQSSs look like stationary states, but only locally. Despite their key role in hydrodynamic descriptions, an unambiguous definition of LQSSs was not given. By solving the dynamics in inhomogeneous noninteracting spin chains, we identify the set of LQSSs as a subspace that is invariant under time evolution, and we explicitly construct the latter in a generalised XY model. As a by-product, we exhibit an exact generalised hydrodynamic theory (including"quantum corrections").


Introduction
The dynamics of integrable quantum many-body systems prepared in inhomogeneous states have attracted an increasing attention for more than two decades , arguably for their potentiality to elucidate basic questions of quantum transport. The last years in particular have seen a lot of progress in the physical and technical understanding of such dynamics. On the one hand, universal properties have been uncovered within conformal field theories or, more generally, studying the time evolution of states sufficiently close to ground states of critical systems [23][24][25][26][27][28][29][30][31][32][33][34]; on the other hand, a generalised hydrodynamic description (GHD) has been shown to capture the behaviour of integrable systems at large scales [35,36]. GHD was originally based on a conjecture on the expectation value of the currents in stationary states and was tested by studying both entanglement entropies [37][38][39][40][41][42][43][44] and correlation functions [67][68][69][70][71][72][73][74]. Remarkably, the theory has been also verified in one experiment [75]. More recently, analytical proofs of its validity in specific integrable models have appeared [76,77]. First-order GHD looks like a collection of standard continuity equations for particle densities; the big step that delayed its discovery for several years was however the identification of the velocities entering the continuity equations. The main complication that arises in the presence of interactions is indeed the impossibility to define particle velocities independently of the state of the system [78]: the triviality of GHD is only apparent, as the velocities satisfy other equations coupled in a nonlinear way to the GHD ones.
Since its discovery, it has been evident that first-order GHD is a very powerful tool. It has however some limitations, and it was often blamed to be valid only in limits in which the system behaves essentially in a classical way (indeed, does not appear in first-order GHD). In addition, the theory is unable to capture diffusive behaviour in interacting systems, a subject that has been attracting more and more attention [66,67,70,71,[79][80][81]. In order to overcame this weakness, there have been proposals to go beyond the infinite time/low inhomogeneity limit of the original formulation in interacting integrable systems [79,80,82,83], but, at the same time, some issues in noninteracting spin chains were uncovered. Specifically, the attempt of Ref. [84] to develop a higher-order GHD in noninteracting spin chains was unsuccessful, notwithstanding the author made use of the exact solution to the dynamics (which, incidentally, will be proven in this paper). We will show that such attempt has been spoiled by an inappropriate definition of space dependent root densities, a concept that was introduced in GHD without a definition that could have been directly lifted into a nonperturbative level. The problem is that a generic homogeneous quantum state can not be characterised solely by the so-called root densities [86], and, in principle, one has also to define additional fields that take care of the off-diagonal matrix elements of the density matrix. In the inhomogeneous setting, the distinction between root densities and auxiliary fields is not transparent, and, without a proper definition, such fields couple to the root densities in an obscure way. This issue is evident in interacting integrable systems, as the aforementioned additional fields do not even appear in the thermodynamic Bethe Ansatz [86,87], which is arguably the framework of GHD. Despite not necessarily hindering the study of second order GHD, as Refs [79,80,82,83] testify, the lack of an explicit connection with the actual quantum state makes it almost impossible to exhibit an inhomogeneous initial state that could be described by such equations: how can one be sure that the additional fields vanish in the initial state and that their contribution die out in a shorter time scale?
The main goal of this work is to show that some of the problems in going beyond first-order GHD can be traced back to the presence of ambiguities in the definitions of the quantities used to describe the dynamics. We propose a resolution of them in noninteracting spin chains, and this will allow us to make the framework of GHD precise and exact. We will comment in Section 6 on the interacting case.

Ambiguities in GHD
Generalised hydrodynamics was introduced as a large-scale description of time evolution in one dimensional integrable systems with inhomogeneities. Roughly speaking, at large times only the slowest degrees of freedom keep varying and the state becomes locally equivalent to a stationary state. Such a quasi-stationary state was called LQSS (locally quasi-stationary state) [69], and it was shown to be well described by the continuity equations satisfied by the charge densities in the limit of low inhomogeneity. GHD was simultaneously developed in quantum field theories [35] and quantum spin chains [36], but we restrict here to spin chains.
The basic elements of GHD are the additive charges Q = a Q , commuting with the Hamiltonian H = a H , and their corresponding currents J [Q j ] = a J [Q j ]. Here a denotes the lattice spacing (in the next section, its definition will be slightly modified) and the operators O are quasi-localised around site 1 . In the Heisenberg picture, charges and currents are connected by a continuity equation of the form Here we made a (asymmetric) choice for the position of the currents; we will reconsider the possibility of other conventions later. It is clear that we are free to add any constant in the definition of J [Q j ]; we choose the convention of imposing null trace, tr[J [Q j ]] = 0, which corresponds to the vanishing of the currents at infinite temperature. These equations are sufficient to determine the current for given charge density, but the definition of the latter is still ambiguous. Indeed, the charge density can be always redefined as follows where G is a quasi-localised operator acting around site . This redefinition produces a change not only in the current but also in the integrated current Because of this, the choice of G = a G and of its density G is not irrelevant, especially if one is interested in going beyond first-order GHD. Depending on G, indeed, the resulting current can be conserved or not. In the former case, if the expectation values of the charges can be expressed in terms of densities, somehow describing the diagonal part of the density matrix, the expectation values of the currents would be as well. Non-conserved currents would have instead also off-diagonal matrix elements, which, through the continuity equation (1), would generate also off-diagonal contributions in the expectation values of the charges in inhomogeneous states.
Assuming that one can find a conserved operator C[Q j ] such that the following operator is quasilocal, we propose to set G equal to (5) so that the integrated current becomes conserved We denote it by J [Q] to emphasise that it is now just a functional of the charge and not of its particular density: whatever the local charge density Q j is defined, we have Note that the current is still defined up to the density of a charge C[Q j ], which will be chosen in such a way to make the current of the current conserved, and so on and so forth 2 . This procedure defines an invariant subspace of operators, linear combinations of charge densities, that we call locally quasi-conserved operators, LQCO. The density matrix of an LQSS can then be defined as the exponential of an LQCO. Generally, in noninteracting spin chains with local Hamiltonians, the standard choice for the local conservation laws [45,88,89] is such that (5) is quasilocal even for C[Q j ] = 0; thus, it is reasonable to expect that (6) could be enforced. We will see that (6) plays a key role in the inhomogeneous generalisation of the so-called root densities, which were originally introduced to characterise stationary states in integrable models that allow for a thermodynamic Bethe Ansatz description [86,87]. Specifically, it allows us to define root densities that depend on the site and that completely determine the expectation values of the charge densities and vice versa.
We mention that similar issues have been already discussed [77,80]. In particular, in the derivation of a diffusive correction to GHD, refs [79,80] have proposed to choose a PTinvariant gauge, which turns out to be unambiguous once requiring the expectation values of the charge densities to be real in stationary states. Our gauge choice, on the other hand, is not driven by uniqueness: we aim at exploiting the degrees of freedom in the definitions of the charge densities in such a way to decouple the dynamical equations describing the time evolution of the locally quasi-stationary states from the rest.
Finally, we stress that the ambiguities pointed out in this section must be lifted in order to compare analytical predictions against numerical (or experimental) data, where following different conventions could result in artificial discrepancies. Besides, unfortunate conventions could introduce cumbersome large-time corrections that could be otherwise reabsorbed in the definitions of the operators. This is becuase, as we will see, the auxiliary fields characterising the state together with the root densities satisfy dynamical equations different from the GHD ones, so their presence in the initial state affects the large time corrections. For example, Ref. [84] could not capture a universal part of the light-cone dynamics within the proposed higher-order GHD. The problem was not in the universality of the phenomenon, while in having defined space-dependent root densities in an inappropriate way. Specifically, the inhomogeneous generalisation of the root density proposed in Ref. [84] is the quantity that, in the next section, will be called "spuriously local root density". While the latter has the correct homogeneous limit, in the presence of inhomogeneities it is not dynamically decoupled by the remaining degrees of freedom. To overcome this problem and write a dynamical equation written solely in terms of the root density, Ref. [84] projected at every time on the subspace of states that were completely characterised by the spuriously local root density. The resulting equation was not supposed to be exact, it was instead questioned whether it could nevertheless be valid at the next order in the inhomogeneity. The answer was disappointedly negative. In this paper we will not make any approximation and show how to fix this issue.

Outline
Section 2 -Preliminaries -introduces the systems under investigation and a convenient representation of states and operators.
Section 3 -Time evolution of inhomogeneous states -collects the main results. In particular, an invariant subspace made of LQSSs is identified, and the equations of motion governing time evolution within that subspace are exhibited (for homogeneous Hamiltonians, they are also solved) and identified with generalised hydrodynamics. We show that the complementary subspace is invariant as well and exhibit the corresponding dynamical equations.
Section 4 -The two-temperature scenario revisited -details the solution to the dynamics in one of the most studied settings. The predictions are checked against numerical data. Section 5 -Continuum scaling limit -compares the lattice dynamics with the dynamics in the quantum field theory emerging in the limit where the lattice spacing a approaches 0.
Section 6 -Conclusion -is a collection of observations.
Appendix A -Non-equilibrium dynamics in inhomogeneous systems -has a proof of the Moyal dynamical equation in noninteracting spin chains.
Appendix B -Weak inhomogeneous limit: an asymptotic expansion -presents a perturbative derivation of the equations of motion through a formal asymptotic expansion of the exact solution in the limit where a parameter characterising the inhomogeneity approaches infinity (homogeneous limit).

List of the notations
[A, B] ± = AB ± BA denote the anticommutator and the commutator, respectively.
Z, N 0 , 1 2 Z, and 1 2 + Z are the set of the integers, the set of the nonnegative integers, the set of the integers and the half-integers, and the set of the half-integers, respectively.
O, in a bold capital letter, is a linear operator acting on the spin chain.
O, in a calligraphic capital letter, is the matrix associated with the noninteracting operator that is indicated with the same capital letter in bold (cf. (13)), i.e., O.
o, with a hat on a lower case letter, denotes the symbol of the matrix indicated with the same letter, in a calligraphic upper case (cf. (13) and (14)), i.e., O. We also call it the symbol of the operator indicated with the corresponding bold capital letter, i.e., O.
o phys is the physical part ofô, which completely characterizes O.
o ± (e ip ) are the π-periodic and the π-anti-periodic parts ofô with respect to ap.
2κ is the dimension of the representation of the symbol.
a denotes κ times the lattice spacing.
denotes the Moyal star product.
x (sometimes y) and t are the position and the time. They generally appear as subscripts. If the quantity of interest does not depend on x and/or t, the letter is simply not shown.

The Hamiltonian
We consider the nonequilibrium time evolution generated by a noninteracting spin-chain Hamiltonian whose density is almost invariant under a shift by κ sites. This class of Hamiltonians includes: a) the quantum Ising model [90,91] (κ = 1) where σ α act like Pauli matrices on site and like the identity elsewhere, and Z is the set of all the integers; b) the XX model [91,92] (κ = 1) c) the (noninteracting) spin-Peierls model [93] (κ = 2) In its most general form, the Hamiltonian is as follows (the homogeneous case was introduced in Ref. [94]) where N 0 is the set of the nonnegative integers. The coupling constants can be arbitrary, but, in almost the entire paper, we will assume that J ;αβ n;i,j and J ;z i do not depend on . Under the Jordan-Wigner transformation where [A, B] + = AB + BA denotes the anticommutator, the Hamiltonian is mapped into a chain of noninteracting Majorana fermions, as follows Here H is an infinite purely imaginary Hermitian matrix (H ji n ) * = H ij n = −H ji n . We point out that, contrary to standard conventions, we call H a matrix despite its indices not being strictly positive.
If the coupling constants are smooth functions of the position, the Hamiltonian is approximately translationally invariant in any sufficiently small region; this property can be exploited by writing H in a "local" Fourier space Here 1 2 Z is the set consisting of the integers and the half-integers; is the reduced Planck constant and a is κ times the lattice spacing. We point out that κ can be any positive integer, but, in the (quasi) homogeneous case, it is conveniently chosen to be proportional to the number of sites the system is (almost) invariant under a shift of.
The symbol. The (2κ)-by-(2κ) matrixĥ x (e iap / ) represents the Hamiltonian in the Fourier space around position x; we call it 'symbol', by analogy with the definition of the symbol of (block-)Toeplitz and (block-)circulant matrices. It will be clearer later that p plays the role of the conjugate variable of x, so we call it momentum.
We note that, for given x, half of the Fourier components ofĥ x (e iap / ) do not enter the expression. In particular, for (half-)integer x/a, the relevant part of the symbol, which we call h phys x (e iap / ), is a π-(anti-)periodic matrix function of ap / . The remaining Fourier components can be chosen arbitrarily, and hence they will be referred to as "unphysical". Nevertheless, we require the symbol to remain Hermitian and to satisfy [ĥ x (e iap / )] t = −ĥ x (e − iap / ), where t denotes transposition. We also defineĥ ± x (e iap / ) as the extensions of the symbols obtained interpolating either the functions evaluated at x a ∈ Z or at x a ∈ 1 2 + Z, the latter being the set of the half-integers. Thus we havê h phys A possible extension to x a ∈ R is the followinĝ Using these definitions,ĥ x (e iap / ) andĥ (±) x (e iap / ) are entire (matrix-)functions of x. The formal expressions that will be derived in this paper (and that can be applied to any noninteracting spin chain) will be explicitly worked out in a generalised XY model with κ = 1. In the homogeneous case, such model describes a generic one-site shift invariant noninteracting spin-chain Hamiltonian. Its symbol will be parametrised as follows: where the subscripts 'e' and 'o' stand for the even and the odd part respectively (with respect to the momentum p); ε(p) is the dispersion relation; θ(p) = −θ(−p) is the Bogoliubov angle and φ(p) = φ(−p) parametrises the angle of the axis of the Bogoliubov rotation. When we want to emphasize the dependence on the phase, we write θ(p) = Θ(e iap / ). For example, in the quantum Ising model (a) we have ε(p) = 2J 1 + g 2 − 2g cos( ap / ), e iθ(p) = 2J(e i ap / −g) , and φ(p) = 0. If not stated otherwise, all the quantities will be assumed to be smooth functions of the momentum, which generally happens in (quasi)local noncritical Hamiltonians.
A convenient parametrization for inhomogeneous Hamiltonians will be shown later -(72).

The state
The systems considered in this paper are prepared in a gaussian state, i.e., a state where the Wick's theorem can be applied to the correlations of the Majorana fermions (12) 'pf' denoting the Pfaffian. This can be the ground state (if a symmetry is not spontaneously broken) or an excited state of a noninteracting spin-chain Hamiltonian, like (11), as well as a thermal state or a generalised Gibbs ensemble [95][96][97]. Since the dynamics are noninteracting, Wick's theorem holds at any time, and the expectation values of the observables can be written in terms of the correlation matrix Here · t denotes the expectation value at time t. As done for the quadratic operators, we define the physical part of the symbolΓ x,t (e iap ) of the correlation matrix around position x as follows [Γ phys where z p = e iap / . At a given time, we can extend (20) so as to allow for real x using, for example, (16); however, we will not impose (16) at a generic time, being more convenient to exploit such degrees of freedom to simplify the equations of motion -Appendix A. Nevertheless, we will assume that the extension remains an entire function of the position at any time.
For κ = 1 (which is a convenient representation when the state is almost invariant under a shift by one site) the state is completely characterised by a spuriously local root density ρ fake x,t (p), and by a spuriously local auxiliary complex field Ψ fake This representation, which at first glance could look bizarre (for the appearance of the angles θ(p) and φ(p), which are properties of the Hamiltonian and not of the state), allows us in fact to attribute a physical meaning to ρ fake x,t (p) and Ψ fake x,t (p): in the homogeneous limit, ρ fake (p) becomes the density in phase space δxdp of the excitations over the ground state of the quasiparticle excitations (provided that ε(p) ≥ 0) [45], and Ψ fake (p) describes the off-diagonal part of the density matrix. On the other hand, such properties are lost in the presence of inhomogeneities. Indeed Ref. [84] showed that a higher-order hydrodynamics based on such spurious quantities can not be exact. The main aim of this paper is to construct a genuinely local version of these quantities, which we will simply call ρ(p) and Ψ(p). They will be chosen so as to accomplish the goal outlined in 1.1; in particular, they will reduce to their global counterparts in the homogeneous limit, and they will transform as simply as possible under time evolution.
We note that, in the homogeneous limit, the root density and the auxiliary field satisfy the following constraints (the maximal eigenvalue ofΓ(e iap ) in absolute value can not exceed 1) These bounds are not always satisfied in inhomogeneous systems.

Expectation values and local operators
The expectation value of a traceless quadratic operator O can be expressed in terms of its symbol and of the symbol of the correlation matrix as follows These equations are exact; the drawback is that they involve the physical part of one of the two symbols; in Section 3.3 we will exhibit an explicitly symmetrical expression. We can express the symbol of operators (quasi)localised around as follows: where e 1 and o 1 select the even and the odd part with respect to the first argument. The symbol of the corresponding translationally invariant operator iŝ For κ = 1 the auxiliary function will be parametrised as followŝ In terms of the spuriously local root density and spuriously local auxiliary field the expectation value reads as One of the goals of this paper is to give a convenient definition of charge density so that its expectation value has a form similar to (27) with the spuriously local root density and auxiliary field replaced by the analogous genuinely local quantities.

Reduced density matrix
The reduced density matrix (RDM) of a connected subsystem A ≡ [x 1 , x 2 ] in a gaussian state is gaussian. Its correlation matrix is obtained by restricting the correlation matrix of the state to the subsystem. Within the Wigner description that we introduced, it is convenient to embed the RDM into the entire system as follows: which has a correlation matrix that zeros outside A. The symbol of ρ ext A is then given bŷ the distribution in convolution with the symbol approaches a Dirac delta function, so, as expected, the bulk is not affected by the boundaries of the subsystem.
We refer the reader to Ref. [45] for more examples and for an overview of the description in terms of symbols in homogeneous noninteracting spin chains.

Time evolution of inhomogeneous states
It is well known that the dynamics generated by quadratic operators have infinitely many invariant subspaces, indeed it conserves the "number" of the Majorana fermions 3 In this section we show that, in fact, it is reducible even more, and one of the invariant subspaces consists of the locally quasi-stationary states, qualitatively introduced in Ref. [69]. In the following the reduced Planck constant and the lattice spacing a will be written explicitly so that the reader can identify quantum and lattice contributions more easily.

Moyal dynamical equation
The symbol of an inhomogeneous state time evolves according to the following Moyal dynamical equation where z p = e i ap / , and f x (p) g x (p) is the Moyal star product [98] A proof of (31) is reported in Appendix A. Using that the star product is associative, the solution to (31) can be written aŝ where This can be simplified if the Hamiltonian is homogeneous (ĥ x (z p ) is independent of x), as followsΓ The solution to (35) is readily obtained and reads aŝ These equations are the fundamentals of a Wigner description [99] of the dynamics in noninteracting spin-chain models, and they generalise the well-known equations describing free particle systems, in the sense that functions are replaced by matrices of functions (a reader could be more comfortable with referring to the symbol of the correlation matrix as a 'block Wigner function'). They have been formerly announced in Ref. [84]. We stress that (31) and (36) apply to spin chains, despite the variable position x not being discrete.

Invariant subspaces
The dynamics described by (31) are reducible. We identify two invariant subspaces: the inhomogeneous generalisation of stationary states, which, following the terminology of Ref. [69], will be called LQSSs, and a subspace of states that we call ODSs (off-diagonal states).

Locally quasi-stationary states
Let the Hamiltonian be translationally invariant. An LQSS has the following properties: -it reduces to a stationary state in the homogeneous limit; -it remains an LQSS under time evolution.
Imposing these conditions in the κ = 1 representation, we find that the state can be characterised by a Wigner function ρ x (p), which we identify with the localised version of the root density, and it appears in the symbol of the correlation matrix as followŝ with In (37), the arrows on top of the derivatives indicate the direction where they act. Roughly speaking, G parametrises a gauge invariance in the definition of the root density. For example, it includes the arbitrariness in associating a position to a quasi-localised operator. Because of G, the same root density can be used to describe different states, therefore, in order to avoid confusion, it must be fixed. If not stated otherwise, G(z, w) will be set to unity and the notations will be eased, dropping the dependence on G in the expressions, e.g.Γ LQSS Crucially, contrary to standard conventions, the Wigner function that we have defined does not depend only on the state, but also on the Hamiltonian; this is one of the reasons why we prefer to call it 'root density', which is instead a quantity that is always defined for given Hamiltonian and reference state.
Note that only integer and half-integer positions appear in (37). In addition, for G = 1 the physical part of the root density is π periodic at integer sites and π anti-periodic at half-integer sites.
Finally, we remind the reader of a very well known property of Wigner functions: even if ρ x (p) is a physical root density at fixed x (namely, 0 ≤ ρ x (p) ≤ 1 2π ), the LQSS associated with ρ x (p) could be unphysical, as its correlation matrix could have eigenvalues with modulus larger than 1. Thus, not every profile of charges is allowed in an LQSS.

Off-diagonal states
An ODS (off-diagonal state) has the following properties: -it reduces to a state with vanishing diagonal elements in a basis of stationary states in the homogeneous limit; -it remains an ODS under time evolution.
Imposing these conditions in the κ = 1 representation, we find that the state is completely characterised by an auxiliary field Ψ x (p), which appears in the symbol of the correlation matrix as followŝ Also in this case, G parametrises a gauge invariance in the definition of the local auxiliary field, with the same caveats pointed out in the LQSS case. If not stated otherwise, G(z, w) will be set to unity.

Dynamics
Assuming a one-site shift invariant Hamiltonian, each state can be conveniently written in a κ = 1 representation, decomposing it in an LQSS part, characterised by a root density ρ x (p), plus an off-diagonal one, characterised by an auxiliary field Ψ x (p); the symbol of the correlation matrix can be written as follows: By definition, the root density degree of freedom remains decoupled from the auxiliary field one also at finite time; the dynamics are described by the decoupled Moyal equations (43) where, in the second line, we used the oddity of the auxiliary field to write the right hand side as an antisymmetric functional of ε and Ψ. Equations (42) and (43) are arguably the main result of this paper for homogeneous Hamiltonians. In particular, they allow us to identify the (localised version of the) root density, as it appears in the symbol of the correlation matrix in (41), with a Wigner function even when there is a nontrivial Bogoliubov angle θ(p) around a direction characterised by φ(p). This means, for example, that they hold true also in the quantum Ising model. We postpone the proof of (42) and (43) to Section 3.4, where we consider the more general case of an inhomogeneous Hamiltonian. Equations (42) and (43) clearly show that the root density ρ and the auxiliary field Ψ characterise two separate subsets of states that are invariant under time evolution. Incidentally, we note that Ψ x,t (p) Ψ * x,t (p) and Ψ * x,t (p) Ψ x,t (p) satisfy the same equation of the root density. The first orders of the asymptotic expansions of (42) and (43) in the limit of low inhomogeneity are where v(p) = dε(p) dp is the velocity of the quasiparticle excitation with momentum p. We note that the first equation truncated at O(∂ x ) is the noninteracting version of first-order GHD [35,36,69]. It also corresponds to the classical limit → 0, in which the root density satisfies a classical continuity equation. On the other hand, as long as the dispersion relation is not odd, Ψ x,t (p) describes purely quantum contributions. It is worth observing that the equations governing time evolution do not depend explicitly on the lattice spacing.
Equations (42) and (43) are solved by Localisation. The root density and the auxiliary field can be expressed in terms of the spurious quantities appearing in (21) by inverting (41). For the sake of simplicity, we do it explicitly only for the generalised XY model with φ(p) = 0. We find where and Θ(e i ap / ) is the Bogoliubov angle as defined below (17). We stress that the local root density at (half-)integer position is a mix of quantities defined at both integer and half-integer positions.
The first orders of the expansion in the limit of weak inhomogeneity read as In the XX model the Bogoliubov angle can be set to zero, therefore the spuriously local root densities are in fact genuinely local.

Locally quasi-conserved operators
If there are no special symmetries, in noninteracting spin chain models the conservation laws are quadratic forms of Majorana fermions, like the Hamiltonian in (14) [45]. In the Heisenberg picture, the matrices associated with quadratic operators time evolve like the correlation matrix, provided to reverse the sign of the time. This allows us to identify an invariant subspace of operators, which we call LQCOs (locally quasi-conserved operators), consisting of the quadratic forms of fermions with symbols of the form of an LQSS (37) The analogue of the root density is the (single particle) eigenvalue's density in phase space, which can be defined as follows where G parametrises the usual gauge invariance. In the Heisenberg picture, q x (p) time evolves according to the Moyal equation which manifests the closure of the LQCOs under time evolution. By definition, the density of a conserved operator is an LQCO. We do not expect the existence of interacting (i.e., not quadratic) charges, which, in turn, rules out the existence of interacting locally quasi-conserved operators; to the best of our knowledge, however, this has not been proved.
Finally, we remind the reader that the set of charges that are obtained by fixing the parameter κ (in our case we have chosen κ = 1) is not always complete -this was originally pointed out in Ref. [100]. In particular, additional charges appear when the symbol of the Hamiltonian in some representation κ is degenerate for generic momentum, which is possible when the dispersion relation in the κ = 1 representation has symmetries like ε(p) = ±ε(p+π).
Charge densities. A charge density can be defined as a locally quasi-conserved operator, associated with a given site, that is conserved when integrated over the position (i.e., summed over all the sites and multiplied by a). By virtue of (64), it is natural to define its density as the quadratic form of fermions with the following symbol where s = 0, 1 2 is such that q(p + π a ) = (−1) 2s q(p).
Currents. The current associated with a given charge can be easily defined from the Moyal dynamical equation (53), which can be rewritten as We can then use that, independently of the exact definition, the translational invariance of the charges implies the charge density q ( ) x,t (p) to be a function of x − a. Thus we have i.e.
Finally, we can define the current as the LQCO with eigenvalue By integrating over the position a, we find the eigenvalue of the total current j(p) = v(p)q(p).
We note that such definition is different from the one exhibited in [45], where the charge densities were defined in a different way (in particular, they were chosen to be local).

Off-diagonal operators
Similarly to what we did for states, we also identify another class of quadratic operators, which we call "off-diagonal operators", that is invariant under time evolution. They have the In the Heisenberg picture, b x (p) time evolves as follows

Expectation values
Formula (23) can be written in an explicitly symmetrical form: For κ = 1, by decomposing a quasi-localised quadratic operator in its quasi-conserved and off-diagonal parts, we can improve the parametrisation (26) as follows: In this way, formula (27) still holds, provided to replace spuriously local root density and auxiliary field by their genuinely local counterparts.
Locally quasi-conserved operators. The expectation value of an LQCO in a gaussian state with correlation matrix as in (41) takes a rather simple form (in the standard gauge G = 1): It is then convenient to classify the LQCOs depending on how their eigenvalues transform under p → p + π a : Charge densities. Using Ansatz (54), the expectation value of a generic charge density reads as

Inhomogeneous Hamiltonians
In this section we extend the previous analysis to inhomogeneous Hamiltonians. If the symbol of the correlation matrix is parametrised as in (41), it could be convenient to express the symbol of the Hamiltonian as followŝ where w x (p) = −w x (−p) is an odd complex field.
Eq. (66) generates the following dynamics This parametrisation clearly manifests that the two invariant subspaces identified before for homogeneous Hamiltonians, i.e., the LQSS and the ODS, are invariant even when time evolution is generated by a locally quasi-conserved operator (w x (p) = 0). We mention that, if time evolution is generated by an off-diagonal operator (ε x (p) = 0), the dynamics can still be written in a decoupled form with the additional (coupled) boundary conditions Equation (66) is a sensible parametrisation when w x (p) can be treated as a small perturbation with respect to ε x (p), for example assuming w x (p) = −i ∂ x c x (p), with c x (p) a weakly inhomogeneous odd complex field with the dimension of a velocity. Such term could arise from a naive definition of an LQCO -for example, overlooking the derivatives in (51). At the first order in the inhomogeneity we then have where the arguments of the functions (x, t and p) are understood.
In the most generic case, it is instead convenient to express the symbol of the Hamiltonian as follows:ĥ whereΘ x (z p ) is a two-by-two Hermitian matrix that generalises the Bogoliubov angle and satisfiesΘ x (1/z p ) = −Θ t x (z p ). The correlation matrix can then be parametrised aŝ In this way the root density and the auxiliary field satisfy a simple generalisation of (42), (43) Proof of (74) and (75). The equations of motions with the parametrisation (72) can be easily worked out by observing that, for any 2-by-2 unitaryÛ (p), the following holds This allows us to reduce the star products between the symbols in (31) to star products between the functions (ρ, Ψ, ε) characterising the symbols, and finally to get (74) and (75).

Integrals of motion
In the inhomogeneous case we can immediately identify a class of commuting quadratic operators, namely, the LQCOs with eigenvalues of the form where Q is a generic -function, i.e., a funcion in which multiplications are replaced by Moyal star products. This class is not expect to be complete, especially if the dispersion relation has particular symmetries. For example, if ε x (p) = ε x (−p), all the charges of the form (77) are even; this contrasts with the homogeneous case where there are infinitely many odd charges (under spatial reflections).

More general invariant subspaces
We still assume κ = 1 and consider time evolution under a locally quasi-conserved operator.
From (67) we can easily infer that there are also invariant subspaces with both root density and complex field different from zero. Specifically, the states lying on surfaces of the form never leave the surfaces. For example, the following state belongs to the invariant subspace on the surface ρ Finally, we note that, if the Hamiltonian is homogeneous, the constant in (78) can be replaced by any function of the momentum.

The two-temperature scenario revisited
In this section we reconsider the two-temperature scenario, which Ref. [61] studied extensively in the transverse-field Ising chain. There, the initial state was the junction of two states prepared at different temperatures in two semi-infinite open chains. Perhaps counterintuitively, that is not an LQSS, and hence the corrections are not purely hydrodynamic. We would like to choose a locally quasi-stationary initial state with a step root density; this would however be inconsistent with our original hypothesis of smoothness of the symbol. Thus, we replace the step function by a smooth function that is equivalent to a step function for x a ∈ 1 2 Z, with a given convention on sgn(0) ∈ {−1, 0, 1}, as follows ρ x,0 (p) = 1 2π The solution to the dynamics is given by (45) From this expression we can readily identify the entire correction to first-order GHD: which generally approaches zero as t − 1 2 in a rapidly oscillating way (for v(p) = x t , it generally does as t −1 ). We remind the reader that the expectation values of local observables have a further integral over the momentum (cf. (61)), therefore, due to the rapidly oscillating terms, the corrections to the first-order result tend to decay to zero faster, generally as t −1 .
For the sake of simplicity, we write down the symbol of the correlation matrix only for reflection symmetric (ε(p) = ε(−p)) generalised XY Hamiltonians with φ(p) = 0 We point out that (83) is not always physical, as some eigenvalues of the corresponding correlation matrix Γ could be outside the physical interval [−1, 1]. In that case, an LQSS with a step root density does not exist, i.e., that state can not be prepared. In the marginal case the correlation matrix has an eigenvalue exactly equal to ±1, and the step profile is an eigenstate of some charge.

Example (reduction to finite chains)
In this subsection we provide a numerical check of our prediction for the time evolution of LQSSs. In order to obtain numerical data by standard means, we are forced to consider finite chains, so we must adapt our description for systems with a finite number of spins. We propose to embed a finite periodic fermionic chain into an infinite one. The first step is to impose periodicity in space ρ x+La,0 (p) = ρ x,0 (p), Ψ x+La,0 (p) = Ψ x,0 (p). For example, (80) could be modified as follows: where the chain is split in 2g parts with alternate temperatures {β L , β R , β L , . . .}; the parts have the same length if L is divisible by 2g. This is not the end of the story: periodicity in space is not sufficient to describe a finite periodic fermionic chain, indeed in that case the momentum is not a continuous variable as in (84) but it is quantised according to e i Lap = 1. This problem can be circumvented by observing that the Moyal star product (32) between any translationally invariant, 2π-periodic function of ap and the root density (84) produces terms where the momentum is shifted by gnπ L a , with n odd; thus, if g is even, the Moyal star product turns out to preserve the Ramond quantisation condition e i Lap = 1. In conclusion, for even g, time evolution in the finite fermionic chain with periodic boundary conditions can be obtained from our results valid in the thermodynamic limit, provided to impose the Ramond quantisation conditions and to replace integrals over momenta with sums. Specifically, the time dependent root density reads For the sake of example, we report the expression for the magnetisation profile which can be written as follows where θ k = θ( k /a), ε k = −1 ε( k /a), and R is the set of all the independent solutions k to e iLk = 1. We warn the reader that, in order to use (87), the Bogoliubov angle must be (defined as) a smooth odd function of k.  (87)). Time evolution is generated by the Hamiltonian

with periodic boundary conditions on the J-W fermions).
We have checked that, for even g, (87) is exactly equal to the local magnetisation in the state obtained by time evolving the initial correlation matrix, which has the symbol σ z σ y sin (4gx−a)(2n−1)π 2La sin (2n−1)π 2L (e iLk = 1) . (88) Figure 1 shows s z (t) when the system is prepared in a 4-partite LQSS 4 .
We conclude with a curiosity, readily recognisable by inspecting (87): despite the total spin in the z direction not being generally conserved, its expectation value S z = s z is stationary. This phenomenon does not actually depend on the observable, but it is a consequence of the state being locally quasi-stationary.
Asymptotic expansion. It can be instructive to carry out the asymptotic expansion of (85), which will give the solution to the generalised hydrodynamic equation at a given order in the inhomogeneity (cf. (44)). We will assume L divisible by 2. Using the parity of the First-order hydrodynamics predicts ρ (1) x,t (p) = ρ x−v(p)t,0 (p), that is to say The solution to third-order hydrodynamics is ρ An attentive reader might have noted that the expansion in the strength of the inhomogeneity is an expansion about p of the argument of the dispersion relations that are multiplied by the time in (89). Appendix B shows that this is indeed general. We can then use this prescription in (87) (with the sum restricted to {1, . . . , L 2 } and multiplied by 2, as in (89)) and check the goodness of higher-order hydrodynamics. Fig. 2 shows that the error of the approximation decreases as the order of GHD is increased. In addition, while the light-cones are clearly visible in the left panel, they are imperceptible in the right panel, showing that third-order hydrodynamics describes the behaviour around the light cones, as conjectured in Ref. [84].

Continuum scaling limit
So far we have considered a spin-chain system at fixed lattice spacing a, providing also perturbative expansions valid when the typical length of the inhomogeneity is much larger than a. In this section we adopt the alternative point of view of keeping the typical length of the inhomogeneity fixed while taking the limit a → 0. The coupling constants can then be chosen to scale with a in such a way to allow for a quantum field theory description of the low-energy degrees of freedom. This is a standard problem, but since our representation in terms of symbols could be unfamiliar to the reader, we describe the procedure in detail.
A preliminary step is to rescale the Majorana operators. We define the continuum Majorana fields as ψ i (x) = lim a→0 ψ i (x; a), where and we used the sinc interpolation to extend the definition of the operators to R; in this way the operators are entire functions of the position even when a is finite (but they satisfy a more complicated algebra). In the following we use the vector notation ψ(x; a) = ψ 1 (x; a) ψ 2 (x; a) . The Hamiltonian can be written as x a , y a ∈ 1 2 Z π a − π a dp 2π e 2i(x−y)p ψ † (2x; a)ĥ x+y (e iap / ) ψ(2y; a) .
By expanding the Majorana fields around x + y, repeated integrations by parts result in We choose the coupling constants in such a way that, in the limit a → 0, the symbol of the Hamiltonian approaches infinity for all momenta but a finite number of them. We then expand the symbol around such points. Physically, the infinities correspond to infinite excitation energies, and the momenta at which the symbol is not divergent describe the lowest (finite) excitations. By definition, the symbol is a 2π-periodic function of ap , therefore the expansion in the momentum is also an expansion in the lattice spacing a. The divergency of some coupling constants will compensate the presence of powers of a in the terms of the expansion. Since generally we consider Hamiltonians with a finite number of coupling constants, only a finite number of the terms of the expansion will remain nonzero in the limit a → 0. In other words, in the scaling limit the symbol becomes a polynomial in the momentum. This is sufficient for the QFT Hamiltonian to be local (in (94), each power of p becomes a derivative with respect to (x − y)).
As a prototypical example, we consider the quantum XY model, whose Hamiltonian's symbol is given bŷ Let us determine in which scaling limits the symbol remains finite for ap ≈ 0. We havê Essentially, there are three options for simplifying the divergencies: where, on the right hand side, we showed the corresponding QFT symbols. Case (97a) is the Ising QFT (see, e.g., [101]); case (97b) describes the multicritical point γ = 0, h = 1; case (97c) describes the limit of infinite anisotropy. In case (97c), the symbol is the direct sum of two terms because, in the scaling limit, there is also another momentum at which the symbol is finite, namely, ap = π; the last term describes that contribution (with the momentum shifted by π). We note that, in (97c), in order to write the corresponding Hamiltonian, one should redefine the Majorana field in such a way to distinguish even and odd sites. In the following, we will restrict to case (97a) with homogeneous couplings, namely Time evolution from a locally quasi-stationary state is then described by the Moyal equation with ε QF T = c m 2 c 2 + p 2 . We aim at comparing this result with the continuum limit of the lattice time evolution. We assume the initial state to be an LQSS. The lattice dynamics are solved by (45); there, the Hamiltonian enters through the quantity where we highlighted the dependence on the lattice spacing. The QFT is the result of the expansion about the zeros of the dispersion relation, but (100) has additional zeros. In our specific case, we can immediately identify a q = 0, a q = 2π, and a p = π, but there can also be further zeros. For example, for γ > √ 2 there are closed curves centred at ( a p, a q) = (π, 0) and ( a p, a q) = (0, 2π) with ω lat (p, q; a) = 0. For the sake of simplicity we assume γ = 1 (quantum Ising model), for which we find (if not stated otherwise, aq 2 ∈ (−π, π) and ap ∈ (−π, π)) lim a→0 ω lat (p, q; a) = So far we have only considered the limit a → 0 in the terms dependent on the Hamiltonian. It is however clear that, if we aim at capturing the continuum scaling limit by the QFT describing the low-energy physics, also the initial state should be described by the QFT. Specifically, excitations that, in the limit a → 0, have infinite energy, should not be present in the state, that is to say lim a→0 ρ x,0 (p) = 0 ap ≈ 0 .
In addition, every nontrivial root density in the lattice for which ∃ lim x→±∞ ρ x (p) = ρ ± (p) will approach a step function in the limit a → 0; let us then assume straight away that ρ x,0 (t) is like in (80). This allows us to simplify the lattice solution in the limit a → 0 valid for ap ≈ 0, otherwise ρ x,t (p) = 0. As expected, the final result is the time evolution in the QFT, which follows equation (99). Importantly, in order to recover the QFT result, we had to assume the initial state to be describable by the QFT. In other words, as testified by (101), in the continuum scaling limit time evolution is not sufficient to extinguish highly energetic degrees of freedom (i.e., in the quantum Ising model, excitations with ap ≈ 0).

Conclusion
In this paper we have studied the fundamentals of generalised hydrodynamics in noninteracting spin chains. GHD was originally developed in a perturbative framework, as the asymptotic solution of time evolution in the limit of large time [35,36] or low inhomogeneity [56,57] in integrable systems. Ref. [69] made a first attempt to lift the theory into a non-perturbative level by conjecturing the existence of so-called "locally quasi-stationary states", which are states completely characterised by a local version of the root densities of the thermodynamic Bethe Ansatz [87]. This suggestion was however not exploited, and the main efforts were rather put in determining the next orders of perturbation theory. There have been indeed proposals to go beyond generalised hydrodynamics, even in the presence of interactions [79,83]. The uncertainty of what should be captured by generalised hydrodynamics has however made it problematic even the determination of the next orders of perturbation theory in noninteracting spin chains without U (1) symmetry [84].
In order to resolve this issue in noninteracting spin chains, we have adopted here the alternative perspective of identifying the invariant subspace made of the locally quasi-stationary states. In one-site shift invariant models with U (1) symmetry, like the XX model, this is a straightforward, intuitive step; this paper provides a solution in the more complicated situation when the number of Jordan-Wigner fermions is not conserved. Once the subspace of locally quasi-stationary states is identified, GHD becomes the theory describing time evolution exactly, indeed the time evolving state is completely determined by the continuity equations satisfied by the charges. This observation has allowed for a non-perturbative derivation of generalised hydrodynamics (cf. Section 3.2.3).
We would like to emphasise that, even if first-order GHD becomes exact only at the Euler scale, the theory does not assume the absence of correlations. Generalised hydrodynamics is a mapping of the initial state into the state at nonzero time; if there are quantum correlations in the initial state, they will be "transported" by GHD. At the first order, the mapping does not introduce new quantum correlations, but the state will still be quantum correlated, even at a finite time. Ref. [74] provides an explicit example of this phenomenon, as it proposes a framework to compute quantum correlations within first-order GHD.
We have also found other invariant subspaces that "cut" the Hilbert space in a different way -Section 3.5. While the existence of so many invariant subspaces could seem extraordinary, Ref. [102] has actually shown that time evolution occurs in a tiny part of the Hilbert space even in generic systems; we then wonder whether approximately invariant subspaces could be defined also in non-integrable models.
Interacting integrable systems. Our findings are restricted to noninteracting spin chains, but we think that a similar change of perspective could be useful also in the presence of interactions. Specifically, our calculations are based on the possibility to define a complete set o quasilocal charge densities in such a way to make all their currents conserved. Such charge densities are then the building blocks for the invariant subspace of locally quasi-stationary states. This very preliminary step is already questionable in interacting integrable systems. For example, in the gapped XXZ spin-1 2 chain, there is evidence that there is a single (quasi-)local charge that is odd under spin flip: the total spin in the direction of the anisotropy. This is enough to exclude the existence of a subspace consisting of LQSSs with a nontrivial profile of the sign of the magnetisation. The results of Ref. [54] go in this direction, indeed the authors had to introduce an auxiliary parameter to describe the time evolution of states without a fixed sign of the local magnetisation. Nevertheless, we do not exclude, and, in fact, we believe, that even in the gapped XXZ model there are two main invariant subspaces consisting of LQSSs with a fixed sign of the local magnetisation. Such invariant subspaces could be investigated within Bethe anstaz, quite independently of the knowledge of the equations of motion describing time evolution within the subspaces, which is instead the program of GHD. As a matter of fact, candidates for LQSSs in interacting integrable systems already exist: for example, they could be stationary states of inhomogeneous integrable models (obtained by introducing inhomogeneities in the spectral parameters of the R-matrix -see, e.g., Ref. [103]). In addition, it is worth observing that, assuming all the quantities used to describe the system to be defined in the appropriate way, all the equations that we obtained could have been guessed by recognising that products in the homogenous case become Moyal star products in the inhomogeneous one. This simple rule could be specific to noninteracting systems, but it is still reasonable to expect that also in interacting integrable models the TBA equations should be modified in some analogous way. More specifically, in the light of our identification of the root density with a Wigner function, the generalisation to the interacting case can be seen as the conception of a phase-space representation of the Bethe Ansatz structure.
Finally, we have observed that the subspace of LQSSs is not only invariant under the effect of a homogeneous Hamiltonian, but also if the Hamiltonian is replaced by a locally quasi-conserved operator (cf. (67)): the algebra of the locally quasi-conserved operators is closed. We wonder whether this property holds true also in the presence of interactions.

A Non-equilibrium dynamics in inhomogeneous systems
In order to ease the notations, in this appendix the lattice spacing a and the reduced Planck constant are set to unity.
In noninteracting models, the correlation matrix (19) satisfies a reduced version of the von Neumann equation (cf. (30)) where H was defined in (13). In homogeneous systems, this can be reduced further to an analogous equation satisfied by the symbolsΓ(e ik ) andĥ(e ik ) (the homogeneous version of (14) and (20)) This equation is extremely powerful: the time evolution of a quantum many-body system is recast into the time evolution of a (2κ)-by-(2κ) matrix, which depends on a parameter, the By repeated integrations by parts, the coefficients of the space derivatives can be rewritten as derivatives with respect to the momentum The sum over j forces q to be either equal to p or to p + π; using the transformation properties of the symbols under a shift of the momentum by π, (15), we then find Putting all the pieces together we get (cf. (108)) This equation describes only the time evolution of the physical part of the correlation matrix. Indeed, strictly speaking, it holds only at the values of x such that p = (−1) 2x . We fix the time evolution of the unphysical part by asking for this equation to hold at any position. Since, in this way, x and p become independent, we can sum the two equations corresponding to p = ±1, obtaining This is the most important result of this section, being the generalization of (105) to the inhomogeneous case. It's a Wigner description of the dynamics, and is valid also in noninteracting systems without a U (1) symmetry (we remind the reader that the standard setting where the so-called Wigner function comes into play is the time evolution of a noninteracting state with a fixed number of particles). The symbol of the correlation matrix could be interpreted as a Wigner function with an internal degree of freedom, but we prefer to identify the Wigner function(s) with the localized version of the root density(ies), as will be discussed in the next section.

B Weak inhomogeneous limit: an asymptotic expansion
In order to ease the notations, in this appendix the lattice spacing a and the reduced Planck constant are set to unity.
If time evolution is generated by a time independent noninteracting Hamiltonian, (116) is solved by (36), which, in an XY-like model with φ(p) = 0, can be written aŝ Here we rescaled position and time back to origin, so that χ has disappeared from the equation. We stress however that the n-th term of the sum is suppressed by a factor χ with respect to the (n − 1)-th term under a rescaling of space and time by χ; the truncation of the series at n = N approximates the correlation matrix up to O(1/χ N +1 ) corrections. We identify two different contributions: s = s: time evolution causes only the mixing of information from different parts of the chain; we will refer to this as the hydrodynamic part of time evolution; s = s: time evolution is characterised by rapidly oscillatory terms; we will refer to this as the purely quantum part of time evolution.
Definition. A locally quasi-stationary state is a state with a null purely quantum part.
Definition. An off-diagonal state is a state with a null hydrodynamic part.
Since in this appendix we only aim at a perturbative solution to the dynamics, we consider the following expansion:Γ x,t (z) . By truncating the number of equations to m = N , the error is O(1/χ N +1 ). A single odd complex field Ψ x (p) completely characterises an off-diagonal state; the first terms of the expansion of its symbol are the following We have defined Ψ x (p) in such a way that, at the leading order, it corresponds to the bare field Ψ bare x (p) as defined in (21).
whereΓ (n) x,t (e ip ) are written as in (131). Up to O(1/ξ 3 ) corrections, the complex field time evolves as follows: